From 8388c320b7a522c2e084a37b45e50cba7add07af Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Thu, 30 May 2024 00:05:48 -0700 Subject: [PATCH 01/13] Update sundials example --- .../SUNDIALS/Exec/inputs_sundials_erk | 35 ++++------ ExampleCodes/SUNDIALS/Source/main.cpp | 70 +++++++++---------- 2 files changed, 49 insertions(+), 56 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk index 59e622aa..f497f742 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk @@ -7,28 +7,21 @@ 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 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 = SUNDIALS -## Native AMReX Explicit Runge-Kutta parameters +# If using the SUNDIALS Submodule, then compile with USE_SUNDIALS=TRUE or +# AMReX_SUNDIALS=ON and set method type here: # -## 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 +# ERK = Explicit Runge-Kutta method +# DIRK = Diagonally Implicit Runge-Kutta method +# IMEX = Implicit-Explicit Additive Runge-Kutta method +# MRI = Multirate Infatesimal Method +integration.sundials.type = ERK -## 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 +# To select a specific SUNDIALS ERK or DIRK method set integration.sundials.method +# to a method name from https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html +integration.sundials.method = ARKODE_FORWARD_EULER_1_1 diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index 71e7ab83..04ba7310 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -155,43 +155,43 @@ void main_main () // fill periodic ghost cells phi_old.FillBoundary(geom.periodicity()); - 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 ) + auto rhs_function = [&](Vector& S_rhs, + 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) { - 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); + 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, time); + integrator.set_rhs(rhs_function); + integrator.set_post_update(post_update_function); + + for (int step = 1; step <= nsteps; ++step) + { // advance from state_old at time to state_new at time + dt integrator.advance(state_old, state_new, time, dt); From 88bfc422113dd101e307e62156180e715768fe55 Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Fri, 31 May 2024 00:50:30 -0700 Subject: [PATCH 02/13] adaptive stepping updates --- .../SUNDIALS/Exec/inputs_sundials_erk | 13 ++++-- ExampleCodes/SUNDIALS/Source/main.cpp | 41 +++++++++++-------- 2 files changed, 32 insertions(+), 22 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk index f497f742..3a5a02c5 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk @@ -1,10 +1,15 @@ n_cell = 32 max_grid_size = 16 -nsteps = 1000 -plot_int = 100 +# time between outputs +dt_out = 1.e-3 -dt = 1.e-5 +# number of output steps +nsteps = 10 +plot_int = 10 + +# final time +t_final = 1.e-2 # INTEGRATION # integration.type can take on the following values: @@ -24,4 +29,4 @@ integration.sundials.type = ERK # To select a specific SUNDIALS ERK or DIRK method set integration.sundials.method # to a method name from https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html -integration.sundials.method = ARKODE_FORWARD_EULER_1_1 +# integration.sundials.method = ARKODE_FORWARD_EULER_1_1 diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index 04ba7310..352d52d1 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -35,7 +35,10 @@ void main_main () int plot_int; // time step - Real dt; + Real dt_out; + + // final time + Real t_final; // inputs parameters { @@ -61,7 +64,10 @@ void main_main () pp.query("plot_int",plot_int); // time step - pp.get("dt",dt); + pp.get("dt_out",dt_out); + + // time step + pp.get("t_final",t_final); } // ********************************** @@ -110,11 +116,9 @@ 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]; + Vector state; + state.push_back(MultiFab(ba, dm, Ncomp, Nghost)); + auto& phi = state[0]; // time = starting time in the simulation Real time = 0.0; @@ -123,11 +127,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& phiOld = phi.array(mfi); // set phi = 1 + e^(-(r-0.5)^2) amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) @@ -149,11 +153,11 @@ 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()); + phi.FillBoundary(geom.periodicity()); auto rhs_function = [&](Vector& S_rhs, Vector& S_data, const Real /* time */) { @@ -186,18 +190,17 @@ void main_main () S_data[0].FillBoundary(geom.periodicity()); }; - TimeIntegrator > integrator(state_old, time); + TimeIntegrator > integrator(state, time); integrator.set_rhs(rhs_function); integrator.set_post_update(post_update_function); for (int step = 1; step <= nsteps; ++step) { - // advance from state_old at time to state_new at time + dt - integrator.advance(state_old, state_new, time, dt); + // Set output time + time += dt_out; - // swap old/new and update time by dt - std::swap(state_old, state_new); - time += dt; + // Advance to output time + integrator.evolve(state, time); // Tell the I/O Processor to write out which step we're doing amrex::Print() << "Advanced step " << step << "\n"; @@ -206,7 +209,9 @@ 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); } + + if (time > t_final) break; } } From 54fefabb2b515ec7696412be49f13ed1e1a6f67e Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Fri, 31 May 2024 07:53:27 -0700 Subject: [PATCH 03/13] add back const --- ExampleCodes/SUNDIALS/Source/main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index 352d52d1..82d74c55 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -160,7 +160,7 @@ void main_main () phi.FillBoundary(geom.periodicity()); auto rhs_function = [&](Vector& S_rhs, - Vector& S_data, const Real /* time */) { + const Vector& S_data, const Real /* time */) { // loop over boxes auto& phi_data = S_data[0]; From fe51b727a65dae17beb1377f10c9c0550e063b15 Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Sat, 1 Jun 2024 21:22:45 -0700 Subject: [PATCH 04/13] post update to pre rhs update --- ExampleCodes/SUNDIALS/Source/main.cpp | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index 82d74c55..c5518cae 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -156,8 +156,10 @@ void main_main () WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, 0); } - // fill periodic ghost cells - phi.FillBoundary(geom.periodicity()); + auto pre_rhs_function = [&](Vector& S_data, const Real /* time */) { + // fill periodic ghost cells + S_data[0].FillBoundary(geom.periodicity()); + }; auto rhs_function = [&](Vector& S_rhs, const Vector& S_data, const Real /* time */) { @@ -185,14 +187,9 @@ void main_main () } }; - auto post_update_function = [&](Vector& S_data, const Real /* time */) { - // fill periodic ghost cells - S_data[0].FillBoundary(geom.periodicity()); - }; - - TimeIntegrator > integrator(state, time); + TimeIntegrator> integrator(state, time); + integrator.set_pre_rhs_update(pre_rhs_function); integrator.set_rhs(rhs_function); - integrator.set_post_update(post_update_function); for (int step = 1; step <= nsteps; ++step) { From 46c01bab1535f9d5d2395bc1e5d21d73019b457d Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Sat, 1 Jun 2024 21:27:41 -0700 Subject: [PATCH 05/13] update input file --- .../SUNDIALS/Exec/inputs_sundials_erk | 31 ++++++++++++++----- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk index 3a5a02c5..c03a8e0a 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk @@ -19,14 +19,29 @@ t_final = 1.e-2 integration.type = SUNDIALS # If using the SUNDIALS Submodule, then compile with USE_SUNDIALS=TRUE or -# AMReX_SUNDIALS=ON and set method type here: -# -# ERK = Explicit Runge-Kutta method -# DIRK = Diagonally Implicit Runge-Kutta method -# IMEX = Implicit-Explicit Additive Runge-Kutta method -# MRI = Multirate Infatesimal Method +# AMReX_SUNDIALS=ON + +# 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 integration.sundials.type = ERK -# To select a specific SUNDIALS ERK or DIRK method set integration.sundials.method -# to a method name from https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html +# Select a specific SUNDIALS method by name, see the SUNDIALS documentation for +# the supported method names +# RK Butcher tables: https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html +# MRI Coupling tables: https://sundials.readthedocs.io/en/latest/arkode/Usage/MRIStep_c_interface/MRIStepCoupling.html#arkode-usage-mristep-mristepcoupling-tables # integration.sundials.method = ARKODE_FORWARD_EULER_1_1 + +# For IMEX-RK methods specify the DIRK and ERK method table names +# integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2 +# integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2 + +# If using an an MRI method, then select the fast method type: ERK or DIRK +# integration.sundials.fast_type = ERK + +# Select a specific fast method by name +# integration.sundials.fast_method = ARKODE_FORWARD_EULER_1_1 From 6ceb9082490c9e270c4595cf19ab4b5665bd9eaa Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Mon, 3 Jun 2024 22:09:47 -0700 Subject: [PATCH 06/13] rename update to action --- ExampleCodes/SUNDIALS/Source/main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index c5518cae..eb4d236b 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -188,7 +188,7 @@ void main_main () }; TimeIntegrator> integrator(state, time); - integrator.set_pre_rhs_update(pre_rhs_function); + integrator.set_pre_rhs_action(pre_rhs_function); integrator.set_rhs(rhs_function); for (int step = 1; step <= nsteps; ++step) From 7fa71d2ff570ecaff056e776026227888e9462ad Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Tue, 4 Jun 2024 00:23:18 -0700 Subject: [PATCH 07/13] remove vector of multifabs --- ExampleCodes/SUNDIALS/Source/main.cpp | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index eb4d236b..d3d9ae26 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -116,9 +116,8 @@ void main_main () DistributionMapping dm(ba); // we allocate two phi multifabs; one will store the old state, the other the new. - Vector state; - state.push_back(MultiFab(ba, dm, Ncomp, Nghost)); - auto& phi = state[0]; + MultiFab state(ba, dm, Ncomp, Nghost); + auto& phi = state; // time = starting time in the simulation Real time = 0.0; @@ -156,17 +155,17 @@ void main_main () WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, 0); } - auto pre_rhs_function = [&](Vector& S_data, const Real /* time */) { + auto pre_rhs_function = [&](MultiFab& S_data, const Real /* time */) { // fill periodic ghost cells - S_data[0].FillBoundary(geom.periodicity()); + S_data.FillBoundary(geom.periodicity()); }; - auto rhs_function = [&](Vector& S_rhs, - const Vector& S_data, const Real /* time */) { + auto rhs_function = [&](MultiFab& S_rhs, const MultiFab& S_data, + const Real /* time */) { // loop over boxes - auto& phi_data = S_data[0]; - auto& phi_rhs = S_rhs[0]; + auto& phi_data = S_data; + auto& phi_rhs = S_rhs; for ( MFIter mfi(phi_data); mfi.isValid(); ++mfi ) { const Box& bx = mfi.validbox(); @@ -187,7 +186,7 @@ void main_main () } }; - TimeIntegrator> integrator(state, time); + TimeIntegrator integrator(state, time); integrator.set_pre_rhs_action(pre_rhs_function); integrator.set_rhs(rhs_function); From e04bfe739fc01c5c57dd46171179092808a1ccad Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Tue, 4 Jun 2024 00:34:27 -0700 Subject: [PATCH 08/13] remove duplicate vars --- ExampleCodes/SUNDIALS/Source/main.cpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index d3d9ae26..60c543e1 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -116,8 +116,7 @@ void main_main () DistributionMapping dm(ba); // we allocate two phi multifabs; one will store the old state, the other the new. - MultiFab state(ba, dm, Ncomp, Nghost); - auto& phi = state; + MultiFab phi(ba, dm, Ncomp, Nghost); // time = starting time in the simulation Real time = 0.0; @@ -130,7 +129,7 @@ void main_main () { const Box& bx = mfi.validbox(); - const Array4& phiOld = phi.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) @@ -143,7 +142,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); }); } @@ -186,7 +185,7 @@ void main_main () } }; - TimeIntegrator integrator(state, time); + TimeIntegrator integrator(phi, time); integrator.set_pre_rhs_action(pre_rhs_function); integrator.set_rhs(rhs_function); @@ -196,7 +195,7 @@ void main_main () time += dt_out; // Advance to output time - integrator.evolve(state, time); + integrator.evolve(phi, time); // Tell the I/O Processor to write out which step we're doing amrex::Print() << "Advanced step " << step << "\n"; From 4f66639c3878f422bd94f855818b0a161c1a26fc Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Tue, 4 Jun 2024 13:49:08 -0700 Subject: [PATCH 09/13] update options --- .../SUNDIALS/Exec/inputs_forward_euler | 4 ++-- ExampleCodes/SUNDIALS/Exec/inputs_rk3 | 4 ++-- .../SUNDIALS/Exec/inputs_sundials_erk | 11 +++++----- ExampleCodes/SUNDIALS/Source/main.cpp | 22 ++++++++++--------- 4 files changed, 21 insertions(+), 20 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler b/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler index eeec44c6..746c08b9 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler +++ b/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler @@ -1,8 +1,8 @@ n_cell = 32 max_grid_size = 16 -nsteps = 1000 -plot_int = 100 +nsteps = 10 +plot_int = 1 dt = 1.e-5 diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_rk3 b/ExampleCodes/SUNDIALS/Exec/inputs_rk3 index dded52b9..b363c7cb 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_rk3 +++ b/ExampleCodes/SUNDIALS/Exec/inputs_rk3 @@ -1,8 +1,8 @@ n_cell = 32 max_grid_size = 16 -nsteps = 1000 -plot_int = 100 +nsteps = 10 +plot_int = 1 dt = 1.e-5 diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk index c03a8e0a..5eb6ab52 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk @@ -1,15 +1,14 @@ n_cell = 32 max_grid_size = 16 -# time between outputs -dt_out = 1.e-3 - # number of output steps nsteps = 10 -plot_int = 10 +plot_int = 1 + +dt = 1.e-5 -# final time -t_final = 1.e-2 +# use adaptive time stepping -- should also set integrator tolerances +# adapt_dt = true # INTEGRATION # integration.type can take on the following values: diff --git a/ExampleCodes/SUNDIALS/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index 60c543e1..a04c62ff 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -35,10 +35,10 @@ void main_main () int plot_int; // time step - Real dt_out; + Real dt; - // final time - Real t_final; + // use adaptive time step (dt used for output time) + bool adapt_dt = false; // inputs parameters { @@ -64,10 +64,10 @@ void main_main () pp.query("plot_int",plot_int); // time step - pp.get("dt_out",dt_out); + pp.get("dt",dt); - // time step - pp.get("t_final",t_final); + // use adaptive step sizes + pp.query("adapt_dt",adapt_dt); } // ********************************** @@ -188,11 +188,15 @@ void main_main () TimeIntegrator integrator(phi, time); integrator.set_pre_rhs_action(pre_rhs_function); integrator.set_rhs(rhs_function); + integrator.set_time_step(dt); + if (adapt_dt) { + integrator.set_adaptive_step(); + } for (int step = 1; step <= nsteps; ++step) { - // Set output time - time += dt_out; + // Set time to evolve to + time += dt; // Advance to output time integrator.evolve(phi, time); @@ -206,7 +210,5 @@ void main_main () const std::string& pltfile = amrex::Concatenate("plt",step,5); WriteSingleLevelPlotfile(pltfile, phi, {"phi"}, geom, time, step); } - - if (time > t_final) break; } } From 70c71f8a8bd1274a7f0a2a00b9c84cecb5d7fde2 Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Mon, 17 Jun 2024 11:48:41 -0700 Subject: [PATCH 10/13] simplify input files --- .../SUNDIALS/Exec/inputs_forward_euler | 20 ------------------- ExampleCodes/SUNDIALS/Exec/inputs_rk3 | 10 ---------- .../SUNDIALS/Exec/inputs_sundials_erk | 6 +++--- 3 files changed, 3 insertions(+), 33 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler b/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler index 746c08b9..68bfab5a 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 b363c7cb..4e38700f 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_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk index 5eb6ab52..d48d36e4 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk @@ -15,10 +15,10 @@ dt = 1.e-5 # 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 = SUNDIALS - +# # 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 @@ -33,7 +33,7 @@ integration.sundials.type = ERK # the supported method names # RK Butcher tables: https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html # MRI Coupling tables: https://sundials.readthedocs.io/en/latest/arkode/Usage/MRIStep_c_interface/MRIStepCoupling.html#arkode-usage-mristep-mristepcoupling-tables -# integration.sundials.method = ARKODE_FORWARD_EULER_1_1 +integration.sundials.method = ARKODE_FORWARD_EULER_1_1 # For IMEX-RK methods specify the DIRK and ERK method table names # integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2 From 9439b54451f8b771b0f8433366fcbbeaf3fd94c7 Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Wed, 19 Jun 2024 09:37:34 -0700 Subject: [PATCH 11/13] revert input files nsteps and plot_int changes --- ExampleCodes/SUNDIALS/Exec/inputs_forward_euler | 4 ++-- ExampleCodes/SUNDIALS/Exec/inputs_rk3 | 4 ++-- ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler b/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler index 68bfab5a..39b09830 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler +++ b/ExampleCodes/SUNDIALS/Exec/inputs_forward_euler @@ -1,8 +1,8 @@ n_cell = 32 max_grid_size = 16 -nsteps = 10 -plot_int = 1 +nsteps = 1000 +plot_int = 100 dt = 1.e-5 diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_rk3 b/ExampleCodes/SUNDIALS/Exec/inputs_rk3 index 4e38700f..6eb60137 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_rk3 +++ b/ExampleCodes/SUNDIALS/Exec/inputs_rk3 @@ -1,8 +1,8 @@ n_cell = 32 max_grid_size = 16 -nsteps = 10 -plot_int = 1 +nsteps = 1000 +plot_int = 100 dt = 1.e-5 diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk index d48d36e4..e9df5617 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk @@ -2,8 +2,8 @@ n_cell = 32 max_grid_size = 16 # number of output steps -nsteps = 10 -plot_int = 1 +nsteps = 1000 +plot_int = 100 dt = 1.e-5 From 2c4ee605aee7640f1c01852d8f0f90de12d4d0cc Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Wed, 19 Jun 2024 09:43:26 -0700 Subject: [PATCH 12/13] remove comment --- ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk | 1 - 1 file changed, 1 deletion(-) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk index e9df5617..bd3a7203 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk @@ -1,7 +1,6 @@ n_cell = 32 max_grid_size = 16 -# number of output steps nsteps = 1000 plot_int = 100 From ac95881a42449bb6c942c87c30952a44b5523979 Mon Sep 17 00:00:00 2001 From: "David J. Gardner" Date: Wed, 19 Jun 2024 10:15:48 -0700 Subject: [PATCH 13/13] updates for adaptive stepping --- .../{inputs_sundials_erk => inputs_sundials} | 29 +++++++++---------- ExampleCodes/SUNDIALS/Source/main.cpp | 14 +++++++-- 2 files changed, 26 insertions(+), 17 deletions(-) rename ExampleCodes/SUNDIALS/Exec/{inputs_sundials_erk => inputs_sundials} (50%) diff --git a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk b/ExampleCodes/SUNDIALS/Exec/inputs_sundials similarity index 50% rename from ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk rename to ExampleCodes/SUNDIALS/Exec/inputs_sundials index bd3a7203..bbdf2c69 100644 --- a/ExampleCodes/SUNDIALS/Exec/inputs_sundials_erk +++ b/ExampleCodes/SUNDIALS/Exec/inputs_sundials @@ -6,8 +6,10 @@ plot_int = 100 dt = 1.e-5 -# use adaptive time stepping -- should also set integrator tolerances +# 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: @@ -26,20 +28,17 @@ integration.type = SUNDIALS # EX-MRI = Explicit Multirate Infatesimal method # IM-MRI = Implicit Multirate Infatesimal method # IMEX-MRI = Implicit-Explicit Multirate Infatesimal method -integration.sundials.type = ERK - -# Select a specific SUNDIALS method by name, see the SUNDIALS documentation for -# the supported method names -# RK Butcher tables: https://sundials.readthedocs.io/en/latest/arkode/Butcher_link.html -# MRI Coupling tables: https://sundials.readthedocs.io/en/latest/arkode/Usage/MRIStep_c_interface/MRIStepCoupling.html#arkode-usage-mristep-mristepcoupling-tables -integration.sundials.method = ARKODE_FORWARD_EULER_1_1 +# +# Optionally select a specific SUNDIALS method by name, see the SUNDIALS +# documentation for the supported method names -# For IMEX-RK methods specify the DIRK and ERK method table names -# integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2 -# integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2 +# Use the SUNDIALS default method for the chosen type (fixed or adaptive step sizes) +integration.sundials.type = ERK -# If using an an MRI method, then select the fast method type: ERK or DIRK -# integration.sundials.fast_type = ERK +# Use forward Euler (fixed step sizes only) +# integration.sundials.type = ERK +# integration.sundials.method = ARKODE_FORWARD_EULER_1_1 -# Select a specific fast method by name -# integration.sundials.fast_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/Source/main.cpp b/ExampleCodes/SUNDIALS/Source/main.cpp index a04c62ff..cd24a35f 100644 --- a/ExampleCodes/SUNDIALS/Source/main.cpp +++ b/ExampleCodes/SUNDIALS/Source/main.cpp @@ -37,9 +37,13 @@ void main_main () // time step Real dt; - // use adaptive time step (dt used for output time) + // 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 @@ -68,6 +72,10 @@ void main_main () // use adaptive step sizes pp.query("adapt_dt",adapt_dt); + + // adaptive step tolerances + pp.query("reltol",reltol); + pp.query("abstol",abstol); } // ********************************** @@ -188,9 +196,11 @@ void main_main () TimeIntegrator integrator(phi, time); integrator.set_pre_rhs_action(pre_rhs_function); integrator.set_rhs(rhs_function); - integrator.set_time_step(dt); 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)