Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add necessary changes to build float simulators #896

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

Add necessary changes to build float simulators #896

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
85689a1
BlackoilMICPModules: fix field prop copy if Scalar is float
akva2 Apr 11, 2024
24e2c47
consistently use Scalar type
akva2 Apr 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
8 changes: 7 additions & 1 deletion opm/models/blackoil/blackoilmicpmodules.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,13 @@ public:
MICPpara.getMaximumUreaConcentration(),
MICPpara.getToleranceBeforeClogging());
// obtain the porosity for the clamp in the blackoilnewtonmethod
params_.phi_ = eclState.fieldProps().get_double("PORO");
if constexpr (std::is_same_v<Scalar, float>) {
const auto phi = eclState.fieldProps().get_double("PORO");
params_.phi_.resize(phi.size());
std::copy(phi.begin(), phi.end(), params_.phi_.begin());
} else {
params_.phi_ = eclState.fieldProps().get_double("PORO");
}
}
#endif

Expand Down
32 changes: 16 additions & 16 deletions opm/models/blackoil/blackoilnewtonmethod.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -407,13 +407,13 @@ protected:
else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
// z fraction updates are also subject to the Appleyard chop
const auto& curr = currentValue[Indices::zFractionIdx]; // or currentValue[pvIdx] given the block condition
delta = std::clamp(delta, curr - 1.0, curr);
delta = std::clamp(delta, curr - Scalar{1.0}, curr);
}
else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
const double sign = delta >= 0. ? 1. : -1.;
// maximum change of polymer molecular weight, the unit is MDa.
// applying this limit to stabilize the simulation. The value itself is still experimental.
const double maxMolarWeightChange = 100.0;
const Scalar maxMolarWeightChange = 100.0;
delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
delta *= satAlpha;
}
Expand All @@ -424,8 +424,8 @@ protected:
else if (enableBrine && pvIdx == Indices::saltConcentrationIdx &&
enableSaltPrecipitation &&
currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
const double maxSaltSaturationChange = 0.1;
const double sign = delta >= 0. ? 1. : -1.;
const Scalar maxSaltSaturationChange = 0.1;
const Scalar sign = delta >= 0. ? 1. : -1.;
delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
}

Expand All @@ -435,35 +435,35 @@ protected:
// keep the solvent saturation between 0 and 1
if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss )
nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
}

// keep the z fraction between 0 and 1
if (enableExtbo && pvIdx == Indices::zFractionIdx)
nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});

// keep the polymer concentration above 0
if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx)
nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});

if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
if (polymerConcentration < 1.e-10)
nextValue[pvIdx] = 0.0;
}

// keep the foam concentration above 0
if (enableFoam && pvIdx == Indices::foamConcentrationIdx)
nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});

if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
// keep the salt concentration above 0
if (!enableSaltPrecipitation || (enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs))
nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
// keep the salt saturation below upperlimit
if ((enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp))
nextValue[pvIdx] = std::min(nextValue[pvIdx], 1.0-1.e-8);
nextValue[pvIdx] = std::min(nextValue[pvIdx], Scalar{1.0-1.e-8});
}

// keep the temperature within given values
Expand All @@ -481,15 +481,15 @@ protected:
// concentration (the urea concentration has been scaled by 10). For
// the biofilm and calcite, we set this value equal to the porosity minus the clogging tolerance.
if (enableMICP && pvIdx == Indices::microbialConcentrationIdx)
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::densityBiofilm());
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::densityBiofilm());
if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx)
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::maximumOxygenConcentration());
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumOxygenConcentration());
if (enableMICP && pvIdx == Indices::ureaConcentrationIdx)
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::maximumUreaConcentration());
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumUreaConcentration());
if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx)
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
if (enableMICP && pvIdx == Indices::calciteConcentrationIdx)
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
}

// switch the new primary variables to something which is physically meaningful.
Expand Down
14 changes: 7 additions & 7 deletions opm/models/blackoil/blackoilprimaryvariables.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -857,7 +857,7 @@ public:
soMax);

Scalar rs = (*this)[Indices::compositionSwitchIdx];
if (rs > std::min(rsMax, rsSat*(1.0 + eps))) {
if (rs > std::min(rsMax, rsSat*(Scalar{1.0} + eps))) {
// the gas phase appears, i.e., switch the primary variables to GasMeaning::Sg
setPrimaryVarsMeaningGas(GasMeaning::Sg);
(*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
Expand All @@ -884,7 +884,7 @@ public:
soMax);

Scalar rv = (*this)[Indices::compositionSwitchIdx];
if (rv > std::min(rvMax, rvSat*(1.0 + eps))) {
if (rv > std::min(rvMax, rvSat*(Scalar{1.0} + eps))) {
// switch to phase equilibrium mode because the oil phase appears. here
// we also need the capillary pressures to calculate the oil phase
// pressure using the gas phase pressure
Expand Down Expand Up @@ -933,10 +933,10 @@ public:
ssol =(*this) [Indices::solventSaturationIdx];

Scalar so = 1.0 - sw - sg - ssol;
sw = std::min(std::max(sw,0.0),1.0);
so = std::min(std::max(so,0.0),1.0);
sg = std::min(std::max(sg,0.0),1.0);
ssol = std::min(std::max(ssol,0.0),1.0);
sw = std::min(std::max(sw, Scalar{0.0}), Scalar{1.0});
so = std::min(std::max(so, Scalar{0.0}), Scalar{1.0});
sg = std::min(std::max(sg, Scalar{0.0}), Scalar{1.0});
ssol = std::min(std::max(ssol, Scalar{0.0}), Scalar{1.0});
Scalar st = sw + so + sg + ssol;
sw = sw/st;
sg = sg/st;
Expand Down Expand Up @@ -989,7 +989,7 @@ public:
template<class Serializer>
void serializeOp(Serializer& serializer)
{
using FV = Dune::FieldVector<double,getPropValue<TypeTag, Properties::NumEq>()>;
using FV = Dune::FieldVector<Scalar, getPropValue<TypeTag, Properties::NumEq>()>;
serializer(static_cast<FV&>(*this));
serializer(primaryVarsMeaningWater_);
serializer(primaryVarsMeaningPressure_);
Expand Down
2 changes: 1 addition & 1 deletion opm/models/discretization/common/fvbasediscretization.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1690,7 +1690,7 @@ public:
}

// do the normalization of the relative error
Scalar alpha = std::max(1e-20,
Scalar alpha = std::max(Scalar{1e-20},
std::max(std::abs(maxRelErr),
std::abs(minRelErr)));
for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx)
Expand Down
4 changes: 2 additions & 2 deletions opm/models/nonlinear/newtonmethod.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -475,13 +475,13 @@ public:
if (numIterations_ > targetIterations_()) {
Scalar percent = Scalar(numIterations_ - targetIterations_())/targetIterations_();
Scalar nextDt = std::max(problem().minTimeStepSize(),
oldDt/(1.0 + percent));
oldDt / (Scalar{1.0} + percent));
return nextDt;
}

Scalar percent = Scalar(targetIterations_() - numIterations_)/targetIterations_();
Scalar nextDt = std::max(problem().minTimeStepSize(),
oldDt*(1.0 + percent/1.2));
oldDt*(Scalar{1.0} + percent / Scalar{1.2}));
return nextDt;
}

Expand Down