diff --git a/Docs/source/SUNDIALS_Tutorial.rst b/Docs/source/SUNDIALS_Tutorial.rst index 1f8ca4e3..2af438ec 100644 --- a/Docs/source/SUNDIALS_Tutorial.rst +++ b/Docs/source/SUNDIALS_Tutorial.rst @@ -22,7 +22,19 @@ user-supplied reaction and diffusion coefficients: The inputs file contains a template for MRI approaches, where the diffusion process can be treated as a "fast" partition relative to the reaction process. -Please see each respective inputs file or + +The third example code at ``amrex-tutorials/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon`` +solves the advection-diffusion equation, where :math:`A` and :math:`D` are +user-supplied advection and diffusion coefficients: + +.. math:: \frac{\partial\phi}{\partial t} = D \nabla^2\phi - A \nabla\phi. + +The inputs file contains a template for IMEX approaches, where diffusion is +treated implicitly with a SUNDIALS preconditioner that uses AMReX MLMG. The +example can also use MLMG's HYPRE bottom solver interface when AMReX is built +with HYPRE support. + +Please see the inputs file or `AMReX User Guide:Time Integration`_ for more details. .. _`AMReX User Guide:Time Integration`: https://amrex-codes.github.io/amrex/docs_html/TimeIntegration_Chapter.html# diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/CMakeLists.txt b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/CMakeLists.txt new file mode 100644 index 00000000..73863d2b --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/CMakeLists.txt @@ -0,0 +1,15 @@ +if (AMReX_SPACEDIM EQUAL 1 OR NOT AMReX_LINEAR_SOLVERS) + return() +endif () + +# List of source files +set(_sources main.cpp myfunc.H) +list(TRANSFORM _sources PREPEND "Source/") + +# List of input files +file( GLOB_RECURSE _input_files LIST_DIRECTORIES false Exec/input* ) + +setup_tutorial(_sources _input_files) + +unset( _sources ) +unset( _input_files ) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/GNUmakefile b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/GNUmakefile new file mode 100644 index 00000000..106c912b --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/GNUmakefile @@ -0,0 +1,23 @@ +# AMREX_HOME defines the directory in which we will find all the AMReX code. +AMREX_HOME ?= ../../../../../amrex + +DEBUG = FALSE +USE_MPI = TRUE +USE_OMP = FALSE +COMP = gnu +DIM = 2 +USE_SUNDIALS = TRUE +USE_HYPRE ?= FALSE +SUNDIALS_HOME ?= ../../../../../sundials/instdir + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ../Source/Make.package +VPATH_LOCATIONS += ../Source +INCLUDE_LOCATIONS += ../Source + +Pdirs := Base Boundary LinearSolvers/MLMG +Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) +include $(Ppack) + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials new file mode 100644 index 00000000..767ba43f --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials @@ -0,0 +1,11 @@ +Refer to https://amrex-codes.github.io/amrex/docs_html/TimeIntegration_Chapter.html +for SUNDIALS installation, build, and usage instructions. + +Make sure that SUNDIALS_HOME in the GNUmakefile points to the installation directory. +This tutorial also requires AMReX to be built with MLMG linear solver support. +The default inputs use a SUNDIALS IMEX-RK method with GMRES and an AMReX MLMG +preconditioner for the implicit diffusion solve. + +Set USE_HYPRE=TRUE when building and mlmg.use_hypre=1 in inputs_sundials to use +MLMG's HYPRE bottom solver interface. The HYPRE interface and options namespace +are controlled by mlmg.hypre_interface and mlmg.hypre_options_namespace. diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/inputs_sundials new file mode 100644 index 00000000..0ff1e019 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/inputs_sundials @@ -0,0 +1,57 @@ +n_cell = 256 +max_grid_size = 128 + +dt = 1.e-5 +nsteps = 100 +plot_int = 100 + +advCoeffx = 100. +advCoeffy = 100. + +diffCoeffx = 1. +diffCoeffy = 1. + +# Use adaptive time stepping (multi-stage integrators only!) and set integrator relative and absolute tolerances +# adapt_dt = true +# reltol = 1.0e-4 +# abstol = 1.0e-9 + +# INTEGRATION +## integration.type can take on the following values: +## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator +## 1 or "RungeKutta" => Native AMReX Explicit Runge Kutta controlled by integration.rk.type +## 2 or "SUNDIALS" => SUNDIALS backend controlled by integration.sundials.type +#integration.type = ForwardEuler +#integration.type = RungeKutta +integration.type = SUNDIALS + +# Set the SUNDIALS method type: +# ERK = Explicit Runge-Kutta method +# DIRK = Diagonally Implicit Runge-Kutta method +# IMEX-RK = Implicit-Explicit Additive Runge-Kutta method +# EX-MRI = Explicit Multirate Infatesimal method +# IM-MRI = Implicit Multirate Infatesimal method +# IMEX-MRI = Implicit-Explicit Multirate Infatesimal method +# +# Optionally select a specific SUNDIALS method by name, see the SUNDIALS +# documentation for the supported method names +integration.sundials.type = IMEX-RK +integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2 +integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2 +integration.sundials.linear_solver = GMRES +integration.sundials.linear_solver_preconditioning = LEFT +integration.sundials.max_linear_iters = 20 + +# MLMG settings for the SUNDIALS preconditioner solve +mlmg.max_iter = 100 +mlmg.max_fmg_iter = 0 +mlmg.verbose = 0 +mlmg.bottom_verbose = 0 +mlmg.use_hypre = 0 +mlmg.hypre_interface = 3 # 1 = structed, 2 = semi_structed, 3 = ij +mlmg.hypre_options_namespace = hypre +mlmg.reltol = 1.e-10 +mlmg.abstol = 0.0 + +# HYPRE settings used when mlmg.use_hypre = 1 and mlmg.hypre_interface = 3 +hypre.hypre_solver = BoomerAMG diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/Make.package b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/Make.package new file mode 100644 index 00000000..fa194a35 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += main.cpp +CEXE_headers += myfunc.H diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp new file mode 100644 index 00000000..2d9693b9 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp @@ -0,0 +1,514 @@ +#include +#include +#include +#include +#include +#include +#ifdef AMREX_USE_HYPRE +#include +#endif + +#include "myfunc.H" + +#include +#include +#include +#include +#include + +using namespace amrex; + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc,argv); + + main_main(); + + amrex::Finalize(); + return 0; +} + +void main_main () +{ + + if (AMREX_SPACEDIM != 2) { + amrex::Abort("Only 2D supported; recompile with DIM=2"); + } + + // ********************************** + // SIMULATION PARAMETERS + + // number of cells on each side of the domain + int n_cell; + + // size of each box (or grid) + int max_grid_size; + + // total steps in simulation + int nsteps; + + // how often to write a plotfile + int plot_int; + + // time step + Real dt; + + // use adaptive time step (dt used to set output times) + bool adapt_dt = false; + + // adaptive time step relative and absolute tolerances + Real reltol = 1.0e-4; + Real abstol = 1.0e-9; + + // Advection and Diffusion Coefficients + Real advCoeffx = 1.0; + Real advCoeffy = 1.0; + Real diffCoeffx = 1.0; + Real diffCoeffy = 1.0; + + // MLMG settings for the SUNDIALS preconditioner solve + int mlmg_max_iter = 100; + int mlmg_max_fmg_iter = 0; + int mlmg_verbose = 0; + int mlmg_bottom_verbose = 0; + bool mlmg_use_hypre = false; + int mlmg_hypre_interface = 3; + std::string mlmg_hypre_options_namespace = "hypre"; + Real mlmg_reltol = 1.e-10; + Real mlmg_abstol = 0.0; + + // inputs parameters + { + // ParmParse is way of reading inputs from the inputs file + // pp.get means we require the inputs file to have it + // pp.query means we optionally need the inputs file to have it - but we must supply a default here + ParmParse pp; + + // We need to get n_cell from the inputs file - this is the number of cells on each side of + // a square (or cubic) domain. + pp.get("n_cell",n_cell); + + // The domain is broken into boxes of size max_grid_size + pp.get("max_grid_size",max_grid_size); + + // Default nsteps to 10, allow us to set it to something else in the inputs file + nsteps = 10; + pp.query("nsteps",nsteps); + + // Default plot_int to -1, allow us to set it to something else in the inputs file + // If plot_int < 0 then no plot files will be written + plot_int = -1; + pp.query("plot_int",plot_int); + + // time step + pp.get("dt",dt); + + // use adaptive step sizes + pp.query("adapt_dt",adapt_dt); + + // adaptive step tolerances + pp.query("reltol",reltol); + pp.query("abstol",abstol); + + pp.query("advCoeffx",advCoeffx); + pp.query("advCoeffy",advCoeffy); + pp.query("diffCoeffx",diffCoeffx); + pp.query("diffCoeffy",diffCoeffy); + + ParmParse pp_mlmg("mlmg"); + pp_mlmg.query("max_iter",mlmg_max_iter); + pp_mlmg.query("max_fmg_iter",mlmg_max_fmg_iter); + pp_mlmg.query("verbose",mlmg_verbose); + pp_mlmg.query("bottom_verbose",mlmg_bottom_verbose); + pp_mlmg.query("use_hypre",mlmg_use_hypre); + pp_mlmg.query("hypre_interface",mlmg_hypre_interface); + pp_mlmg.query("hypre_options_namespace",mlmg_hypre_options_namespace); + pp_mlmg.query("reltol",mlmg_reltol); + pp_mlmg.query("abstol",mlmg_abstol); + +#ifndef AMREX_USE_HYPRE + if (mlmg_use_hypre) { + amrex::Abort("mlmg.use_hypre requires AMReX to be built with HYPRE support"); + } +#endif + } + + // ********************************** + // SIMULATION SETUP + + // AMREX_D_DECL means "do the first X of these, where X is the dimensionality of the simulation" + IntVect dom_lo(AMREX_D_DECL( 0, 0, 0)); + IntVect dom_hi(AMREX_D_DECL(n_cell-1, n_cell-1, n_cell-1)); + + // Make a single box that is the entire domain + Box domain(dom_lo, dom_hi); + + // ba will contain a list of boxes that cover the domain + // Initialize the boxarray "ba" from the single box "domain" + BoxArray ba(domain); + + // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction + ba.maxSize(max_grid_size); + + // This defines the physical box, [0,1] in each direction. + RealBox real_box({AMREX_D_DECL(-1.,-1.,-1.)}, + {AMREX_D_DECL( 1., 1., 1.)}); + + // periodic in all direction + Array is_periodic{AMREX_D_DECL(1,1,1)}; + + // geom contains information such as the physical domain size, + // number of points in the domain, and periodicity + // This defines a Geometry object + Geometry geom(domain, real_box, CoordSys::cartesian, is_periodic); + + // extract dx from the geometry object + GpuArray dx = geom.CellSizeArray(); + GpuArray prob_lo = geom.ProbLoArray(); + GpuArray prob_hi = geom.ProbHiArray(); + + // Nghost = number of ghost cells for each array + int Nghost = 1; + + // Ncomp = number of components for each array + int Ncomp = 1; + + // How Boxes are distrubuted among MPI processes + DistributionMapping dm(ba); + + // allocate phi MultiFab + MultiFab phi(ba, dm, Ncomp, Nghost); + + // time = starting time in the simulation + Real time = 0.0; + + // ********************************** + // INITIALIZE DATA + + InitializeData(phi,dx,prob_lo,prob_hi,time,advCoeffx,advCoeffy); + + // Write a plotfile of the initial data if plot_int > 0 + if (plot_int > 0) + { + int step = 0; + const std::string& pltfile = amrex::Concatenate("plt",step,5); + WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, 0); + } + + auto rhs_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) { + + // fill periodic ghost cells + S_data.FillBoundary(geom.periodicity()); + + S_rhs.setVal(0.); + + ComputeDiffusion(S_rhs, S_data, diffCoeffx, diffCoeffy, dx); + ComputeAdvection(S_rhs, S_data, advCoeffx, advCoeffy, dx); + }; + + auto rhs_im_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) { + + // fill periodic ghost cells + S_data.FillBoundary(geom.periodicity()); + + S_rhs.setVal(0.); + + ComputeDiffusion(S_rhs, S_data, diffCoeffx, diffCoeffy, dx); + }; + + auto rhs_ex_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) { + + // fill periodic ghost cells + S_data.FillBoundary(geom.periodicity()); + + S_rhs.setVal(0.); + + ComputeAdvection(S_rhs, S_data, advCoeffx, advCoeffy, dx); + }; + + struct MLMGPreconditioner { + std::unique_ptr linop; + std::unique_ptr solver; + std::array face_bcoef; + }; + + std::unique_ptr mlmg_preconditioner; + Real mlmg_gamma = std::numeric_limits::quiet_NaN(); + + auto gamma_has_changed = [&] (Real gamma) + { + if (!std::isfinite(mlmg_gamma)) { + return true; + } + const Real scale = std::max({Real(1.0), std::abs(gamma), std::abs(mlmg_gamma)}); + return std::abs(gamma - mlmg_gamma) > + Real(10.0) * std::numeric_limits::epsilon() * scale; + }; + + auto update_mlmg_gamma = [&](MLMGPreconditioner& preconditioner, Real gamma) + { + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const Real diff_coeff = (idim == 0) ? diffCoeffx : ((idim == 1) ? diffCoeffy : 0.0); + preconditioner.face_bcoef[idim].setVal(gamma * diff_coeff); + } + preconditioner.linop->setBCoeffs(0, amrex::GetArrOfConstPtrs(preconditioner.face_bcoef)); + mlmg_gamma = gamma; + }; + + auto build_mlmg_preconditioner = [&](Real gamma) + { + auto preconditioner = std::make_unique(); + + LPInfo info; + preconditioner->linop = std::make_unique( + Vector{geom}, Vector{ba}, Vector{dm}, info); + preconditioner->linop->setMaxOrder(2); + preconditioner->linop->setDomainBC( + {AMREX_D_DECL(LinOpBCType::Periodic, LinOpBCType::Periodic, LinOpBCType::Periodic)}, + {AMREX_D_DECL(LinOpBCType::Periodic, LinOpBCType::Periodic, LinOpBCType::Periodic)}); + preconditioner->linop->setLevelBC(0, nullptr); + preconditioner->linop->setScalars(1.0, 1.0); + preconditioner->linop->setACoeffs(0, 1.0); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const BoxArray& face_ba = amrex::convert(ba, IntVect::TheDimensionVector(idim)); + preconditioner->face_bcoef[idim].define(face_ba, dm, 1, 0); + } + update_mlmg_gamma(*preconditioner, gamma); + + preconditioner->solver = std::make_unique(*preconditioner->linop); + preconditioner->solver->setMaxIter(mlmg_max_iter); + preconditioner->solver->setMaxFmgIter(mlmg_max_fmg_iter); + preconditioner->solver->setVerbose(mlmg_verbose); + preconditioner->solver->setBottomVerbose(mlmg_bottom_verbose); +#ifdef AMREX_USE_HYPRE + if (mlmg_use_hypre) { + Hypre::Interface hypre_interface = Hypre::Interface::ij; + if (mlmg_hypre_interface == 1) { + hypre_interface = Hypre::Interface::structed; + } else if (mlmg_hypre_interface == 2) { + hypre_interface = Hypre::Interface::semi_structed; + } else if (mlmg_hypre_interface == 3) { + hypre_interface = Hypre::Interface::ij; + } else { + amrex::Abort("mlmg.hypre_interface must be 1 (structed), 2 (semi_structed), or 3 (ij)"); + } + preconditioner->solver->setBottomSolver(MLMG::BottomSolver::hypre); + preconditioner->solver->setHypreInterface(hypre_interface); + preconditioner->solver->setHypreOptionsNamespace(mlmg_hypre_options_namespace); + } +#endif + + mlmg_preconditioner = std::move(preconditioner); + mlmg_gamma = gamma; + }; + + auto precond_setup = [&](MultiFab& /* S_data */, MultiFab& /* S_rhs */, const Real /* time */, + bool jok, bool& jcur, const Real gamma) + { + // The implicit diffusion Jacobian is constant for this problem, so a + // SUNDIALS refresh only needs to update the gamma-scaled operator. + if (mlmg_preconditioner == nullptr) { + build_mlmg_preconditioner(gamma); + jcur = true; + return; + } + + if (!jok) { + if (gamma_has_changed(gamma)) { + update_mlmg_gamma(*mlmg_preconditioner, gamma); + } + jcur = true; + } else if (gamma_has_changed(gamma)) { + update_mlmg_gamma(*mlmg_preconditioner, gamma); + } + }; + + auto precond_solve = [&](MultiFab& S_soln, MultiFab& S_rhs, MultiFab& /* S_data */, + MultiFab& /* S_state_rhs */, const Real /* time */, + const Real /* gamma */, const Real /* delta */, + int /* lr */) + { + AMREX_ALWAYS_ASSERT(mlmg_preconditioner != nullptr); + S_soln.setVal(0.0); + mlmg_preconditioner->solver->solve({&S_soln}, {&S_rhs}, mlmg_reltol, mlmg_abstol); + }; + + TimeIntegrator integrator(phi, time); + integrator.set_rhs(rhs_function); + integrator.set_imex_rhs(rhs_im_function, rhs_ex_function); + integrator.set_preconditioner(precond_setup, precond_solve); + + if (adapt_dt) { + integrator.set_adaptive_step(); + integrator.set_tolerances(reltol, abstol); + } else { + integrator.set_time_step(dt); + } + + Real evolution_start_time = ParallelDescriptor::second(); + + for (int step = 1; step <= nsteps; ++step) + { + // Set time to evolve to + time += dt; + + Real step_start_time = ParallelDescriptor::second(); + + // Advance to output time + integrator.evolve(phi, time); + + Real step_stop_time = ParallelDescriptor::second() - step_start_time; + ParallelDescriptor::ReduceRealMax(step_stop_time); + + // Tell the I/O Processor to write out which step we're doing + amrex::Print() << "Advanced step " << step << " in " << step_stop_time << " seconds; dt = " << dt << " time = " << time << "\n"; + + // Write a plotfile of the current data (plot_int was defined in the inputs file) + if (plot_int > 0 && step%plot_int == 0) + { + const std::string& pltfile = amrex::Concatenate("plt",step,5); + WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, step); + } + } + + Real evolution_stop_time = ParallelDescriptor::second() - evolution_start_time; + ParallelDescriptor::ReduceRealMax(evolution_stop_time); + amrex::Print() << "Total evolution time = " << evolution_stop_time << " seconds\n"; + + // exact solution + BoxArray ba_exact(domain); + DistributionMapping dm_exact(ba_exact); + MultiFab phi_exact(ba_exact, dm_exact, Ncomp, n_cell); + InitializeData(phi_exact,dx,prob_lo,prob_hi,time,advCoeffx,advCoeffy); + phi_exact.SumBoundary(geom.periodicity()); + + { + const std::string& pltfile = amrex::Concatenate("exact",nsteps,5); + WriteSingleLevelPlotfile(pltfile, phi_exact, {"phi"}, geom, time, nsteps); + } + + MultiFab phi_exact_dist(ba,dm,Ncomp,0); + phi_exact_dist.ParallelCopy(phi_exact,0,0,1); + + MultiFab::Subtract(phi_exact_dist,phi,0,0,1,0); + + { + const std::string& pltfile = amrex::Concatenate("diff",nsteps,5); + WriteSingleLevelPlotfile(pltfile, phi_exact_dist, {"phi"}, geom, time, nsteps); + } + + Real error = phi_exact_dist.norm1(0,geom.periodicity()); + amrex::Print() << "L1 error = " << error << std::endl; +} + +void InitializeData(MultiFab& phi, + const GpuArray dx, + const GpuArray prob_lo, + const GpuArray prob_hi, + const Real& time, + const Real& Ax, + const Real& Ay) { + + int ng = phi.nGrow(); + + GpuArray L; + for (int d=0; d& phi_array = phi.array(mfi); + + Real sigma = 0.1; + Real a = 1.0/(sigma*sqrt(2*M_PI)); + Real b = -0.5/(sigma*sigma); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + Real x = prob_lo[0] + (((Real) i) + 0.5) * dx[0]; + Real y = prob_lo[1] + (((Real) j) + 0.5) * dx[1]; + Real r = (x-Ax*time) * (x-Ax*time) + (y-Ay*time) * (y-Ay*time); + phi_array(i,j,k) = a * std::exp(b * r); + }); + } + +} + +void ComputeDiffusion(MultiFab& S_rhs, + MultiFab& S_data, + const Real& Dx, + const Real& Dy, + const GpuArray dx) { + + for ( MFIter mfi(S_data,TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + const Box& bx = mfi.tilebox(); + + const Array4& phi_array = S_data.array(mfi); + const Array4< Real>& rhs_array = S_rhs.array(mfi); + + // fill the right-hand-side for phi + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + rhs_array(i,j,k) += Dx * ( (phi_array(i+1,j,k) - 2.*phi_array(i,j,k) + phi_array(i-1,j,k)) / (dx[0]*dx[0]) ) + + Dy * ( (phi_array(i,j+1,k) - 2.*phi_array(i,j,k) + phi_array(i,j-1,k)) / (dx[1]*dx[1]) ); + }); + } +} + + +void ComputeAdvection(MultiFab& S_rhs, + MultiFab& S_data, + const Real& Ax, + const Real& Ay, + const GpuArray dx) { + + Real dxInv = 1.0 / dx[0]; + Real dyInv = 1.0 / dx[1]; + Real sideCoeffx = Ax * dxInv; + Real sideCoeffy = Ay * dyInv; + + for (MFIter mfi(S_data,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + const Array4& phi_array = S_data.array(mfi); + const Array4< Real>& rhs_array = S_rhs.array(mfi); + + // x-direction + if (Ax > 0) + { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + rhs_array(i,j,k) -= sideCoeffx * (phi_array(i,j,k) - phi_array(i-1,j,k)); + }); + } + else + { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + rhs_array(i,j,k) -= sideCoeffx * (phi_array(i+1,j,k) - phi_array(i,j,k)); + }); + } + + // y-direction + if (Ay > 0) + { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + rhs_array(i,j,k) -= sideCoeffy * (phi_array(i,j,k) - phi_array(i,j-1,k)); + }); + } + else + { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + rhs_array(i,j,k) -= sideCoeffy * (phi_array(i,j+1,k) - phi_array(i,j,k)); + }); + } + } +} diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/myfunc.H b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/myfunc.H new file mode 100644 index 00000000..0879e55d --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/myfunc.H @@ -0,0 +1,30 @@ +#ifndef MYFUNC_H_ +#define MYFUNC_H_ + +#include + +using namespace amrex; + +void main_main (); + +void InitializeData(MultiFab& phi, + const GpuArray dx, + const GpuArray prob_lo, + const GpuArray prob_hi, + const Real& time, + const Real& Ax, + const Real& Ay); + +void ComputeDiffusion(MultiFab& S_rhs, + MultiFab& S_data, + const Real& Dx, + const Real& Dy, + const GpuArray dx); + +void ComputeAdvection(MultiFab& S_rhs, + MultiFab& S_data, + const Real& Ax, + const Real& Ay, + const GpuArray dx); + +#endif