diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/CMakeLists.txt b/ExampleCodes/SUNDIALS/Reaction-Diffusion/CMakeLists.txt new file mode 100644 index 00000000..71eb16ae --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/CMakeLists.txt @@ -0,0 +1,15 @@ +if (AMReX_SPACEDIM EQUAL 1) + 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/Reaction-Diffusion/Exec/GNUmakefile b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/GNUmakefile new file mode 100644 index 00000000..7947b0b7 --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/GNUmakefile @@ -0,0 +1,51 @@ +# 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 = 3 +USE_RPATH = TRUE +USE_SUNDIALS = TRUE + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ../Source/Make.package +VPATH_LOCATIONS += ../Source +INCLUDE_LOCATIONS += ../Source + +ifeq ($(USE_SUNDIALS),TRUE) + ifeq ($(USE_CUDA),TRUE) + SUNDIALS_ROOT ?= $(TOP)../../../../../sundials/instdir_cuda + else + SUNDIALS_ROOT ?= $(TOP)../../../../../sundials/instdir + endif + ifeq ($(NERSC_HOST),perlmutter) + SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib64 + else + SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib + endif + + USE_CVODE_LIBS ?= TRUE + USE_ARKODE_LIBS ?= TRUE + + DEFINES += -DAMREX_USE_SUNDIALS + INCLUDE_LOCATIONS += $(SUNDIALS_ROOT)/include + LIBRARY_LOCATIONS += $(SUNDIALS_LIB_DIR) + + LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_cvode + LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_arkode + LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_nvecmanyvector + LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_core + + ifeq ($(USE_CUDA),TRUE) + LIBRARIES += -L$(SUNDIALS_LIB_DIR) -lsundials_nveccuda + endif + +endif + +include $(AMREX_HOME)/Src/Base/Make.package + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules + diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/README_sundials b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/README_sundials new file mode 100644 index 00000000..c0f315aa --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/README_sundials @@ -0,0 +1,85 @@ +SUNDIALS installation guide: +https://computing.llnl.gov/projects/sundials/faq#inst + +Installation + +# Download sundials-x.y.z.tar.gz file for SUNDIALS and extract it at the same level as amrex. +# Update 6/11/24 - v7.0.0. tested +https://computing.llnl.gov/projects/sundials/sundials-software + +>> tar -xzvf sundials-x.y.z.tar.gz # where x.y.z is the version of sundials + +# at the same level that amrex is cloned, do: + +>> mkdir sundials +>> cd sundials + +###################### +HOST BUILD +###################### + +>> mkdir instdir +>> mkdir builddir +>> cd builddir + +>> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir/examples -DENABLE_MPI=ON ../../sundials-x.y.z + +###################### +NVIDIA/CUDA BUILD +###################### + +# Navigate back to the 'sundials' directory and do: + +>> mkdir instdir_cuda +>> mkdir builddir_cuda +>> cd builddir_cuda + +>> cmake -DCMAKE_INSTALL_PREFIX=/pathto/sundials/instdir_cuda -DEXAMPLES_INSTALL_PATH=/pathto/sundials/instdir_cuda/examples -DENABLE_CUDA=ON -DENABLE_MPI=ON ../../sundials-x.y.z + +###################### + +>> make -j4 +>> make install + +# in your .bashrc or preferred configuration file, add the following (and then "source ~/.bashrc") + +# If you have a CPU build: + +>> export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/pathto/sundials/instdir/lib/ + +# If you have a NVIDIA/CUDA build: + +>> export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/pathto/sundials/instdir_cuda/lib/ + +# now you are ready to compile amrex-tutorials/ExampleCodes/SUNDIALS/Exec with: + +>> make -j4 # optional to have 'USE_CUDA=TRUE' as well + +# in your inputs file, you will need to have: + +# 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.strategy +integration.type = + +## Native AMReX Explicit Runge-Kutta parameters +# +## integration.rk.type can take the following values: +### 0 = User-specified Butcher Tableau +### 1 = Forward Euler +### 2 = Trapezoid Method +### 3 = SSPRK3 Method +### 4 = RK4 Method +integration.rk.type = 3 + +## If using the SUNDIALS Submodule, then +## compile with USE_SUNDIALS=TRUE or AMReX_SUNDIALS=ON and +## set strategy here: +# +## integration.sundials.strategy can take the following values: +### ERK = ERKStep from ARKode with SSPRK3 Method +### MRI = MRIStep from ARKode with Explict Trapezoid Method +### MRITEST = MRIStep from ARKode modified to use no-op inner f0 +integration.sundials.strategy = ERK diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_forward_euler b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_forward_euler new file mode 100644 index 00000000..39b09830 --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_forward_euler @@ -0,0 +1,14 @@ +n_cell = 32 +max_grid_size = 16 + +nsteps = 1000 +plot_int = 100 + +dt = 1.e-5 + +# 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.strategy +integration.type = ForwardEuler diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_rk3 b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_rk3 new file mode 100644 index 00000000..6eb60137 --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_rk3 @@ -0,0 +1,24 @@ +n_cell = 32 +max_grid_size = 16 + +nsteps = 1000 +plot_int = 100 + +dt = 1.e-5 + +# 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.strategy +integration.type = RungeKutta + +## Native AMReX Explicit Runge-Kutta parameters +# +## integration.rk.type can take the following values: +### 0 = User-specified Butcher Tableau +### 1 = Forward Euler +### 2 = Trapezoid Method +### 3 = SSPRK3 Method +### 4 = RK4 Method +integration.rk.type = 3 diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials new file mode 100644 index 00000000..bbdf2c69 --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials @@ -0,0 +1,44 @@ +n_cell = 32 +max_grid_size = 16 + +nsteps = 1000 +plot_int = 100 + +dt = 1.e-5 + +# Use adaptive time stepping 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 +# +# If using the SUNDIALS Submodule, then compile with USE_SUNDIALS=TRUE or +# AMReX_SUNDIALS=ON +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 + +# Use the SUNDIALS default method for the chosen type (fixed or adaptive step sizes) +integration.sundials.type = ERK + +# Use forward Euler (fixed step sizes only) +# integration.sundials.type = ERK +# integration.sundials.method = ARKODE_FORWARD_EULER_1_1 + +# Use backward Euler (fixed step sizes only) +# integration.sundials.type = DIRK +# integration.sundials.method = ARKODE_BACKWARD_EULER_1_1 diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/Make.package b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/Make.package new file mode 100644 index 00000000..fa194a35 --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += main.cpp +CEXE_headers += myfunc.H diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp new file mode 100644 index 00000000..4eab09d4 --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp @@ -0,0 +1,233 @@ +#include +#include +#include + +#include "myfunc.H" + +using namespace amrex; + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc,argv); + + main_main(); + + amrex::Finalize(); + return 0; +} + +void main_main () +{ + + // ********************************** + // 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; + + // 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); + + } + + // ********************************** + // SIMULATION SETUP + + // make BoxArray and Geometry + // ba will contain a list of boxes that cover the domain + // geom contains information such as the physical domain size, + // number of points in the domain, and periodicity + BoxArray ba; + Geometry geom; + + // 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); + + // Initialize the boxarray "ba" from the single box "domain" + ba.define(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( 0., 0., 0.)}, + {AMREX_D_DECL( 1., 1., 1.)}); + + // periodic in all direction + Array is_periodic{AMREX_D_DECL(1,1,1)}; + + // This defines a Geometry object + geom.define(domain, real_box, CoordSys::cartesian, is_periodic); + + // extract dx from the geometry object + GpuArray dx = geom.CellSizeArray(); + + // 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); + + // we allocate two phi multifabs; one will store the old state, the other the new. + MultiFab phi(ba, dm, Ncomp, Nghost); + + // time = starting time in the simulation + Real time = 0.0; + + // ********************************** + // INITIALIZE DATA + + // loop over boxes + for (MFIter mfi(phi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + + const Array4& phi_array = phi.array(mfi); + + // set phi = 1 + e^(-(r-0.5)^2) + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + Real x = (i+0.5) * dx[0]; + Real y = (j+0.5) * dx[1]; +#if (AMREX_SPACEDIM == 2) + Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5))/0.01; +#elif (AMREX_SPACEDIM == 3) + Real z= (k+0.5) * dx[2]; + Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/0.01; +#endif + phi_array(i,j,k) = 1. + std::exp(-rsquared); + }); + } + + // 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 pre_rhs_function = [&](MultiFab& S_data, const Real /* time */) { + // fill periodic ghost cells + S_data.FillBoundary(geom.periodicity()); + }; + + auto rhs_function = [&](MultiFab& S_rhs, const MultiFab& S_data, + const Real /* time */) { + + // loop over boxes + auto& phi_data = S_data; + auto& phi_rhs = S_rhs; + + //Reaction and Diffusion Coefficients + Real reaction_coef = 0.0; + Real diffusion_coef = 1.0; + ParmParse pp; + pp.query("reaction_coef",reaction_coef); + pp.query("diffusion_coef",diffusion_coef); + + for ( MFIter mfi(phi_data); mfi.isValid(); ++mfi ) + { + const Box& bx = mfi.validbox(); + + const Array4& phi_array = phi_data.array(mfi); + const Array4& phi_rhs_array = phi_rhs.array(mfi); + + // fill the right-hand-side for phi + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_rhs_array(i,j,k) = diffusion_coef*( (phi_array(i+1,j,k) - 2.*phi_array(i,j,k) + phi_array(i-1,j,k)) / (dx[0]*dx[0]) + +(phi_array(i,j+1,k) - 2.*phi_array(i,j,k) + phi_array(i,j-1,k)) / (dx[1]*dx[1]) +#if (AMREX_SPACEDIM == 3) + +(phi_array(i,j,k+1) - 2.*phi_array(i,j,k) + phi_array(i,j,k-1)) / (dx[2]*dx[2]) ) +#endif + + reaction_coef*phi_array(i,j,k); + }); + } + }; + + TimeIntegrator integrator(phi, time); + integrator.set_pre_rhs_action(pre_rhs_function); + integrator.set_rhs(rhs_function); + if (adapt_dt) { + integrator.set_adaptive_step(); + integrator.set_tolerances(reltol, abstol); + } else { + integrator.set_time_step(dt); + } + + for (int step = 1; step <= nsteps; ++step) + { + // Set time to evolve to + time += dt; + + // Advance to output time + integrator.evolve(phi, time); + + // Tell the I/O Processor to write out which step we're doing + amrex::Print() << "Advanced step " << step << "\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); + } + } +} diff --git a/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/myfunc.H b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/myfunc.H new file mode 100644 index 00000000..0f828c4e --- /dev/null +++ b/ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/myfunc.H @@ -0,0 +1,6 @@ +#ifndef MYFUNC_H_ +#define MYFUNC_H_ + +void main_main (); + +#endif