Skip to content

Commit

Permalink
sets operator split option for reactive+transport runs
Browse files Browse the repository at this point in the history
  • Loading branch information
@ecoon
ecoon committed Mar 4, 2025
1 parent 72c7b7e commit f37124b
Show file tree
Hide file tree
Showing 8 changed files with 153 additions and 29 deletions.
2 changes: 1 addition & 1 deletion src/pks/chem_pk_helpers.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ convertConcentrationToATS(const Epetra_MultiVector& mol_dens,
Epetra_MultiVector& tcc_ats)
{
// convert from [mol C / L] to mol fraction [mol C / mol H20]
tcc_ats.ReciprocalMultiply(1.e3, tcc_amanzi, mol_dens, 0.);
tcc_ats.ReciprocalMultiply(1.e3, mol_dens, tcc_amanzi, 0.);
}


Expand Down
4 changes: 4 additions & 0 deletions src/pks/flow/overland_pressure_pk.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,10 @@ OverlandPressureFlow::parseParameterList()
patm_hard_limit_ = plist_->get<bool>("allow no negative ponded depths", false);
min_vel_ponded_depth_ = plist_->get<double>("min ponded depth for velocity calculation", 1e-2);
min_tidal_bc_ponded_depth_ = plist_->get<double>("min ponded depth for tidal bc", 0.02);

// require a few primary variable keys now to set the leaf node in the dep graph
requireAtCurrent(pd_key_, tag_current_, *S_, name_);

}


Expand Down
3 changes: 3 additions & 0 deletions src/pks/flow/richards_pk.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,9 @@ Richards::parseParameterList()
ss_primary_key_ = Keys::readKey(*plist_, domain_surf, "pressure", "pressure");
}

// require a few primary variable keys now to set the leaf node in the dep graph
requireAtCurrent(sat_key_, tag_current_, *S_, name_);

// parse inherited lists
PK_PhysicalBDF_Default::parseParameterList();
}
Expand Down
12 changes: 10 additions & 2 deletions src/pks/mpc/mpc_coupled_reactivetransport.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,13 +52,21 @@ MPCCoupledReactiveTransport::parseParameterList()
"primary variable key", "total_component_concentration");

// chemistry and transport share the same primary variable
pks_list_->sublist(chem_names[0]).set<std::string>("primary variable key", tcc_surf_key_);
pks_list_->sublist(chem_names[1]).set<std::string>("primary variable key", tcc_key_);
pks_list_->sublist(chem_names[0]).set<std::string>("primary variable key", tcc_key_);
pks_list_->sublist(chem_names[1]).set<std::string>("primary variable key", tcc_surf_key_);
pks_list_->sublist(chem_names[0]).set<std::string>("primary variable password", name_);
pks_list_->sublist(chem_names[1]).set<std::string>("primary variable password", name_);
pks_list_->sublist(transport_names[0]).set<std::string>("primary variable password", name_);
pks_list_->sublist(transport_names[1]).set<std::string>("primary variable password", name_);

// tell chemistry to operator split
pks_list_->sublist(chem_names[0]).set("operator split", true);
pks_list_->sublist(chem_names[1]).set("operator split", true);

// Only reaction PKs set IC, but all need the list to be PK_Physical.
pks_list_->sublist(transport_names[0]).sublist("initial conditions");
pks_list_->sublist(transport_names[1]).sublist("initial conditions");

cast_sub_pks_();
coupled_chemistry_pk_->parseParameterList();

Expand Down
146 changes: 122 additions & 24 deletions src/pks/mpc/mpc_flow_transport.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ MPCFlowTransport::parseParameterList()
// are we doing surface, subsurface, or integrated transport?
auto flow_plist = getSubPKPlist_(0);
bool surface(false), subsurface(false);
bool chemistry(false);
Key surf_lwc_key, sub_lwc_key;

