Merge pull request #127 from ajnonaka/sundials_work

SUNDIALS Work

ExampleCodes/SUNDIALS/Reaction-Diffusion/Exec/inputs_sundials_mri

@@ -12,13 +12,12 @@ max_grid_size = 16
 nsteps = 12 # 60 # 600
 plot_int = 2 # 10 # 100
 
-dt = 8.5e-3 # 1.7e-3 #1.7e-4
-
 # 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
@@ -29,19 +28,18 @@ use_MRI = true
 fast_dt_ratio = 0.02
 
 
-# Use adaptive time stepping and set integrator relative and absolute tolerances
+# 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
-#
-# If using the SUNDIALS Submodule, then compile with USE_SUNDIALS=TRUE or
-# AMReX_SUNDIALS=ON
+## 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:

ExampleCodes/SUNDIALS/Reaction-Diffusion/Source/main.cpp

@@ -48,6 +48,10 @@ void main_main ()
  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
@@ -85,6 +89,8 @@ void main_main ()
  pp.query("reltol",reltol);
  pp.query("abstol",abstol);
 
+ pp.query("reaction_coef",reaction_coef);
+ pp.query("diffusion_coef",diffusion_coef);
  }
 
  // **********************************
@@ -132,7 +138,7 @@ void main_main ()
  // How Boxes are distrubuted among MPI processes
  DistributionMapping dm(ba);
 
- // we allocate two phi multifabs; one will store the old state, the other the new.
+ // allocate phi MultiFab
  MultiFab phi(ba, dm, Ncomp, Nghost);
 
  // time = starting time in the simulation
@@ -171,8 +177,7 @@ void main_main ()
  WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, 0);
  }
 
- auto rhs_fast_function = [&](MultiFab& S_rhs, MultiFab& S_data,
- const Real /* time */) {
+ auto rhs_fast_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) {
 
  // fill periodic ghost cells
  S_data.FillBoundary(geom.periodicity());
@@ -181,12 +186,6 @@ void main_main ()
  auto& phi_data = S_data;
  auto& phi_rhs = S_rhs;
 
- //Reaction and Diffusion Coefficients
- Real reaction_coef = 0.0;
- Real diffusion_coef = 1.0;
- ParmParse pp;
- pp.query("diffusion_coef",diffusion_coef);
-
  for ( MFIter mfi(phi_data); mfi.isValid(); ++mfi )
  {
  const Box& bx = mfi.validbox();
@@ -200,14 +199,14 @@ void main_main ()
  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]) );
+ +(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 */) {
+ auto rhs_function = [&](MultiFab& S_rhs, MultiFab& S_data, const Real /* time */) {
 
  // fill periodic ghost cells
  S_data.FillBoundary(geom.periodicity());
@@ -216,13 +215,6 @@ void main_main ()
  auto& phi_data = S_data;
  auto& phi_rhs = S_rhs;
 
- //Reaction and Diffusion Coefficients
- Real reaction_coef = 0.0;
- Real diffusion_coef = 1.0;
- ParmParse pp;
- pp.query("reaction_coef",reaction_coef);
- pp.query("diffusion_coef",diffusion_coef);
-
  for ( MFIter mfi(phi_data); mfi.isValid(); ++mfi )
  {
  const Box& bx = mfi.validbox();
@@ -239,15 +231,14 @@ void main_main ()
  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]) )
+ +(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);
+ ) - 1.*reaction_coef*phi_array(i,j,k);
  }
  });
  }
  };
 
-
  TimeIntegrator<MultiFab> integrator(phi, time);
  integrator.set_rhs(rhs_function);
  if (use_MRI) {

ExampleCodes/SUNDIALS/Exec/GNUmakefile → ExampleCodes/SUNDIALS/Single-Rate/Exec/GNUmakefile

@@ -1,5 +1,5 @@
 # AMREX_HOME defines the directory in which we will find all the AMReX code.
-AMREX_HOME ?= ../../../../amrex
+AMREX_HOME ?= ../../../../../amrex
 
 DEBUG = FALSE
 USE_MPI = TRUE
@@ -17,9 +17,9 @@ INCLUDE_LOCATIONS += ../Source
 
 ifeq ($(USE_SUNDIALS),TRUE)
  ifeq ($(USE_CUDA),TRUE)
- SUNDIALS_ROOT ?= $(TOP)../../../../sundials/instdir_cuda
+ SUNDIALS_ROOT ?= $(TOP)../../../../../sundials/instdir_cuda
  else
- SUNDIALS_ROOT ?= $(TOP)../../../../sundials/instdir
+ SUNDIALS_ROOT ?= $(TOP)../../../../../sundials/instdir
  endif
  ifeq ($(NERSC_HOST),perlmutter)
  SUNDIALS_LIB_DIR ?= $(SUNDIALS_ROOT)/lib64

ExampleCodes/SUNDIALS/Single-Rate/Exec/inputs

@@ -0,0 +1,53 @@
+n_cell = 32
+max_grid_size = 16
+
+nsteps = 5
+plot_int = 5
+dt = 1.e-4
+
+#nsteps = 10
+#plot_int = 10
+#dt = 5.e-5
+
+#nsteps = 20
+#plot_int = 20
+#dt = 2.5e-5
+
+# 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
+
+## Native AMReX Explicit Runge-Kutta parameters
+#
+## integration.rk.type can take the following values:
+### 0 = User-specified Butcher Tableau
+### 1 = Forward Euler
+### 2 = Trapezoid Method
+### 3 = SSPRK3 Method
+### 4 = RK4 Method
+integration.rk.type = 1
+
+# Set the SUNDIALS method type:
+# ERK = Explicit Runge-Kutta method
+# DIRK = Diagonally Implicit Runge-Kutta method
+#
+# Optionally select a specific SUNDIALS method by name, see the SUNDIALS
+# documentation for the supported method names
+
+# Use forward Euler (fixed step sizes only)
+integration.sundials.type = ERK
+integration.sundials.method = ARKODE_FORWARD_EULER_1_1
+
+# Use backward Euler (fixed step sizes only)
+#integration.sundials.type = DIRK
+#integration.sundials.method = ARKODE_BACKWARD_EULER_1_1

ExampleCodes/SUNDIALS/Source/main.cpp → ExampleCodes/SUNDIALS/Single-Rate/Source/main.cpp

@@ -123,7 +123,7 @@ void main_main ()
  // How Boxes are distrubuted among MPI processes
  DistributionMapping dm(ba);
 
- // we allocate two phi multifabs; one will store the old state, the other the new.
+ // allocate phi MultiFab
  MultiFab phi(ba, dm, Ncomp, Nghost);
 
  // time = starting time in the simulation
@@ -162,15 +162,15 @@ void main_main ()
  WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, 0);
  }
 
- auto rhs_function = [&](MultiFab& S_rhs, MultiFab& S_data,
- const Real /* time */) {
+ 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();
@@ -181,7 +181,7 @@ void main_main ()
  // 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) = ( (phi_array(i+1,j,k) - 2.*phi_array(i,j,k) + phi_array(i-1,j,k)) / (dx[0]*dx[0])
+ phi_rhs_array(i,j,k) = (  (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])
@@ -200,16 +200,23 @@ void main_main ()
  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 << "\n";
+ 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)
@@ -218,4 +225,8 @@ void main_main ()
  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";
 }