Skip to content

Commit

Permalink
make Omega a persistent variable and update project_initial_velocity …
Browse files Browse the repository at this point in the history 
…flag (#1940)
  • Loading branch information
@asalmgren
asalmgren authored Nov 8, 2024
1 parent c0983e7 commit f531433
Show file tree
Hide file tree
Showing 16 changed files with 1,054 additions and 64 deletions.
74 changes: 39 additions & 35 deletions Docs/sphinx_doc/Inputs.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1165,41 +1165,41 @@ In addition, there is a run-time option to project the initial velocity field to
List of Parameters
------------------

+----------------------------------+-------------------+--------------------+------------------+
| Parameter | Definition | Acceptable | Default |
| | | Values | |
+==================================+===================+====================+==================+
| **erf.init_type** | Initialization | “custom”, | “*custom*” |
| | type | “ideal”, | |
| | | "real", | |
| | |"input_sounding" | |
+----------------------------------+-------------------+--------------------+------------------+
| **erf.input_sounding_file** | Path to WRF-style | String | "input_sounding" |
| | input sounding | | |
| | file | | |
+----------------------------------+-------------------+--------------------+------------------+
| **erf.init_sounding_ideal** | Perform | true or false | false |
| | initialization | | |
| | like WRF's | | |
| | ideal.exe | | |
+----------------------------------+-------------------+--------------------+------------------+
| **erf.use_real_bcs** | If init_type is | true or false | true if |
| | real or metgrid, | | if init_type |
| | do we want to use | | is real or |
| | these bcs? | | metgrid; |
| | | | else false |
+----------------------------------+-------------------+--------------------+------------------+
| **erf.nc_init_file** | NetCDF file with | String | NONE |
| | initial mesoscale | | |
| | data | | |
+----------------------------------+-------------------+--------------------+------------------+
| **erf.nc_bdy_file** | NetCDF file with | String | NONE |
| | mesoscale data at | | |
| | lateral boundaries| | |
+----------------------------------+-------------------+--------------------+------------------+
| **erf.project_initial_velocity** | project initial | Integer | 1 |
| | velocity? | | |
+----------------------------------+-------------------+--------------------+------------------+
+----------------------------------+-------------------+--------------------+-----------------------+
| Parameter | Definition | Acceptable | Default |
| | | Values | |
+==================================+===================+====================+=======================+
| **erf.init_type** | Initialization | “custom”, | “*custom*” |
| | type | “ideal”, | |
| | | "real", | |
| | |"input_sounding" | |
+----------------------------------+-------------------+--------------------+-----------------------+
| **erf.input_sounding_file** | Path to WRF-style | String | "input_sounding" |
| | input sounding | | |
| | file | | |
+----------------------------------+-------------------+--------------------+-----------------------+
| **erf.init_sounding_ideal** | Perform | true or false | false |
| | initialization | | |
| | like WRF's | | |
| | ideal.exe | | |
+----------------------------------+-------------------+--------------------+-----------------------+
| **erf.use_real_bcs** | If init_type is | true or false | true if |
| | real or metgrid, | | if init_type |
| | do we want to use | | is real or |
| | these bcs? | | metgrid; |
| | | | else false |
+----------------------------------+-------------------+--------------------+-----------------------+
| **erf.nc_init_file** | NetCDF file with | String | NONE |
| | initial mesoscale | | |
| | data | | |
+----------------------------------+-------------------+--------------------+-----------------------+
| **erf.nc_bdy_file** | NetCDF file with | String | NONE |
| | mesoscale data at | | |
| | lateral boundaries| | |
+----------------------------------+-------------------+--------------------+-----------------------+
| **erf.project_initial_velocity** | project initial | true or false | true if anelastic; |
| | velocity? | | false if compressible |
+----------------------------------+-------------------+--------------------+-----------------------+

Notes
-----------------
Expand Down Expand Up @@ -1227,6 +1227,10 @@ integrating the hydrostatic equation from the surface.

If **erf.init_type = custom** or **erf.init_type = input_sounding**, ``erf.nc_init_file`` and ``erf.nc_bdy_file`` do not need to be set.

Note that the **erf.project_initial_velocity** option is available for all **init_type** options. If using the anelastic
formulation this will be true regardless of the input; if using the compressible formulation the default is false but
that can be over-written.

Map Scale Factors
=================

Expand Down
9 changes: 4 additions & 5 deletions Source/DataStructs/ERF_DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -193,11 +193,6 @@ struct SolverChoice {
}
// *******************************************************************************

// Should we project the initial velocity field to make it divergence-free?
pp.query("project_initial_velocity", project_initial_velocity);

// *******************************************************************************

bool any_anelastic = false;
for (int i = 0; i <= max_level; ++i) {
if (anelastic[i] == 1) any_anelastic = true;
Expand All @@ -209,10 +204,14 @@ struct SolverChoice {
project_initial_velocity = true;
buoyancy_type = 2; // uses Tprime
} else {
pp.query("project_initial_velocity", project_initial_velocity);

constant_density = false; // We default to false but allow the user to set it
pp.query("constant_density", constant_density);
}

// *******************************************************************************

pp.query("project_every_stage", project_every_stage);
pp.query("ncorr", ncorr);
pp.query("poisson_abstol", poisson_abstol);
Expand Down
9 changes: 7 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,10 +164,12 @@ public:
#endif

// Project the velocities to be divergence-free -- this is only relevant if anelastic == 1
void project_velocities (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars, amrex::MultiFab& p);
void project_velocities (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars,
amrex::MultiFab& Omega, amrex::MultiFab& p);

// Project the velocities to be divergence-free with a thin body
void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars, amrex::MultiFab& p);
void project_velocities_tb (int lev, amrex::Real dt, amrex::Vector<amrex::MultiFab >& vars,
amrex::MultiFab& Omega, amrex::MultiFab& p);

// Define the projection bc's based on the domain bc types
amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM>
Expand Down Expand Up @@ -731,6 +733,9 @@ private:
amrex::Vector<amrex::MultiFab> rW_old;
amrex::Vector<amrex::MultiFab> rW_new;

// Normal momentum on z-face
amrex::Vector<amrex::MultiFab> Omega;

// amrex::Vector<amrex::MultiFab> xmom_crse_rhs;
// amrex::Vector<amrex::MultiFab> ymom_crse_rhs;
amrex::Vector<amrex::MultiFab> zmom_crse_rhs;
Expand Down
9 changes: 8 additions & 1 deletion Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -197,6 +197,8 @@ ERF::ERF_shared ()
rV_old.resize(nlevs_max);
rW_old.resize(nlevs_max);

Omega.resize(nlevs_max);

// xmom_crse_rhs.resize(nlevs_max);
// ymom_crse_rhs.resize(nlevs_max);
zmom_crse_rhs.resize(nlevs_max);
Expand Down Expand Up @@ -856,6 +858,11 @@ ERF::InitData_post ()
}
}