if (flow_plist->isParameter("domain name")) {
Expand All @@ -40,7 +41,9 @@ MPCFlowTransport::parseParameterList()
AMANZI_ASSERT(transport_names.size() == 2);

if (!pks_list_->sublist(transport_names[0]).isParameter("domain name")) {
// reactions first
chemistry = true;

// chemistry first
transport_names = pks_list_->sublist(transport_names[1]).get<Teuchos::Array<std::string>>("PKs order");
}

Expand All @@ -64,25 +67,58 @@ MPCFlowTransport::parseParameterList()
// hard-code this as the default. If this breaks in the future it can be
// fixed. --ETC
Teuchos::ParameterList& flux_list = S_->GetEvaluatorList(Keys::getKey("water_flux", transport_next_tag));
flux_list.set<std::string>("evaluator type", "alias");
flux_list.set<std::string>("target", Keys::getKey("water_flux", flow_next_tag, true));
if (!flux_list.isParameter("evaluator type")) {
flux_list.set<std::string>("evaluator type", "alias");
flux_list.set<std::string>("target", Keys::getKey("water_flux", flow_next_tag, true));
}

// velocity for dispersivity
Teuchos::ParameterList& velo_list = S_->GetEvaluatorList(Keys::getKey("darcy_velocity", transport_next_tag));
velo_list.set<std::string>("evaluator type", "alias");
velo_list.set<std::string>("target", Keys::getKey("darcy_velocity", flow_next_tag, true));
if (!velo_list.isParameter("evaluator type")) {
velo_list.set<std::string>("evaluator type", "alias");
velo_list.set<std::string>("target", Keys::getKey("darcy_velocity", flow_next_tag, true));
}

// now set the liquid water content as an interpolated field at next
// note that current copy is kept by transport PK
Teuchos::ParameterList& lwc_list_next = S_->GetEvaluatorList(Keys::getKey(sub_lwc_key, transport_next_tag));
lwc_list_next.set<std::string>("evaluator type", "temporal interpolation");
lwc_list_next.set<std::string>("current tag", flow_current_tag.get());
lwc_list_next.set<std::string>("next tag", flow_next_tag.get());
if (!lwc_list_next.isParameter("evaluator type")) {
lwc_list_next.set<std::string>("evaluator type", "temporal interpolation");
lwc_list_next.set<std::string>("current tag", flow_current_tag.get());
lwc_list_next.set<std::string>("next tag", flow_next_tag.get());
}

// porosity used with velocity to compute particle velocity when dispersion is on
Teuchos::ParameterList& poro_list_next = S_->GetEvaluatorList(Keys::getKey("porosity", transport_next_tag));
poro_list_next.set<std::string>("evaluator type", "alias");
poro_list_next.set<std::string>("target", Keys::getKey("porosity", flow_next_tag, true));
if (!poro_list_next.isParameter("evaluator type")) {
poro_list_next.set<std::string>("evaluator type", "temporal interpolation");
poro_list_next.set<std::string>("current tag", flow_current_tag.get());
poro_list_next.set<std::string>("next tag", flow_next_tag.get());
}

// chemistry uses (independently) density, saturation, and porosity
if (chemistry) {
Teuchos::ParameterList& dens_list_current = S_->GetEvaluatorList(Keys::getKey("mass_density_liquid", transport_current_tag));
if (!dens_list_current.isParameter("evaluator type")) {
dens_list_current.set<std::string>("evaluator type", "temporal interpolation");
dens_list_current.set<std::string>("current tag", flow_current_tag.get());
dens_list_current.set<std::string>("next tag", flow_next_tag.get());
}

Teuchos::ParameterList& sat_list_current = S_->GetEvaluatorList(Keys::getKey("saturation_liquid", transport_current_tag));
if (!sat_list_current.isParameter("evaluator type")) {
sat_list_current.set<std::string>("evaluator type", "temporal interpolation");
sat_list_current.set<std::string>("current tag", flow_current_tag.get());
sat_list_current.set<std::string>("next tag", flow_next_tag.get());
}

Teuchos::ParameterList& poro_list_current = S_->GetEvaluatorList(Keys::getKey("porosity", transport_current_tag));
if (!poro_list_current.isParameter("evaluator type")) {
poro_list_current.set<std::string>("evaluator type", "temporal interpolation");
poro_list_current.set<std::string>("current tag", flow_current_tag.get());
poro_list_current.set<std::string>("next tag", flow_next_tag.get());
}
}
}
}

