diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler b/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler index eeec44c6..39b09830 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler +++ b/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler @@ -12,23 +12,3 @@ dt = 1.e-5 ## 1 or "RungeKutta" => Native AMReX Explicit Runge Kutta controlled by integration.rk.type ## 2 or "SUNDIALS" => SUNDIALS backend controlled by integration.sundials.strategy integration.type = ForwardEuler - -## Native AMReX Explicit Runge-Kutta parameters -# -## integration.rk.type can take the following values: -### 0 = User-specified Butcher Tableau -### 1 = Forward Euler -### 2 = Trapezoid Method -### 3 = SSPRK3 Method -### 4 = RK4 Method -integration.rk.type = 3 - -## If using the SUNDIALS Submodule, then -## compile with USE_SUNDIALS=TRUE or AMReX_SUNDIALS=ON and -## set strategy here: -# -## integration.sundials.strategy can take the following values: -### ERK = ERKStep from ARKode with SSPRK3 Method -### MRI = MRIStep from ARKode with Explict Trapezoid Method -### MRITEST = MRIStep from ARKode modified to use no-op inner f0 -integration.sundials.strategy = ERK diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_rk3 b/ExampleCodes/SUNDIALS/Exec/inputs_rk3 index dded52b9..6eb60137 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_rk3 +++ b/ExampleCodes/SUNDIALS/Exec/inputs_rk3 @@ -22,13 +22,3 @@ integration.type = RungeKutta ### 3 = SSPRK3 Method ### 4 = RK4 Method integration.rk.type = 3 - -## If using the SUNDIALS Submodule, then -## compile with USE_SUNDIALS=TRUE or AMReX_SUNDIALS=ON and -## set strategy here: -# -## integration.sundials.strategy can take the following values: -### ERK = ERKStep from ARKode with SSPRK3 Method -### MRI = MRIStep from ARKode with Explict Trapezoid Method -### MRITEST = MRIStep from ARKode modified to use no-op inner f0 -integration.sundials.strategy = ERK diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials b/ExampleCodes/SUNDIALS/Exec/inputs_sundials new file mode 100644 index 00000000..bbdf2c69 --- /dev/null +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials @@ -0,0 +1,44 @@ +n_cell = 32 +max_grid_size = 16 + +nsteps = 1000 +plot_int = 100 + +dt = 1.e-5 + +# Use adaptive time stepping and set integrator relative and absolute tolerances +# adapt_dt = true +# reltol = 1.0e-4 +# abstol = 1.0e-9 + +# INTEGRATION +# integration.type can take on the following values: +# 0 or "ForwardEuler" => Native AMReX Forward Euler integrator +# 1 or "RungeKutta" => Native AMReX Explicit Runge Kutta controlled by integration.rk.type +# 2 or "SUNDIALS" => SUNDIALS backend controlled by integration.sundials.type +# +# If using the SUNDIALS Submodule, then compile with USE_SUNDIALS=TRUE or +# AMReX_SUNDIALS=ON +integration.type = SUNDIALS + +# Set the SUNDIALS method type: +# ERK = Explicit Runge-Kutta method +# DIRK = Diagonally Implicit Runge-Kutta method +# IMEX-RK = Implicit-Explicit Additive Runge-Kutta method +# EX-MRI = Explicit Multirate Infatesimal method +# IM-MRI = Implicit Multirate Infatesimal method +# IMEX-MRI = Implicit-Explicit Multirate Infatesimal method +# +# Optionally select a specific SUNDIALS method by name, see the SUNDIALS +# documentation for the supported method names + +# Use the SUNDIALS default method for the chosen type (fixed or adaptive step sizes) +integration.sundials.type = ERK + +# Use forward Euler (fixed step sizes only) +# integration.sundials.type = ERK +# integration.sundials.method = ARKODE_FORWARD_EULER_1_1 + +# Use backward Euler (fixed step sizes only) +# integration.sundials.type = DIRK +# integration.sundials.method = ARKODE_BACKWARD_EULER_1_1 diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk deleted file mode 100644 index 59e622aa..00000000 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ /dev/null @@ -1,34 +0,0 @@ -n_cell = 32 -max_grid_size = 16 - -nsteps = 1000 -plot_int = 100 - -dt = 1.e-5 - -# INTEGRATION -## integration.type can take on the following values: -## 0 or "ForwardEuler" => Native AMReX Forward Euler integrator -## 1 or "RungeKutta" => Native AMReX Explicit Runge Kutta controlled by integration.rk.type -## 2 or "SUNDIALS" => SUNDIALS backend controlled by integration.sundials.strategy -integration.type = 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 = 3 - -## If using the SUNDIALS Submodule, then -## compile with USE_SUNDIALS=TRUE or AMReX_SUNDIALS=ON and -## set strategy here: -# -## integration.sundials.strategy can take the following values: -### ERK = ERKStep from ARKode with SSPRK3 Method -### MRI = MRIStep from ARKode with Explict Trapezoid Method -### MRITEST = MRIStep from ARKode modified to use no-op inner f0 -integration.sundials.strategy = ERK \ No newline at end of file diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index 71e7ab83..cd24a35f 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -37,6 +37,13 @@ void main_main () // time step Real dt; + // use adaptive time step (dt used to set output times) + bool adapt_dt = false; + + // adaptive time step relative and absolute tolerances + Real reltol = 1.0e-4; + Real abstol = 1.0e-9; + // inputs parameters { // ParmParse is way of reading inputs from the inputs file @@ -62,6 +69,13 @@ void main_main () // 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); } // ********************************** @@ -110,11 +124,7 @@ void main_main () DistributionMapping dm(ba); // we allocate two phi multifabs; one will store the old state, the other the new. - Vector state_old, state_new; - state_old.push_back(MultiFab(ba, dm, Ncomp, Nghost)); - state_new.push_back(MultiFab(ba, dm, Ncomp, Nghost)); - auto& phi_old = state_old[0]; - auto& phi_new = state_new[0]; + MultiFab phi(ba, dm, Ncomp, Nghost); // time = starting time in the simulation Real time = 0.0; @@ -123,11 +133,11 @@ void main_main () // INITIALIZE DATA // loop over boxes - for (MFIter mfi(phi_old); mfi.isValid(); ++mfi) + for (MFIter mfi(phi); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); - const Array4& phiOld = phi_old.array(mfi); + 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) @@ -140,7 +150,7 @@ void main_main () 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 - phiOld(i,j,k) = 1. + std::exp(-rsquared); + phi_array(i,j,k) = 1. + std::exp(-rsquared); }); } @@ -149,56 +159,58 @@ void main_main () { int step = 0; const std::string& pltfile = amrex::Concatenate("plt",step,5); - WriteSingleLevelPlotfile(pltfile, phi_old, {"phi"}, geom, time, 0); + WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, 0); } - // fill periodic ghost cells - phi_old.FillBoundary(geom.periodicity()); + auto pre_rhs_function = [&](MultiFab& S_data, const Real /* time */) { + // fill periodic ghost cells + S_data.FillBoundary(geom.periodicity()); + }; + + auto rhs_function = [&](MultiFab& S_rhs, const MultiFab& S_data, + const Real /* time */) { + + // loop over boxes + auto& phi_data = S_data; + auto& phi_rhs = S_rhs; + 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) = ( (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 + ); + }); + } + }; + + TimeIntegrator integrator(phi, time); + integrator.set_pre_rhs_action(pre_rhs_function); + integrator.set_rhs(rhs_function); + if (adapt_dt) { + integrator.set_adaptive_step(); + integrator.set_tolerances(reltol, abstol); + } else { + integrator.set_time_step(dt); + } for (int step = 1; step <= nsteps; ++step) { - auto rhs_function = [&](Vector& S_rhs, - const Vector& S_data, const Real /* time */) { - // loop over boxes - auto& phi_data = S_data[0]; - auto& phi_rhs = S_rhs[0]; - 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) = ( (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 post_update_function = [&](Vector& S_data, const Real /* time */) { - // fill periodic ghost cells - S_data[0].FillBoundary(geom.periodicity()); - }; - - TimeIntegrator > integrator(state_old); - - integrator.set_rhs(rhs_function); - integrator.set_post_update(post_update_function); - - // advance from state_old at time to state_new at time + dt - integrator.advance(state_old, state_new, time, dt); - - // swap old/new and update time by dt - std::swap(state_old, state_new); + // Set time to evolve to time += dt; + // Advance to output time + integrator.evolve(phi, time); + // Tell the I/O Processor to write out which step we're doing amrex::Print() << "Advanced step " << step << "\n"; @@ -206,7 +218,7 @@ void main_main () if (plot_int > 0 && step%plot_int == 0) { const std::string& pltfile = amrex::Concatenate("plt",step,5); - WriteSingleLevelPlotfile(pltfile, phi_new, {"phi"}, geom, time, step); + WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, step); } } }