From dcd04a4851680457af2382a8bc8632d9b59ddda9 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Fri, 14 Nov 2025 14:04:55 -0800 Subject: [PATCH 01/15] initial commit, copying reaction-diffusion example as starting template --- .../AdvDiff-HyprePrecon/CMakeLists.txt | 15 + .../AdvDiff-HyprePrecon/Exec/GNUmakefile | 20 ++ .../AdvDiff-HyprePrecon/Exec/README_sundials | 4 + .../Exec/inputs_sundials_mri | 68 +++++ .../AdvDiff-HyprePrecon/Source/Make.package | 2 + .../AdvDiff-HyprePrecon/Source/main.cpp | 285 ++++++++++++++++++ .../AdvDiff-HyprePrecon/Source/myfunc.H | 6 + 7 files changed, 400 insertions(+) create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt new file mode 100644 index 00000000..71eb16ae --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/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/AdvDiff-HyprePrecon/Exec/GNUmakefile b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile new file mode 100644 index 00000000..c8feaf99 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile @@ -0,0 +1,20 @@ +# 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_SUNDIALS = TRUE +SUNDIALS_HOME ?= ../../../../../sundials/instdir + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ../Source/Make.package +VPATH_LOCATIONS += ../Source +INCLUDE_LOCATIONS += ../Source + +include $(AMREX_HOME)/Src/Base/Make.package + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials new file mode 100644 index 00000000..3b6830a6 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials @@ -0,0 +1,4 @@ +Refer to https://amrex-codes.github.io/amrex/docs_html/TimeIntegration_Chapter.html +for SNUDIALS installation, build, and usage instructions. + +Make sure that SUNDIALS_HOME in the GNUmakefile points to the installation directory. diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri new file mode 100644 index 00000000..92d09735 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri @@ -0,0 +1,68 @@ +## +# Theoretical forward euler explicit time step constraint for each process is based +# on the diffusion and reaction coefficients. +# For diffusion term +# dt < dx^2 /(2*DIM*diffusion_coef) = (1/32)^2/(2*3*1) = 1.63e-4 +# For reaction term +# dt <= 2/reaction_coef + +n_cell = 32 +max_grid_size = 16 + +nsteps = 12 # 60 # 600 +plot_int = 2 # 10 # 100 + +# With integration.sundials.type = ERK, dt = 1.6e4 is stable, dt = 1.7e-4 is unstable +# With integration.sundials.type = EX-MRI, dt = 1.7e-4 with fast_dt_ratio = 0.1 is stable +# With integration.sundials.type = EX-MRI, dt = 1.7e-3 with fast_dt_ratio = 0.1 is stable +# With integration.sundials.type = EX-MRI, dt = 8.5e-3 with fast_dt_ratio = 0.1 FAILS +# With integration.sundials.type = EX-MRI, dt = 8.5e-3 with fast_dt_ratio = 0.02 is stable +dt = 8.5e-3 + +# To replicate heat equation +diffusion_coef = 1.0 +reaction_coef = 1.e-4 + +# MRI parameters +use_MRI = true +fast_dt_ratio = 0.02 + + +# 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 = EX-MRI +integration.sundials.fast_type = ERK + +## *** Select a specific SUNDIALS ERK method *** +#integration.sundials.method = ARKODE_BOGACKI_SHAMPINE_4_2_3 +# +## *** Select a specific SUNDIALS ImEx method *** +#integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2 +#integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2 + +# *** Select a specific SUNDIALS MRI method *** +integration.sundials.method = ARKODE_MIS_KW3 +integration.sundials.fast_method = ARKODE_KNOTH_WOLKE_3_3 + diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package new file mode 100644 index 00000000..fa194a35 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += main.cpp +CEXE_headers += myfunc.H diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp new file mode 100644 index 00000000..0f78e831 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -0,0 +1,285 @@ +#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; + + // use MRI + bool use_MRI = false; + Real fast_dt_ratio = 0.1; + + // adaptive time step relative and absolute tolerances + Real reltol = 1.0e-4; + Real abstol = 1.0e-9; + + // Reaction and Diffusion Coefficients + Real reaction_coef = 0.0; + Real diffusion_coef = 1.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); + + // use MRI + pp.query("use_MRI",use_MRI); + pp.query("fast_dt_ratio",fast_dt_ratio); + + // adaptive step tolerances + pp.query("reltol",reltol); + pp.query("abstol",abstol); + + pp.query("reaction_coef",reaction_coef); + pp.query("diffusion_coef",diffusion_coef); + } + + // ********************************** + // 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); + + // allocate phi MultiFab + 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 rhs_fast_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) { + + // fill periodic ghost cells + S_data.FillBoundary(geom.periodicity()); + + // loop over boxes + auto& phi_data = S_data; + auto& phi_rhs = S_rhs; + + 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 + ); + }); + } + }; + + auto rhs_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) { + + // fill periodic ghost cells + S_data.FillBoundary(geom.periodicity()); + + // loop over boxes + auto& phi_data = S_data; + auto& phi_rhs = S_rhs; + + 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) + { + if (use_MRI) { + phi_rhs_array(i,j,k) = -1.*reaction_coef*phi_array(i,j,k); + } else { + 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 + ) - 1.*reaction_coef*phi_array(i,j,k); + } + }); + } + }; + + TimeIntegrator integrator(phi, time); + integrator.set_rhs(rhs_function); + if (use_MRI) { + integrator.set_fast_rhs(rhs_fast_function); + integrator.set_fast_time_step(fast_dt_ratio*dt); + } + + if (adapt_dt && !use_MRI) { + 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"; +} diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H new file mode 100644 index 00000000..0f828c4e --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H @@ -0,0 +1,6 @@ +#ifndef MYFUNC_H_ +#define MYFUNC_H_ + +void main_main (); + +#endif From 8ae3f7258dfa947b4e56f062ba05f41d9a4708a6 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Fri, 14 Nov 2025 15:52:11 -0800 Subject: [PATCH 02/15] updates --- .../AdvDiff-HyprePrecon/Exec/GNUmakefile | 2 +- .../AdvDiff-HyprePrecon/Source/main.cpp | 96 ++++++------------- 2 files changed, 28 insertions(+), 70 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile index c8feaf99..657fa127 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile @@ -5,7 +5,7 @@ DEBUG = FALSE USE_MPI = TRUE USE_OMP = FALSE COMP = gnu -DIM = 3 +DIM = 2 USE_SUNDIALS = TRUE SUNDIALS_HOME ?= ../../../../../sundials/instdir diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index 0f78e831..261f9390 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -19,6 +19,10 @@ int main (int argc, char* argv[]) void main_main () { + if (AMREX_SPACEDIM != 2) { + amrex::Abort("Only 2D supported; recompile with DIM=2"); + } + // ********************************** // SIMULATION PARAMETERS @@ -40,17 +44,15 @@ void main_main () // use adaptive time step (dt used to set output times) bool adapt_dt = false; - // use MRI - bool use_MRI = false; - Real fast_dt_ratio = 0.1; - // adaptive time step relative and absolute tolerances Real reltol = 1.0e-4; Real abstol = 1.0e-9; - // Reaction and Diffusion Coefficients - Real reaction_coef = 0.0; - Real diffusion_coef = 1.0; + // Advection and Diffusion Coefficients + Real advCoeffx = 1.0; + Real advCoeffy = 1.0; + Real diffCoeffx = 1.0; + Real diffCoeffy = 1.0; // inputs parameters { @@ -81,16 +83,14 @@ void main_main () // use adaptive step sizes pp.query("adapt_dt",adapt_dt); - // use MRI - pp.query("use_MRI",use_MRI); - pp.query("fast_dt_ratio",fast_dt_ratio); - // adaptive step tolerances pp.query("reltol",reltol); pp.query("abstol",abstol); - pp.query("reaction_coef",reaction_coef); - pp.query("diffusion_coef",diffusion_coef); + pp.query("advCoeffx",advCoeffx); + pp.query("advCoeffx",advCoeffy); + pp.query("diffCoeffx",diffCoeffx); + pp.query("diffCoeffx",diffCoeffy); } // ********************************** @@ -117,7 +117,7 @@ void main_main () ba.maxSize(max_grid_size); // This defines the physical box, [0,1] in each direction. - RealBox real_box({AMREX_D_DECL( 0., 0., 0.)}, + RealBox real_box({AMREX_D_DECL(-1.,-1.,-1.)}, {AMREX_D_DECL( 1., 1., 1.)}); // periodic in all direction @@ -128,6 +128,8 @@ void main_main () // 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; @@ -151,21 +153,18 @@ void main_main () 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) + 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 = (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); + Real y = prob_lo[1] + (((Real) j) + 0.5) * dx[1]; + Real x = prob_lo[0] + (((Real) i) + 0.5) * dx[0]; + Real r = x * x + y * y; + phi_array(i,j,k) = a * std::exp(b * r); }); } @@ -177,35 +176,6 @@ void main_main () WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, 0); } - auto rhs_fast_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) { - - // fill periodic ghost cells - S_data.FillBoundary(geom.periodicity()); - - // loop over boxes - auto& phi_data = S_data; - auto& phi_rhs = S_rhs; - - 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 - ); - }); - } - }; - auto rhs_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) { // fill periodic ghost cells @@ -225,28 +195,16 @@ void main_main () // fill the right-hand-side for phi amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (use_MRI) { - phi_rhs_array(i,j,k) = -1.*reaction_coef*phi_array(i,j,k); - } else { - 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 - ) - 1.*reaction_coef*phi_array(i,j,k); - } + phi_rhs_array(i,j,k) = diffCoeffx * ( (phi_array(i+1,j,k) - 2.*phi_array(i,j,k) + phi_array(i-1,j,k)) / (dx[0]*dx[0]) ) + + diffCoeffy * ( (phi_array(i,j+1,k) - 2.*phi_array(i,j,k) + phi_array(i,j-1,k)) / (dx[1]*dx[1]) ); }); } }; TimeIntegrator integrator(phi, time); integrator.set_rhs(rhs_function); - if (use_MRI) { - integrator.set_fast_rhs(rhs_fast_function); - integrator.set_fast_time_step(fast_dt_ratio*dt); - } - if (adapt_dt && !use_MRI) { + if (adapt_dt) { integrator.set_adaptive_step(); integrator.set_tolerances(reltol, abstol); } else { From 6d57dea067c334aa1ad7942e5ccfaa637db0d2e3 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Sat, 15 Nov 2025 07:09:10 -0800 Subject: [PATCH 03/15] code refactor --- .../AdvDiff-HyprePrecon/Exec/inputs_sundials | 38 +++++++++++ .../Exec/inputs_sundials_mri | 68 ------------------- .../AdvDiff-HyprePrecon/Source/main.cpp | 44 +++++++----- .../AdvDiff-HyprePrecon/Source/myfunc.H | 10 +++ 4 files changed, 74 insertions(+), 86 deletions(-) create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials new file mode 100644 index 00000000..d1c8f921 --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials @@ -0,0 +1,38 @@ +n_cell = 32 +max_grid_size = 16 + +nsteps = 10 +plot_int = 10 + +dt = 1.e-4 + +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 = ERK + + diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri deleted file mode 100644 index 92d09735..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials_mri +++ /dev/null @@ -1,68 +0,0 @@ -## -# Theoretical forward euler explicit time step constraint for each process is based -# on the diffusion and reaction coefficients. -# For diffusion term -# dt < dx^2 /(2*DIM*diffusion_coef) = (1/32)^2/(2*3*1) = 1.63e-4 -# For reaction term -# dt <= 2/reaction_coef - -n_cell = 32 -max_grid_size = 16 - -nsteps = 12 # 60 # 600 -plot_int = 2 # 10 # 100 - -# With integration.sundials.type = ERK, dt = 1.6e4 is stable, dt = 1.7e-4 is unstable -# With integration.sundials.type = EX-MRI, dt = 1.7e-4 with fast_dt_ratio = 0.1 is stable -# With integration.sundials.type = EX-MRI, dt = 1.7e-3 with fast_dt_ratio = 0.1 is stable -# With integration.sundials.type = EX-MRI, dt = 8.5e-3 with fast_dt_ratio = 0.1 FAILS -# With integration.sundials.type = EX-MRI, dt = 8.5e-3 with fast_dt_ratio = 0.02 is stable -dt = 8.5e-3 - -# To replicate heat equation -diffusion_coef = 1.0 -reaction_coef = 1.e-4 - -# MRI parameters -use_MRI = true -fast_dt_ratio = 0.02 - - -# 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 = EX-MRI -integration.sundials.fast_type = ERK - -## *** Select a specific SUNDIALS ERK method *** -#integration.sundials.method = ARKODE_BOGACKI_SHAMPINE_4_2_3 -# -## *** Select a specific SUNDIALS ImEx method *** -#integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2 -#integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2 - -# *** Select a specific SUNDIALS MRI method *** -integration.sundials.method = ARKODE_MIS_KW3 -integration.sundials.fast_method = ARKODE_KNOTH_WOLKE_3_3 - diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index 261f9390..9c7f8766 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -181,24 +181,7 @@ void main_main () // fill periodic ghost cells S_data.FillBoundary(geom.periodicity()); - // loop over boxes - auto& phi_data = S_data; - auto& phi_rhs = S_rhs; - - 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) = diffCoeffx * ( (phi_array(i+1,j,k) - 2.*phi_array(i,j,k) + phi_array(i-1,j,k)) / (dx[0]*dx[0]) ) - + diffCoeffy * ( (phi_array(i,j+1,k) - 2.*phi_array(i,j,k) + phi_array(i,j-1,k)) / (dx[1]*dx[1]) ); - }); - } + ComputeDiffusion(S_rhs, S_data, diffCoeffx, diffCoeffy, dx); }; TimeIntegrator integrator(phi, time); @@ -241,3 +224,28 @@ void main_main () ParallelDescriptor::ReduceRealMax(evolution_stop_time); amrex::Print() << "Total evolution time = " << evolution_stop_time << " seconds\n"; } + +void ComputeDiffusion(MultiFab& S_rhs, + MultiFab& S_data, + const Real& Dx, + const Real& Dy, + const GpuArray dx) { + + auto& phi_data = S_data; + auto& phi_rhs = S_rhs; + + 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) = 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]) ); + }); + } +} diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H index 0f828c4e..d4a13855 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H @@ -1,6 +1,16 @@ #ifndef MYFUNC_H_ #define MYFUNC_H_ +#include + +using namespace amrex; + void main_main (); +void ComputeDiffusion(MultiFab& S_rhs, + MultiFab& S_data, + const Real& Dx, + const Real& Dy, + const GpuArray dx); + #endif From 7bc2b266fa708802befe9d88f92e639a84e73208 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Sat, 15 Nov 2025 07:28:08 -0800 Subject: [PATCH 04/15] advDiff WIP --- .../AdvDiff-HyprePrecon/Exec/inputs_sundials | 9 ++- .../AdvDiff-HyprePrecon/Source/main.cpp | 75 ++++++++++++++++--- .../AdvDiff-HyprePrecon/Source/myfunc.H | 6 ++ 3 files changed, 77 insertions(+), 13 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials index d1c8f921..ebfc683b 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials @@ -1,13 +1,16 @@ n_cell = 32 max_grid_size = 16 -nsteps = 10 +nsteps = 100 plot_int = 10 dt = 1.e-4 -diffCoeffx = 1. -diffCoeffy = 1. +advCoeffx = 1. +advCoeffy = 1. + +diffCoeffx = 0.1 +diffCoeffy = 0.1 # Use adaptive time stepping (multi-stage integrators only!) and set integrator relative and absolute tolerances # adapt_dt = true diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index 9c7f8766..a4790c9a 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -180,8 +180,9 @@ void main_main () // fill periodic ghost cells S_data.FillBoundary(geom.periodicity()); - + ComputeDiffusion(S_rhs, S_data, diffCoeffx, diffCoeffy, dx); + ComputeAdvection(S_rhs, S_data, advCoeffx, advCoeffy, dx); }; TimeIntegrator integrator(phi, time); @@ -231,21 +232,75 @@ void ComputeDiffusion(MultiFab& S_rhs, const Real& Dy, const GpuArray dx) { - auto& phi_data = S_data; - auto& phi_rhs = S_rhs; - - for ( MFIter mfi(phi_data); mfi.isValid(); ++mfi ) + for ( MFIter mfi(S_data,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - const Box& bx = mfi.validbox(); + const Box& bx = mfi.tilebox(); - const Array4& phi_array = phi_data.array(mfi); - const Array4& phi_rhs_array = phi_rhs.array(mfi); + 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) { - phi_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]) ); + 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-HyprePrecon/Source/myfunc.H b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H index d4a13855..1addc9e3 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H @@ -13,4 +13,10 @@ void ComputeDiffusion(MultiFab& S_rhs, const Real& Dy, const GpuArray dx); +void ComputeAdvection(MultiFab& S_rhs, + MultiFab& S_data, + const Real& Ax, + const Real& Ay, + const GpuArray dx); + #endif From 1acc339e86d3ed9434a098ad83105bf472cc603a Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Sat, 15 Nov 2025 07:28:23 -0800 Subject: [PATCH 05/15] tweak parameteres --- .../SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials index ebfc683b..95299e51 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials @@ -6,11 +6,11 @@ plot_int = 10 dt = 1.e-4 -advCoeffx = 1. -advCoeffy = 1. +advCoeffx = 10. +advCoeffy = 10. -diffCoeffx = 0.1 -diffCoeffy = 0.1 +diffCoeffx = 1. +diffCoeffy = 1. # Use adaptive time stepping (multi-stage integrators only!) and set integrator relative and absolute tolerances # adapt_dt = true From 0130ccd4f2507725638700972bde152001ff87c1 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 17 Nov 2025 09:15:34 -0800 Subject: [PATCH 06/15] fix bugs --- .../AdvDiff-HyprePrecon/Exec/inputs_sundials | 8 ++++---- .../SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp | 14 ++++++-------- 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials index 95299e51..ca8dc102 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials @@ -6,8 +6,8 @@ plot_int = 10 dt = 1.e-4 -advCoeffx = 10. -advCoeffy = 10. +advCoeffx = 100. +advCoeffy = 100. diffCoeffx = 1. diffCoeffy = 1. @@ -22,9 +22,9 @@ diffCoeffy = 1. ## 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 = ForwardEuler #integration.type = RungeKutta -integration.type = SUNDIALS +#integration.type = SUNDIALS # Set the SUNDIALS method type: # ERK = Explicit Runge-Kutta method diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index a4790c9a..ce176509 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -180,6 +180,8 @@ void main_main () // 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); @@ -272,16 +274,14 @@ void ComputeAdvection(MultiFab& S_rhs, { 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)); + 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)); + rhs_array(i,j,k) -= sideCoeffx * (phi_array(i+1,j,k) - phi_array(i,j,k)); }); } @@ -290,16 +290,14 @@ void ComputeAdvection(MultiFab& S_rhs, { 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)); + 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)); + rhs_array(i,j,k) -= sideCoeffy * (phi_array(i,j+1,k) - phi_array(i,j,k)); }); } } From 7e15f698ff431242e9dd597b5b1e14b8a28e2fa9 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 17 Nov 2025 09:27:59 -0800 Subject: [PATCH 07/15] add imex --- .../AdvDiff-HyprePrecon/Source/main.cpp | 64 +++++++++++++------ .../AdvDiff-HyprePrecon/Source/myfunc.H | 4 ++ 2 files changed, 50 insertions(+), 18 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index ce176509..9f348365 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -149,24 +149,7 @@ void main_main () // ********************************** // INITIALIZE DATA - // loop over boxes - for (MFIter mfi(phi); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.validbox(); - const Array4& 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 y = prob_lo[1] + (((Real) j) + 0.5) * dx[1]; - Real x = prob_lo[0] + (((Real) i) + 0.5) * dx[0]; - Real r = x * x + y * y; - phi_array(i,j,k) = a * std::exp(b * r); - }); - } + InitializeData(phi,dx,prob_lo); // Write a plotfile of the initial data if plot_int > 0 if (plot_int > 0) @@ -187,8 +170,29 @@ void main_main () 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); + }; + TimeIntegrator integrator(phi, time); integrator.set_rhs(rhs_function); + integrator.set_imex_rhs(rhs_im_function, rhs_ex_function); if (adapt_dt) { integrator.set_adaptive_step(); @@ -228,6 +232,30 @@ void main_main () amrex::Print() << "Total evolution time = " << evolution_stop_time << " seconds\n"; } +void InitializeData(MultiFab& phi, + const GpuArray dx, + const GpuArray prob_lo) { + + // loop over boxes + for (MFIter mfi(phi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + const Array4& 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 y = prob_lo[1] + (((Real) j) + 0.5) * dx[1]; + Real x = prob_lo[0] + (((Real) i) + 0.5) * dx[0]; + Real r = x * x + y * y; + phi_array(i,j,k) = a * std::exp(b * r); + }); + } +} + void ComputeDiffusion(MultiFab& S_rhs, MultiFab& S_data, const Real& Dx, diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H index 1addc9e3..7f21b4ed 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H @@ -7,6 +7,10 @@ using namespace amrex; void main_main (); +void InitializeData(MultiFab& phi, + const GpuArray dx, + const GpuArray prob_lo); + void ComputeDiffusion(MultiFab& S_rhs, MultiFab& S_data, const Real& Dx, From fa807755b277ab5c4e28f8482ed448c7bb76c1cc Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 17 Nov 2025 09:28:45 -0800 Subject: [PATCH 08/15] turn ERK back on --- .../SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials index ca8dc102..dadf9588 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials @@ -22,9 +22,9 @@ diffCoeffy = 1. ## 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 = ForwardEuler #integration.type = RungeKutta -#integration.type = SUNDIALS +integration.type = SUNDIALS # Set the SUNDIALS method type: # ERK = Explicit Runge-Kutta method @@ -37,5 +37,3 @@ integration.type = ForwardEuler # Optionally select a specific SUNDIALS method by name, see the SUNDIALS # documentation for the supported method names integration.sundials.type = ERK - - From d897279b464a66fc7bf593fe5ed066e0688179cc Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Thu, 18 Dec 2025 16:10:53 -0800 Subject: [PATCH 09/15] error WIP, needs periodic wrap --- .../AdvDiff-HyprePrecon/Source/main.cpp | 32 ++++++++++++++++--- .../AdvDiff-HyprePrecon/Source/myfunc.H | 6 +++- 2 files changed, 32 insertions(+), 6 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index 9f348365..8eb79f36 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -149,7 +149,7 @@ void main_main () // ********************************** // INITIALIZE DATA - InitializeData(phi,dx,prob_lo); + 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) @@ -230,11 +230,32 @@ void main_main () 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 + MultiFab phi_exact(ba, dm, Ncomp, 0); + InitializeData(phi_exact,dx,prob_lo,prob_hi,time,advCoeffx,advCoeffy); + const std::string& pltfile = amrex::Concatenate("exact",nsteps,5); + WriteSingleLevelPlotfile(pltfile, phi_exact, {"phi"}, geom, time, nsteps); + + + MultiFab::Subtract(phi_exact,phi,0,0,1,0); + Real error = phi_exact.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_lo, + const GpuArray prob_hi, + const Real& time, + const Real& Ax, + const Real& Ay) { + + GpuArray L; + for (int d=0; d dx, - const GpuArray prob_lo); + 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, From afe2aab9d09af7ced24ee2ff9815d302ea3e3889 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Thu, 18 Dec 2025 16:34:37 -0800 Subject: [PATCH 10/15] progress toward exact solution. Error rates in time not coming out correct. Diffusion not implemented in exact solution --- .../AdvDiff-HyprePrecon/Exec/inputs_sundials | 9 ++-- .../AdvDiff-HyprePrecon/Source/main.cpp | 41 ++++++++++++------- 2 files changed, 30 insertions(+), 20 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials index dadf9588..1f1a4cd1 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials @@ -1,10 +1,9 @@ -n_cell = 32 -max_grid_size = 16 +n_cell = 256 +max_grid_size = 128 +dt = 1.e-5 nsteps = 100 -plot_int = 10 - -dt = 1.e-4 +plot_int = 100 advCoeffx = 100. advCoeffy = 100. diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index 8eb79f36..a5f2279a 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -96,13 +96,6 @@ void main_main () // ********************************** // 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)); @@ -110,8 +103,9 @@ void main_main () // 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" - ba.define(domain); + BoxArray ba(domain); // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction ba.maxSize(max_grid_size); @@ -123,8 +117,10 @@ void main_main () // 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 - geom.define(domain, real_box, CoordSys::cartesian, is_periodic); + Geometry geom(domain, real_box, CoordSys::cartesian, is_periodic); // extract dx from the geometry object GpuArray dx = geom.CellSizeArray(); @@ -232,15 +228,28 @@ void main_main () amrex::Print() << "Total evolution time = " << evolution_stop_time << " seconds\n"; // exact solution - MultiFab phi_exact(ba, dm, Ncomp, 0); + 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); - const std::string& pltfile = amrex::Concatenate("exact",nsteps,5); - WriteSingleLevelPlotfile(pltfile, phi_exact, {"phi"}, geom, time, nsteps); + 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,phi,0,0,1,0); - Real error = phi_exact.norm1(0,geom.periodicity()); + 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; } @@ -251,6 +260,8 @@ void InitializeData(MultiFab& phi, 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; From 315721ea502cdfd31120bd27089e31cf047b8f0f Mon Sep 17 00:00:00 2001 From: Cody Balos Date: Wed, 24 Jun 2026 17:05:02 -0700 Subject: [PATCH 11/15] use new preconditioner routines; add MLMG copy of example --- Docs/source/SUNDIALS_Tutorial.rst | 9 + .../AdvDiff-HyprePrecon/CMakeLists.txt | 2 +- .../AdvDiff-HyprePrecon/Exec/GNUmakefile | 1 + .../AdvDiff-HyprePrecon/Exec/README_sundials | 3 +- .../AdvDiff-HyprePrecon/Exec/inputs_sundials | 10 +- .../AdvDiff-HyprePrecon/Source/main.cpp | 103 +++- .../AdvDiff-MLMGPrecon/CMakeLists.txt | 15 + .../AdvDiff-MLMGPrecon/Exec/GNUmakefile | 23 + .../AdvDiff-MLMGPrecon/Exec/README_sundials | 6 + .../AdvDiff-MLMGPrecon/Exec/inputs_sundials | 57 ++ .../AdvDiff-MLMGPrecon/Source/Make.package | 2 + .../AdvDiff-MLMGPrecon/Source/main.cpp | 491 ++++++++++++++++++ .../AdvDiff-MLMGPrecon/Source/myfunc.H | 30 ++ 13 files changed, 747 insertions(+), 5 deletions(-) create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/CMakeLists.txt create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/GNUmakefile create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/inputs_sundials create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/Make.package create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp create mode 100644 ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/myfunc.H diff --git a/Docs/source/SUNDIALS_Tutorial.rst b/Docs/source/SUNDIALS_Tutorial.rst index 1f8ca4e3..6ac4a72c 100644 --- a/Docs/source/SUNDIALS_Tutorial.rst +++ b/Docs/source/SUNDIALS_Tutorial.rst @@ -22,6 +22,15 @@ 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. + +The advection-diffusion examples at +``amrex-tutorials/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon`` and +``amrex-tutorials/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon`` demonstrate +SUNDIALS with preconditioning for an implicit diffusion solve using HYPRE and +AMReX MLMG, respectively. The MLMG example can also use MLMG's HYPRE bottom +solver interface when AMReX is built with HYPRE support. Their inputs also +show how to select the SUNDIALS SPGMR +preconditioning side with ``integration.sundials.linear_solver_preconditioning``. Please see each respective inputs file or `AMReX User Guide:Time Integration`_ for more details. diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt index 71eb16ae..7814cb07 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt @@ -1,4 +1,4 @@ -if (AMReX_SPACEDIM EQUAL 1) +if (AMReX_SPACEDIM EQUAL 1 OR NOT AMReX_HYPRE) return() endif () diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile index 657fa127..af311ed7 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile @@ -7,6 +7,7 @@ USE_OMP = FALSE COMP = gnu DIM = 2 USE_SUNDIALS = TRUE +USE_HYPRE = TRUE SUNDIALS_HOME ?= ../../../../../sundials/instdir include $(AMREX_HOME)/Tools/GNUMake/Make.defs diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials index 3b6830a6..abbd9116 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials @@ -1,4 +1,5 @@ Refer to https://amrex-codes.github.io/amrex/docs_html/TimeIntegration_Chapter.html -for SNUDIALS installation, build, and usage instructions. +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 HYPRE support. diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials index 1f1a4cd1..3b135e2d 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials @@ -35,4 +35,12 @@ integration.type = SUNDIALS # # Optionally select a specific SUNDIALS method by name, see the SUNDIALS # documentation for the supported method names -integration.sundials.type = ERK +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 + +# HYPRE settings for the SUNDIALS preconditioner solve +hypre.hypre_solver = BoomerAMG diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index a5f2279a..3fa5d63d 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -1,9 +1,15 @@ #include #include #include +#include #include "myfunc.H" +#include +#include +#include +#include + using namespace amrex; int main (int argc, char* argv[]) @@ -88,9 +94,9 @@ void main_main () pp.query("abstol",abstol); pp.query("advCoeffx",advCoeffx); - pp.query("advCoeffx",advCoeffy); + pp.query("advCoeffy",advCoeffy); pp.query("diffCoeffx",diffCoeffx); - pp.query("diffCoeffx",diffCoeffy); + pp.query("diffCoeffy",diffCoeffy); } // ********************************** @@ -186,9 +192,102 @@ void main_main () ComputeAdvection(S_rhs, S_data, advCoeffx, advCoeffy, dx); }; + constexpr int hypre_stencil_size = 2 * AMREX_SPACEDIM + 1; + using HyprePreconditioner = HypreSolver; + + std::unique_ptr hypre_preconditioner; + Real hypre_gamma = std::numeric_limits::quiet_NaN(); + + auto build_hypre_preconditioner = [&](Real gamma) + { + const IntVect domain_lo = geom.Domain().smallEnd(); + const IntVect domain_hi = geom.Domain().bigEnd(); + const HYPRE_Real hx = static_cast(gamma * diffCoeffx / (dx[0] * dx[0])); + const HYPRE_Real hy = static_cast(gamma * diffCoeffy / (dx[1] * dx[1])); + const HYPRE_Real diag = static_cast(1.0) + 2.0 * (hx + hy); + + auto marker = [=] AMREX_GPU_DEVICE (int /* boxno */, int /* i */, int /* j */, + int /* k */, int /* n */) -> bool + { + return true; + }; + + auto filler = [=] AMREX_GPU_DEVICE (int /* boxno */, int i, int j, int k, int n, + Array4 const* gid, + HYPRE_Int& ncols, HYPRE_Int* cols, + HYPRE_Real* mat) + { + const int ilo = domain_lo[0]; + const int ihi = domain_hi[0]; + const int jlo = domain_lo[1]; + const int jhi = domain_hi[1]; + + const int im = (i == ilo) ? ihi : i - 1; + const int ip = (i == ihi) ? ilo : i + 1; + const int jm = (j == jlo) ? jhi : j - 1; + const int jp = (j == jhi) ? jlo : j + 1; + + ncols = 0; + + cols[ncols] = gid[n](im,j,k); + mat [ncols] = -hx; + ++ncols; + + cols[ncols] = gid[n](ip,j,k); + mat [ncols] = -hx; + ++ncols; + + cols[ncols] = gid[n](i,jm,k); + mat [ncols] = -hy; + ++ncols; + + cols[ncols] = gid[n](i,jp,k); + mat [ncols] = -hy; + ++ncols; + + cols[ncols] = gid[n](i,j,k); + mat [ncols] = diag; + ++ncols; + }; + + hypre_preconditioner = std::make_unique( + Vector{IndexType::TheCellType()}, + IntVect(1), geom, ba, dm, marker, filler); + hypre_gamma = gamma; + }; + + auto precond_setup = [&](MultiFab& /* S_data */, MultiFab& /* S_rhs */, const Real /* time */, + bool jok, bool& jcur, const Real gamma) + { + const bool same_gamma = hypre_preconditioner && + std::abs(gamma - hypre_gamma) <= + (10.0 * std::numeric_limits::epsilon() * + std::max(Real(1.0), std::max(std::abs(gamma), std::abs(hypre_gamma)))); + + if (!(jok && same_gamma)) { + build_hypre_preconditioner(gamma); + jcur = true; + } else { + jcur = false; + } + }; + + 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(hypre_preconditioner != nullptr); + S_soln.setVal(0.0); + hypre_preconditioner->solve(Vector{&S_soln}, + Vector{&S_rhs}, + 0.0_rt, 0.0_rt, 1); + }; + 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(); 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..02cc12df --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials @@ -0,0 +1,6 @@ +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. +Set USE_HYPRE=TRUE and mlmg.use_hypre=1 to use MLMG's HYPRE bottom solver interface. 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..ef2a377e --- /dev/null +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp @@ -0,0 +1,491 @@ +#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::unique_ptr mlmg_preconditioner; + Real mlmg_gamma = std::numeric_limits::quiet_NaN(); + + 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); + + std::array face_bcoef; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + const BoxArray& face_ba = amrex::convert(ba, IntVect::TheDimensionVector(idim)); + face_bcoef[idim].define(face_ba, dm, 1, 0); + const Real diff_coeff = (idim == 0) ? diffCoeffx : ((idim == 1) ? diffCoeffy : 0.0); + face_bcoef[idim].setVal(gamma * diff_coeff); + } + preconditioner->linop->setBCoeffs(0, amrex::GetArrOfConstPtrs(face_bcoef)); + + 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) + { + const bool same_gamma = mlmg_preconditioner && + std::abs(gamma - mlmg_gamma) <= + (10.0 * std::numeric_limits::epsilon() * + std::max(Real(1.0), std::max(std::abs(gamma), std::abs(mlmg_gamma)))); + + if (!(jok && same_gamma)) { + build_mlmg_preconditioner(gamma); + jcur = true; + } else { + jcur = false; + } + }; + + 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 From a61bd2fc5f118f7e7098c0b4301325f1033f2290 Mon Sep 17 00:00:00 2001 From: Cody Balos Date: Thu, 25 Jun 2026 10:59:35 -0700 Subject: [PATCH 12/15] only rebuild if no jok --- ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp | 4 +--- ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp | 4 +--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp index 3fa5d63d..e9fa225b 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp @@ -264,11 +264,9 @@ void main_main () (10.0 * std::numeric_limits::epsilon() * std::max(Real(1.0), std::max(std::abs(gamma), std::abs(hypre_gamma)))); - if (!(jok && same_gamma)) { + if (!jok) { build_hypre_preconditioner(gamma); jcur = true; - } else { - jcur = false; } }; diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp index ef2a377e..3100ed92 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp @@ -293,11 +293,9 @@ void main_main () (10.0 * std::numeric_limits::epsilon() * std::max(Real(1.0), std::max(std::abs(gamma), std::abs(mlmg_gamma)))); - if (!(jok && same_gamma)) { + if (!jok) { build_mlmg_preconditioner(gamma); jcur = true; - } else { - jcur = false; } }; From af255468b6f4af5794ff9e4f80223518792f0115 Mon Sep 17 00:00:00 2001 From: Cody Balos Date: Tue, 30 Jun 2026 16:16:29 -0700 Subject: [PATCH 13/15] remove hypre direct interface example --- .../AdvDiff-HyprePrecon/CMakeLists.txt | 15 - .../AdvDiff-HyprePrecon/Exec/GNUmakefile | 21 - .../AdvDiff-HyprePrecon/Exec/README_sundials | 5 - .../AdvDiff-HyprePrecon/Exec/inputs_sundials | 46 -- .../AdvDiff-HyprePrecon/Source/Make.package | 2 - .../AdvDiff-HyprePrecon/Source/main.cpp | 462 ------------------ .../AdvDiff-HyprePrecon/Source/myfunc.H | 30 -- 7 files changed, 581 deletions(-) delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp delete mode 100644 ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt deleted file mode 100644 index 7814cb07..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/CMakeLists.txt +++ /dev/null @@ -1,15 +0,0 @@ -if (AMReX_SPACEDIM EQUAL 1 OR NOT AMReX_HYPRE) - 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-HyprePrecon/Exec/GNUmakefile b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile deleted file mode 100644 index af311ed7..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/GNUmakefile +++ /dev/null @@ -1,21 +0,0 @@ -# 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 = TRUE -SUNDIALS_HOME ?= ../../../../../sundials/instdir - -include $(AMREX_HOME)/Tools/GNUMake/Make.defs - -include ../Source/Make.package -VPATH_LOCATIONS += ../Source -INCLUDE_LOCATIONS += ../Source - -include $(AMREX_HOME)/Src/Base/Make.package - -include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials deleted file mode 100644 index abbd9116..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/README_sundials +++ /dev/null @@ -1,5 +0,0 @@ -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 HYPRE support. diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials deleted file mode 100644 index 3b135e2d..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Exec/inputs_sundials +++ /dev/null @@ -1,46 +0,0 @@ -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 - -# HYPRE settings for the SUNDIALS preconditioner solve -hypre.hypre_solver = BoomerAMG diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package deleted file mode 100644 index fa194a35..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/Make.package +++ /dev/null @@ -1,2 +0,0 @@ -CEXE_sources += main.cpp -CEXE_headers += myfunc.H diff --git a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp deleted file mode 100644 index e9fa225b..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/main.cpp +++ /dev/null @@ -1,462 +0,0 @@ -#include -#include -#include -#include - -#include "myfunc.H" - -#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; - - // 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); - } - - // ********************************** - // 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); - }; - - constexpr int hypre_stencil_size = 2 * AMREX_SPACEDIM + 1; - using HyprePreconditioner = HypreSolver; - - std::unique_ptr hypre_preconditioner; - Real hypre_gamma = std::numeric_limits::quiet_NaN(); - - auto build_hypre_preconditioner = [&](Real gamma) - { - const IntVect domain_lo = geom.Domain().smallEnd(); - const IntVect domain_hi = geom.Domain().bigEnd(); - const HYPRE_Real hx = static_cast(gamma * diffCoeffx / (dx[0] * dx[0])); - const HYPRE_Real hy = static_cast(gamma * diffCoeffy / (dx[1] * dx[1])); - const HYPRE_Real diag = static_cast(1.0) + 2.0 * (hx + hy); - - auto marker = [=] AMREX_GPU_DEVICE (int /* boxno */, int /* i */, int /* j */, - int /* k */, int /* n */) -> bool - { - return true; - }; - - auto filler = [=] AMREX_GPU_DEVICE (int /* boxno */, int i, int j, int k, int n, - Array4 const* gid, - HYPRE_Int& ncols, HYPRE_Int* cols, - HYPRE_Real* mat) - { - const int ilo = domain_lo[0]; - const int ihi = domain_hi[0]; - const int jlo = domain_lo[1]; - const int jhi = domain_hi[1]; - - const int im = (i == ilo) ? ihi : i - 1; - const int ip = (i == ihi) ? ilo : i + 1; - const int jm = (j == jlo) ? jhi : j - 1; - const int jp = (j == jhi) ? jlo : j + 1; - - ncols = 0; - - cols[ncols] = gid[n](im,j,k); - mat [ncols] = -hx; - ++ncols; - - cols[ncols] = gid[n](ip,j,k); - mat [ncols] = -hx; - ++ncols; - - cols[ncols] = gid[n](i,jm,k); - mat [ncols] = -hy; - ++ncols; - - cols[ncols] = gid[n](i,jp,k); - mat [ncols] = -hy; - ++ncols; - - cols[ncols] = gid[n](i,j,k); - mat [ncols] = diag; - ++ncols; - }; - - hypre_preconditioner = std::make_unique( - Vector{IndexType::TheCellType()}, - IntVect(1), geom, ba, dm, marker, filler); - hypre_gamma = gamma; - }; - - auto precond_setup = [&](MultiFab& /* S_data */, MultiFab& /* S_rhs */, const Real /* time */, - bool jok, bool& jcur, const Real gamma) - { - const bool same_gamma = hypre_preconditioner && - std::abs(gamma - hypre_gamma) <= - (10.0 * std::numeric_limits::epsilon() * - std::max(Real(1.0), std::max(std::abs(gamma), std::abs(hypre_gamma)))); - - if (!jok) { - build_hypre_preconditioner(gamma); - jcur = true; - } - }; - - 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(hypre_preconditioner != nullptr); - S_soln.setVal(0.0); - hypre_preconditioner->solve(Vector{&S_soln}, - Vector{&S_rhs}, - 0.0_rt, 0.0_rt, 1); - }; - - 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-HyprePrecon/Source/myfunc.H b/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H deleted file mode 100644 index 0879e55d..00000000 --- a/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon/Source/myfunc.H +++ /dev/null @@ -1,30 +0,0 @@ -#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 From 3bb3160c684223436db7dd58e10048ea1f5b1cbf Mon Sep 17 00:00:00 2001 From: Cody Balos Date: Fri, 26 Jun 2026 15:53:10 -0700 Subject: [PATCH 14/15] finer control of precond update --- .../AdvDiff-MLMGPrecon/Source/main.cpp | 45 ++++++++++++++----- 1 file changed, 35 insertions(+), 10 deletions(-) diff --git a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp index 3100ed92..2d9693b9 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Source/main.cpp @@ -229,11 +229,32 @@ void main_main () 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(); @@ -249,14 +270,11 @@ void main_main () preconditioner->linop->setScalars(1.0, 1.0); preconditioner->linop->setACoeffs(0, 1.0); - std::array face_bcoef; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { const BoxArray& face_ba = amrex::convert(ba, IntVect::TheDimensionVector(idim)); - face_bcoef[idim].define(face_ba, dm, 1, 0); - const Real diff_coeff = (idim == 0) ? diffCoeffx : ((idim == 1) ? diffCoeffy : 0.0); - face_bcoef[idim].setVal(gamma * diff_coeff); + preconditioner->face_bcoef[idim].define(face_ba, dm, 1, 0); } - preconditioner->linop->setBCoeffs(0, amrex::GetArrOfConstPtrs(face_bcoef)); + update_mlmg_gamma(*preconditioner, gamma); preconditioner->solver = std::make_unique(*preconditioner->linop); preconditioner->solver->setMaxIter(mlmg_max_iter); @@ -288,14 +306,21 @@ void main_main () auto precond_setup = [&](MultiFab& /* S_data */, MultiFab& /* S_rhs */, const Real /* time */, bool jok, bool& jcur, const Real gamma) { - const bool same_gamma = mlmg_preconditioner && - std::abs(gamma - mlmg_gamma) <= - (10.0 * std::numeric_limits::epsilon() * - std::max(Real(1.0), std::max(std::abs(gamma), std::abs(mlmg_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) { - build_mlmg_preconditioner(gamma); + 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); } }; From 70bc2feb110b5c0c36ee35e603d0aabe0b569442 Mon Sep 17 00:00:00 2001 From: Cody Balos Date: Tue, 30 Jun 2026 16:49:54 -0700 Subject: [PATCH 15/15] udpate docs --- Docs/source/SUNDIALS_Tutorial.rst | 21 +++++++++++-------- .../AdvDiff-MLMGPrecon/Exec/README_sundials | 7 ++++++- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/Docs/source/SUNDIALS_Tutorial.rst b/Docs/source/SUNDIALS_Tutorial.rst index 6ac4a72c..2af438ec 100644 --- a/Docs/source/SUNDIALS_Tutorial.rst +++ b/Docs/source/SUNDIALS_Tutorial.rst @@ -23,15 +23,18 @@ 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. -The advection-diffusion examples at -``amrex-tutorials/ExampleCodes/SUNDIALS/AdvDiff-HyprePrecon`` and -``amrex-tutorials/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon`` demonstrate -SUNDIALS with preconditioning for an implicit diffusion solve using HYPRE and -AMReX MLMG, respectively. The MLMG example can also use MLMG's HYPRE bottom -solver interface when AMReX is built with HYPRE support. Their inputs also -show how to select the SUNDIALS SPGMR -preconditioning side with ``integration.sundials.linear_solver_preconditioning``. -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/Exec/README_sundials b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials index 02cc12df..767ba43f 100644 --- a/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials +++ b/ExampleCodes/SUNDIALS/AdvDiff-MLMGPrecon/Exec/README_sundials @@ -3,4 +3,9 @@ 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. -Set USE_HYPRE=TRUE and mlmg.use_hypre=1 to use MLMG's HYPRE bottom solver interface. +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.