diff --git a/integration/BackwardEuler/Make.package b/integration/BackwardEuler/Make.package index 3f96bef756..ca0dd265d2 100644 --- a/integration/BackwardEuler/Make.package +++ b/integration/BackwardEuler/Make.package @@ -1,4 +1,10 @@ -CEXE_headers += actual_integrator.H +ifeq ($(USE_SIMPLIFIED_SDC), TRUE) + CEXE_headers += actual_integrator_simplified_sdc.H +else + ifneq ($(USE_TRUE_SDC), TRUE) + CEXE_headers += actual_integrator.H + endif +endif + CEXE_headers += be_integrator.H CEXE_headers += be_type.H -CEXE_headers += be_type_strang.H diff --git a/integration/BackwardEuler/actual_integrator_simplified_sdc.H b/integration/BackwardEuler/actual_integrator_simplified_sdc.H new file mode 100644 index 0000000000..7d34ea837e --- /dev/null +++ b/integration/BackwardEuler/actual_integrator_simplified_sdc.H @@ -0,0 +1,181 @@ +#ifndef actual_integrator_H +#define actual_integrator_H + +#include + +#include +#include +#include + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void actual_integrator (burn_t& state, const Real dt) +{ + + be_t be; + + // set the Jacobian type + be.jacobian_type = jacobian; + + // Start off by assuming a successful burn. + + state.success = true; + + // Initialize the integration time. + + be.t = 0.0; + be.tout = dt; + + // Fill in the initial integration state. + + burn_to_int(state, be); + + // Save the initial composition and temperature for our later diagnostics. + +#ifndef AMREX_USE_GPU + Real xn_in[NumSpec]; + for (int n = 0; n < NumSpec; ++n) { + xn_in[n] = state.y[SFS+n] / state.y[SRHO]; + } + // we are assuming that the temperature was valid on input + Real T_in = state.T; +#ifdef AUX_THERMO + Real aux_in[NumAux]; + for (int n = 0; n < NumAux; ++n) { + aux_in[n] = state.y[SFX+n] / state.y[SRHO]; + } +#endif + Real rhoe_in = state.y[SEINT]; +#endif + + + // Set the tolerances. + + Real sdc_tol_fac = std::pow(sdc_burn_tol_factor, state.num_sdc_iters - state.sdc_iter - 1); + + // we use 1-based indexing inside of BackwardEuler, so we need to shift the + // indices SRHO, SFS, etc by 1 + + Real sdc_min_density = amrex::min(state.rho, state.rho_orig + state.ydot_a[SRHO] * dt); + + be.atol_enuc = sdc_min_density * atol_enuc * sdc_tol_fac; + be.rtol_enuc = rtol_enuc * sdc_tol_fac; + + // Note: we define the input atol for species to refer only to the + // mass fraction part, and we multiply by a representative density + // so that atol becomes an absolutely tolerance on (rho X) + + be.atol_spec = sdc_min_density * atol_spec * sdc_tol_fac; + be.rtol_spec = rtol_spec * sdc_tol_fac; + + // Call the integration routine. + + int istate = be_integrator(state, be); + + // Get the number of RHS and Jacobian evaluations. + + state.n_rhs = be.n_rhs; + state.n_jac = be.n_jac; + state.n_step = be.n_step; + + // Copy the integration data back to the burn state. + // This will also update the aux state from X if we are using NSE + + int_to_burn(be.t, be, state); + + // we only evolved (rho e), not (rho E), so we need to update the + // total energy now to ensure we are conservative + + Real rho_Sdot = 0.0_rt; + if (state.time > 0) { + rho_Sdot = (be.y(SEINT+1) - state.rhoe_orig) / state.time - state.ydot_a[SEINT]; + } + + state.y[SEDEN] += state.time * (state.ydot_a[SEDEN] + rho_Sdot); + + // also momentum + + state.y[SMX] += state.time * state.ydot_a[SMX]; + state.y[SMY] += state.time * state.ydot_a[SMY]; + state.y[SMZ] += state.time * state.ydot_a[SMZ]; + + // normalize the abundances on exit. We'll assume that the driver + // calling this is making use of the conserved state (state.y[]), + // so that is what will be normalized. + + normalize_abundances_sdc_burn(state); + + if (istate < 0) { + state.success = false; + } + + +#ifndef AMREX_USE_GPU + if (burner_verbose) { + // Print out some integration statistics, if desired. + std::cout << "integration summary: " << std::endl; + std::cout << "dens: " << state.rho << " temp: " << state.T << std::endl; + std::cout << "energy released: " << state.e << std::endl; + std::cout << "number of steps taken: " << state.n_step << std::endl; + std::cout << "number of f evaluations: " << state.n_rhs << std::endl; + } +#endif + + // If we failed, print out the current state of the integration. + + if (!state.success) { +#ifndef AMREX_USE_GPU + std::cout << Font::Bold << FGColor::Red << "[ERROR] integration failed in net" << ResetDisplay << std::endl; + std::cout << "istate = " << istate << std::endl; + std::cout << "zone = (" << state.i << ", " << state.j << ", " << state.k << ")" << std::endl; + std::cout << "time = " << be.t << std::endl; + std::cout << "dt = " << std::setprecision(16) << dt << std::endl; + std::cout << "dens start = " << std::setprecision(16) << state.rho_orig << std::endl; + std::cout << "temp start = " << std::setprecision(16) << T_in << std::endl; + std::cout << "rhoe start = " << std::setprecision(16) << rhoe_in << std::endl; + std::cout << "xn start = "; + for (int n = 0; n < NumSpec; ++n) { + std::cout << std::setprecision(16) << xn_in[n] << " "; + } + std::cout << std::endl; +#ifdef AUX_THERMO + std::cout << "aux start = "; + for (int n = 0; n < NumAux; ++n) { + std::cout << std::setprecision(16) << aux_in[n] << " "; + } + std::cout << std::endl; +#endif + std::cout << "dens current = " << std::setprecision(16) << state.rho << std::endl; + std::cout << "temp current = " << std::setprecision(16) << state.T << std::endl; + std::cout << "xn current = "; + for (int n = 0; n < NumSpec; ++n) { + std::cout << std::setprecision(16) << state.xn[n] << " "; + } + std::cout << std::endl; +#ifdef AUX_THERMO + std::cout << "aux current = "; + for (int n = 0; n < NumAux; ++n) { + std::cout << std::setprecision(16) << state.aux[n] << " "; + } + std::cout << std::endl; +#endif + std::cout << "A(rho) = " << std::setprecision(16) << state.ydot_a[SRHO] << std::endl; + std::cout << "A(rho e) = " << std::setprecision(16) << state.ydot_a[SEINT] << std::endl; + std::cout << "A(rho X_k) = "; + for (int n = 0; n < NumSpec; n++) { + std::cout << std::setprecision(16) << state.ydot_a[SFS+n] << " "; + } + std::cout << std::endl; +#ifdef AUX_THERMO + std::cout << "A(rho aux_k) = "; + for (int n = 0; n < NumAux; n++) { + std::cout << std::setprecision(16) << state.ydot_a[SFX+n] << " "; + } + std::cout << std::endl; +#endif +#endif + } + + +} + +#endif diff --git a/integration/BackwardEuler/be_integrator.H b/integration/BackwardEuler/be_integrator.H index 4ec50629d8..d98edb4510 100644 --- a/integration/BackwardEuler/be_integrator.H +++ b/integration/BackwardEuler/be_integrator.H @@ -2,14 +2,18 @@ #define BE_INTEGRATOR_H #include -#include #include #include #include #include #include #include - +#ifdef STRANG +#include +#endif +#ifdef SIMPLIFIED_SDC +#include +#endif AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real initial_dt (burn_t& state, be_t& be, @@ -48,9 +52,9 @@ Real initial_dt (burn_t& state, be_t& be, // routine for (int ii = 1; ii <= NumSpec; ii++) { - ewt(ii) = rtol_spec * std::abs(y_old(ii)) + atol_spec; + ewt(ii) = be.rtol_spec * std::abs(y_old(ii)) + be.atol_spec; } - ewt(net_ienuc) = rtol_enuc * std::abs(y_old(net_ienuc)) + atol_enuc; + ewt(net_ienuc) = be.rtol_enuc * std::abs(y_old(net_ienuc)) + be.atol_enuc; // Construct the trial point. @@ -143,13 +147,15 @@ int single_step (burn_t& state, be_t& be, const Real dt) if (be.jacobian_type == 1) { jac(be.t, state, be, be.jac); - } else { jac_info_t jac_info; jac_info.h = dt; numerical_jac(state, jac_info, be.jac); + be.n_rhs += (NumSpec+1); } + be.n_jac++; + if (!integrate_energy) { for (int j = 1; j <= BE_NEQS; j++) { be.jac(net_ienuc, j) = 0.0_rt; @@ -242,6 +248,7 @@ int be_integrator (burn_t& state, be_t& be) { be.n_rhs = 0; + be.n_jac = 0; be.n_step = 0; int ierr; @@ -304,9 +311,9 @@ int be_integrator (burn_t& state, be_t& be) Array1D w; for (int n = 1; n <= NumSpec; n++) { - w(n) = 1.0_rt / (rtol_spec * std::abs(y_fine(n)) + atol_spec); + w(n) = 1.0_rt / (be.rtol_spec * std::abs(y_fine(n)) + be.atol_spec); } - w(net_ienuc) = 1.0_rt / (rtol_enuc * std::abs(y_fine(net_ienuc)) + atol_enuc); + w(net_ienuc) = 1.0_rt / (be.rtol_enuc * std::abs(y_fine(net_ienuc)) + be.atol_enuc); // now look for w |y_fine - y_coarse| < 1 diff --git a/integration/BackwardEuler/be_type.H b/integration/BackwardEuler/be_type.H index 994c45a6d5..129d8a2b05 100644 --- a/integration/BackwardEuler/be_type.H +++ b/integration/BackwardEuler/be_type.H @@ -9,7 +9,9 @@ #include #ifdef STRANG #include -#include +#endif +#ifdef SIMPLIFIED_SDC +#include #endif #include @@ -50,7 +52,6 @@ struct be_t { JacNetArray2D jac; short jacobian_type; - }; #endif diff --git a/integration/Make.package b/integration/Make.package index fe34b19b3c..920367dbc0 100644 --- a/integration/Make.package +++ b/integration/Make.package @@ -14,13 +14,15 @@ endif CEXE_headers += integrator.H CEXE_headers += integrator_data.H -ifneq ($(USE_TRUE_SDC), TRUE) - CEXE_headers += integrator_type_strang.H - CEXE_headers += integrator_rhs_strang.H -endif + ifeq ($(USE_SIMPLIFIED_SDC), TRUE) - CEXE_headers += integrator_type_simplified_sdc.H CEXE_headers += integrator_rhs_simplified_sdc.H + CEXE_headers += integrator_type_simplified_sdc.H +else + ifneq ($(USE_TRUE_SDC), TRUE) + CEXE_headers += integrator_type_strang.H + CEXE_headers += integrator_rhs_strang.H + endif endif INCLUDE_LOCATIONS += $(MICROPHYSICS_HOME)/integration/utils diff --git a/integration/VODE/_parameters b/integration/VODE/_parameters index e7f6dd5b33..7bfebf611a 100644 --- a/integration/VODE/_parameters +++ b/integration/VODE/_parameters @@ -1,8 +1,5 @@ @namespace: integrator -# SDC iteration tolerance adjustment factor -sdc_burn_tol_factor real 1.d0 - # for the step rejection logic on mass fractions, we only consider # species that are > X_reject_buffer * atol_spec X_reject_buffer real 1.0 diff --git a/integration/_parameters b/integration/_parameters index 3be0be2d1b..4080c6da4c 100644 --- a/integration/_parameters +++ b/integration/_parameters @@ -72,3 +72,7 @@ use_number_densities logical .false. # flag for tuning on the subtraction of internal energy subtract_internal_energy logical .true. + +# SDC iteration tolerance adjustment factor +sdc_burn_tol_factor real 1.d0 + diff --git a/integration/integrator_rhs_simplified_sdc.H b/integration/integrator_rhs_simplified_sdc.H index 664275fe52..f58d502cb7 100644 --- a/integration/integrator_rhs_simplified_sdc.H +++ b/integration/integrator_rhs_simplified_sdc.H @@ -16,7 +16,7 @@ #include #endif -// The f_rhs routine provides the right-hand-side for the integrator. +// The f_rhs routine provides the right-hand-side for the integration solver. // This is a generic interface that calls the specific RHS routine in the // network you're actually using. diff --git a/integration/integrator_type_simplified_sdc.H b/integration/integrator_type_simplified_sdc.H index 412fadd32c..5a58efab81 100644 --- a/integration/integrator_type_simplified_sdc.H +++ b/integration/integrator_type_simplified_sdc.H @@ -95,7 +95,7 @@ void clean_state(const Real time, burn_t& state, T& int_state) for (int n = 1; n <= NumSpec; ++n) { // we use 1-based indexing, so we need to offset SFS int_state.y(SFS+n) = amrex::max(amrex::min(int_state.y(SFS+n), state.rho), - state.rho * SMALL_X_SAFE); + state.rho * SMALL_X_SAFE); } }