//
// If we are starting from scratch, we have the option to project the initial velocity field
// regardless of how we initialized.
// pp_inc is used as scratch space here; we zero it out after the projection
//
if (restart_chkfile == "")
{
// Note -- this projection is only defined for no terrain
Expand All @@ -864,7 +871,7 @@ ERF::InitData_post ()
Real dummy_dt = 1.0;
for (int lev = 0; lev <= finest_level; ++lev)
{
project_velocities(lev, dummy_dt, vars_new[lev], pp_inc[lev]);
project_velocities(lev, dummy_dt, vars_new[lev], Omega[lev], pp_inc[lev]);
pp_inc[lev].setVal(0.);
}
}
Expand Down
5 changes: 5 additions & 0 deletions Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,11 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
pp_inc[lev].setVal(0.0);
}

// We keep Omega as a persistent variable now in order to use the projected Omega
// directly rather than calculating it again from (u,v,w)
Omega[lev].define(convert(ba, IntVect(0,0,1)), dm, 1, 1);
Omega[lev].setVal(5.6e23);

// ********************************************************************************************
// These are just used for scratch in the time integrator but we might as well define them here
// ********************************************************************************************
Expand Down
2 changes: 2 additions & 0 deletions Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -448,6 +448,8 @@ ERF::ClearLevel (int lev)
rW_new[lev].clear();
rW_old[lev].clear();

Omega[lev].clear();

if (lev > 0) {
zmom_crse_rhs[lev].clear();
}
Expand Down
4 changes: 2 additions & 2 deletions Source/TimeIntegration/ERF_TI_no_substep_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,9 +138,9 @@
bool have_tb = (thin_xforce[0] || thin_yforce[0] || thin_zforce[0]);
if (solverChoice.project_every_stage || (nrk==1)) {
if (have_tb) {
project_velocities_tb(level, slow_dt, S_sum, pp_inc[level]);
project_velocities_tb(level, slow_dt, S_sum, Omega[level], pp_inc[level]);
} else {
project_velocities(level, slow_dt, S_sum, pp_inc[level]);
project_velocities(level, slow_dt, S_sum, Omega[level], pp_inc[level]);
}
}
}
Expand Down
6 changes: 3 additions & 3 deletions Source/TimeIntegration/ERF_TI_slow_rhs_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@

erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch,
xvel_new, yvel_new, zvel_new,
z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src,
z_t_rk[level], Omega[level], cc_src, xmom_src, ymom_src, zmom_src,
(level > 0) ? &zmom_crse_rhs[level] : nullptr,
Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(),
Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(),
Expand Down Expand Up @@ -242,7 +242,7 @@

erf_slow_rhs_pre(level, finest_level, nrk, slow_dt, S_rhs, S_old, S_data, S_prim, S_scratch,
xvel_new, yvel_new, zvel_new,
z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src,
z_t_rk[level], Omega[level], cc_src, xmom_src, ymom_src, zmom_src,
(level > 0) ? &zmom_crse_rhs[level] : nullptr,
Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(),
Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(),
Expand Down Expand Up @@ -448,7 +448,7 @@
erf_slow_rhs_pre(level, finest_level, nrk, slow_dt,
S_rhs, S_old, S_data, S_prim, S_scratch,
xvel_new, yvel_new, zvel_new,
z_t_rk[level], Omega, cc_src, xmom_src, ymom_src, zmom_src,
z_t_rk[level], Omega[level], cc_src, xmom_src, ymom_src, zmom_src,
(level > 0) ? &zmom_crse_rhs[level] : nullptr,
Tau11_lev[level].get(), Tau22_lev[level].get(), Tau33_lev[level].get(), Tau12_lev[level].get(),
Tau13_lev[level].get(), Tau21_lev[level].get(), Tau23_lev[level].get(), Tau31_lev[level].get(),
Expand Down
8 changes: 4 additions & 4 deletions Source/TimeIntegration/ERF_TI_substep_fun.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
S_slow_rhs, S_old, S_stage, S_prim, pi_stage, fast_coeffs,
S_data, S_scratch, fine_geom,
solverChoice.gravity, solverChoice.use_lagged_delta_rt,
Omega, z_t_rk[level], z_t_pert.get(),
Omega[level], z_t_rk[level], z_t_pert.get(),
z_phys_nd[level], z_phys_nd_new[level], z_phys_nd_src[level],
detJ_cc[level], detJ_cc_new[level], detJ_cc_src[level],
dtau, beta_s, inv_fac,
Expand All @@ -114,7 +114,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs,
S_data, S_scratch, fine_geom,
solverChoice.gravity, solverChoice.use_lagged_delta_rt,
Omega, z_t_rk[level], z_t_pert.get(),
Omega[level], z_t_rk[level], z_t_pert.get(),
z_phys_nd[level], z_phys_nd_new[level], z_phys_nd_src[level],
detJ_cc[level], detJ_cc_new[level], detJ_cc_src[level],
dtau, beta_s, inv_fac,
Expand All @@ -133,7 +133,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
// and S_data is the new-time solution to be defined here
erf_fast_rhs_T(fast_step, nrk, level, finest_level,
S_slow_rhs, S_old, S_stage, S_prim, pi_stage, fast_coeffs,
S_data, S_scratch, fine_geom, solverChoice.gravity, Omega,
S_data, S_scratch, fine_geom, solverChoice.gravity, Omega[level],
z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac,
mapfac_m[level], mapfac_u[level], mapfac_v[level],
fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
Expand All @@ -142,7 +142,7 @@ auto fast_rhs_fun = [&](int fast_step, int /*n_sub*/, int nrk,
// and as the new-time solution to be defined here
erf_fast_rhs_T(fast_step, nrk, level, finest_level,
S_slow_rhs, S_data, S_stage, S_prim, pi_stage, fast_coeffs,
S_data, S_scratch, fine_geom, solverChoice.gravity, Omega,
S_data, S_scratch, fine_geom, solverChoice.gravity, Omega[level],
z_phys_nd[level], detJ_cc[level], dtau, beta_s, inv_fac,
mapfac_m[level], mapfac_u[level], mapfac_v[level],
fr_as_crse, fr_as_fine, l_use_moisture, l_reflux, l_implicit_substepping);
Expand Down
2 changes: 0 additions & 2 deletions Source/TimeIntegration/ERF_advance_dycore.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,8 +183,6 @@ void ERF::advance_dycore(int level,
} // l_use_diff
} // profile

MultiFab Omega (state_old[IntVars::zmom].boxArray(),dm,1,1);

#include "ERF_TI_utils.H"

// Additional SFS quantities, calculated once per timestep
Expand Down
Loading

0 comments on commit f531433

Please sign in to comment.