Expand All @@ -95,43 +131,105 @@ MPCFlowTransport::parseParameterList()
// hard-code this as the default. If this breaks in the future it can be
// fixed. --ETC
Teuchos::ParameterList& flux_list = S_->GetEvaluatorList(Keys::getKey("surface-water_flux", transport_next_tag));
flux_list.set<std::string>("evaluator type", "alias");
flux_list.set<std::string>("target", Keys::getKey("surface-water_flux", flow_next_tag, true));
if (!flux_list.isParameter("evaluator type")) {
flux_list.set<std::string>("evaluator type", "alias");
flux_list.set<std::string>("target", Keys::getKey("surface-water_flux", flow_next_tag, true));
}

// velocity for dispersivity
Teuchos::ParameterList& velo_list = S_->GetEvaluatorList(Keys::getKey("surface-water_velocity", transport_next_tag));
velo_list.set<std::string>("evaluator type", "alias");
velo_list.set<std::string>("target", Keys::getKey("surface-water_velocity", flow_next_tag, true));
if (!velo_list.isParameter("evaluator type")) {
velo_list.set<std::string>("evaluator type", "alias");
velo_list.set<std::string>("target", Keys::getKey("surface-water_velocity", flow_next_tag, true));
}

// set the liquid water content as an interpolated field
// note that current copy is kept by transport PK
Teuchos::ParameterList& lwc_list_next = S_->GetEvaluatorList(Keys::getKey(surf_lwc_key, transport_next_tag));
lwc_list_next.set<std::string>("evaluator type", "temporal interpolation");
lwc_list_next.set<std::string>("current tag", flow_current_tag.get());
lwc_list_next.set<std::string>("next tag", flow_next_tag.get());
if (!lwc_list_next.isParameter("evaluator type")) {
lwc_list_next.set<std::string>("evaluator type", "temporal interpolation");
lwc_list_next.set<std::string>("current tag", flow_current_tag.get());
lwc_list_next.set<std::string>("next tag", flow_next_tag.get());
}

// porosity used with velocity to compute particle velocity when dispersion is on
Teuchos::ParameterList& poro_list_next = S_->GetEvaluatorList(Keys::getKey("surface-porosity", transport_next_tag));
if (!poro_list_next.isParameter("evaluator type")) {
poro_list_next.set<std::string>("evaluator type", "temporal interpolation");
poro_list_next.set<std::string>("current tag", flow_current_tag.get());
poro_list_next.set<std::string>("next tag", flow_next_tag.get());
}

// chemistry uses (independently) density, saturation, and porosity at the current tag
if (chemistry) {
Teuchos::ParameterList& dens_list_current = S_->GetEvaluatorList(Keys::getKey("surface-mass_density_liquid", transport_current_tag));
if (!dens_list_current.isParameter("evaluator type")) {
dens_list_current.set<std::string>("evaluator type", "temporal interpolation");
dens_list_current.set<std::string>("current tag", flow_current_tag.get());
dens_list_current.set<std::string>("next tag", flow_next_tag.get());
}

Teuchos::ParameterList& sat_list_current = S_->GetEvaluatorList(Keys::getKey("surface-ponded_depth", transport_current_tag));
if (!sat_list_current.isParameter("evaluator type")) {
sat_list_current.set<std::string>("evaluator type", "temporal interpolation");
sat_list_current.set<std::string>("current tag", flow_current_tag.get());
sat_list_current.set<std::string>("next tag", flow_next_tag.get());
}

Teuchos::ParameterList& poro_list_current = S_->GetEvaluatorList(Keys::getKey("surface-porosity", transport_current_tag));
if (!poro_list_current.isParameter("evaluator type")) {
poro_list_current.set<std::string>("evaluator type", "temporal interpolation");
poro_list_current.set<std::string>("current tag", flow_current_tag.get());
poro_list_current.set<std::string>("next tag", flow_next_tag.get());
}
}

}
}

MPCSubcycled::parseParameterList();

