Skip to content

Commit

Permalink
Merge pull request #896 from akva2/floats
Browse files Browse the repository at this point in the history 
Add necessary changes to build float simulators
  • Loading branch information
@bska
bska authored Aug 30, 2024
2 parents 629d790 + 3e96dea commit b7e9476
Show file tree
Hide file tree
Showing 6 changed files with 50 additions and 41 deletions.
31 changes: 17 additions & 14 deletions opm/models/blackoil/blackoilconvectivemixingmodule.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,15 +28,16 @@
#ifndef EWOMS_CONVECTIVEMIXING_MODULE_HH
#define EWOMS_CONVECTIVEMIXING_MODULE_HH

#include <opm/models/discretization/common/fvbaseproperties.hh>
#include "opm/material/common/MathToolbox.hpp"
#include <dune/common/fvector.hh>

#include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp>
#include <opm/input/eclipse/Schedule/Schedule.hpp>
#include <opm/material/common/Valgrind.hpp>

#include <dune/common/fvector.hh>
#include <opm/material/common/Valgrind.hpp>

#include <stdexcept>
#include <iostream>
#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>

namespace Opm {

Expand Down Expand Up @@ -119,8 +120,8 @@ public:
const Scalar,
const ConvectiveMixingModuleParam&)
{}

};

template <class TypeTag>
class BlackOilConvectiveMixingModule<TypeTag, /*enableConvectiveMixing=*/true>
{
Expand Down Expand Up @@ -150,7 +151,10 @@ public:
};

#if HAVE_ECL_INPUT
static void beginEpisode(const EclipseState& eclState, const Schedule& schedule, const int episodeIdx, ConvectiveMixingModuleParam& info)
static void beginEpisode(const EclipseState& eclState,
const Schedule& schedule,
const int episodeIdx,
ConvectiveMixingModuleParam& info)
{
// check that Xhi and Psi didn't change
std::size_t numRegions = eclState.runspec().tabdims().getNumPVTTables();
Expand Down Expand Up @@ -211,8 +215,10 @@ public:
const auto salt_ex = Toolbox::value(intQuantsEx.fluidState().saltConcentration());

const auto bLiquidEx = (FluidSystem::phaseIsActive(waterPhaseIdx)) ?
FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(), t_ex, p_ex, 0.0, salt_ex) :
FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(), t_ex, p_ex, 0.0);
FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
t_ex, p_ex, Scalar{0.0}, salt_ex) :
FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsEx.pvtRegionIndex(),
t_ex, p_ex, Scalar{0.0});

const auto& refDensityLiquidEx = (FluidSystem::phaseIsActive(waterPhaseIdx)) ?
FluidSystem::waterPvt().waterReferenceDensity(intQuantsEx.pvtRegionIndex()) :
Expand Down Expand Up @@ -373,12 +379,9 @@ public:
}
}
}
};

}
};

}
#endif



#endif
8 changes: 7 additions & 1 deletion opm/models/blackoil/blackoilmicpmodules.hh
View file
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
Original file line number Diff line number Diff line change
Expand Up @@ -324,13 +324,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 @@ -341,8 +341,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 @@ -352,35 +352,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 @@ -398,15 +398,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
Original file line number Diff line number Diff line change
Expand Up @@ -849,7 +849,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 @@ -876,7 +876,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 @@ -925,10 +925,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 @@ -981,7 +981,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
Original file line number Diff line number Diff line change
Expand Up @@ -1633,7 +1633,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
Original file line number Diff line number Diff line change
Expand Up @@ -454,13 +454,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

0 comments on commit b7e9476

Please sign in to comment.