if (subsurface) {
// first require lwc@NEXT -- otherwise it will get called to be an alias
// lwc@TRANSPORT_NEXT, which is an interpolator that depends upon this.
// for interpolated evaluators, we must first require key@NEXT -- otherwise
// it will get called to be an alias key@TRANSPORT_NEXT, which is an
// interpolator that depends upon this.
requireAtNext(sub_lwc_key, Tags::NEXT, *S_);
requireAtNext("porosity", Tags::NEXT, *S_);
if (chemistry) {
requireAtNext("mass_density_liquid", Tags::NEXT, *S_);
requireAtNext("saturation_liquid", Tags::NEXT, *S_);
}

// now call at lwc@transport_next, which will be the interpolant
// now require key@transport_next, which will be the interpolant
requireAtNext(sub_lwc_key, transport_next_tag, *S_);
requireAtNext("porosity", transport_next_tag, *S_);
if (chemistry) {
requireAtCurrent("mass_density_liquid", transport_current_tag, *S_);
requireAtCurrent("saturation_liquid", transport_current_tag, *S_);
requireAtCurrent("porosity", transport_current_tag, *S_);
}
}

if (surface) {
// first require lwc@NEXT -- otherwise it will get called to be an alias
// lwc@TRANSPORT_NEXT, which is an interpolator that depends upon this.
// for interpolated evaluators, we must first require key@NEXT -- otherwise
// it will get called to be an alias key@TRANSPORT_NEXT, which is an
// interpolator that depends upon this.
requireAtNext(surf_lwc_key, Tags::NEXT, *S_);
requireAtNext("surface-porosity", Tags::NEXT, *S_);
if (chemistry) {
requireAtNext("surface-mass_density_liquid", Tags::NEXT, *S_);
requireAtNext("surface-ponded_depth", Tags::NEXT, *S_);
}

// now call at lwc@transport_next, which will be the interpolant
// now require key@transport_next, which will be the interpolant
requireAtNext(surf_lwc_key, transport_next_tag, *S_);
requireAtNext("surface-porosity", transport_next_tag, *S_);
if (chemistry) {
requireAtCurrent("surface-mass_density_liquid", transport_current_tag, *S_);
requireAtCurrent("surface-ponded_depth", transport_current_tag, *S_);
requireAtCurrent("surface-porosity", transport_current_tag, *S_);
}
}

}


Expand Down
12 changes: 11 additions & 1 deletion src/pks/mpc/mpc_reactivetransport.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,13 +40,23 @@ MPCReactiveTransport::parseParameterList()
Teuchos::Array<std::string> names = plist_->get<Teuchos::Array<std::string>>("PKs order");
domain_ = getSubPKPlist_(1)->get<std::string>("domain name", "domain");

// chemistry and transport share the same primary variable
// chemistry and transport share the same primary variable, we take it from transport
tcc_key_ = Keys::readKey(*getSubPKPlist_(1), domain_, "primary variable key",
"total_component_concentration");
getSubPKPlist_(0)->set<std::string>("primary variable key", tcc_key_);

// both pks need access to the primary variable
getSubPKPlist_(0)->set<std::string>("primary variable password", name_);
getSubPKPlist_(1)->set<std::string>("primary variable password", name_);

// tell chemistry to operator split
getSubPKPlist_(0)->set("operator split", true);

// Only one PK needs to set the IC -- we use chemistry to process geochemical
// initial conditions -- but the "initial conditions" list must be present
// for all PK_Physical PKs, so just touch it here to make sure it exists.
getSubPKPlist_(1)->sublist("initial conditions");

mol_dens_key_ = Keys::readKey(*plist_, domain_, "molar density liquid", "molar_density_liquid");

cast_sub_pks_();
Expand Down
2 changes: 1 addition & 1 deletion testing/ats-regression-tests
Submodule ats-regression-tests updated 184 files
1 change: 1 addition & 0 deletions tools/utils/plot_lines.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ def plotLines(x, y, t, ax=None, colorbar=True, colorbar_ticks=True, colorbar_lab
if colorbar_ticks:
xticks = axcb.get_ticks()
new_ticks = [t * (t_max - t_min) + t_min for t in xticks]
axcb.set_ticks(new_ticks)
axcb.set_ticklabels([str(np.round(t, 2)) for t in new_ticks])
else:
axcb.set_ticks(list())
Expand Down

0 comments on commit f37124b

Please sign in to comment.