Added option for projection saturation. It do not change the result for default setting. To be able to do this it

is slightly unlogical for surfactant runs which do projection two times.
This commit is contained in:
hnil 2020-05-07 20:44:56 +02:00
parent 4ae340b953
commit 09d58b770a
2 changed files with 85 additions and 0 deletions

View File

@ -39,10 +39,12 @@ BEGIN_PROPERTIES
NEW_PROP_TAG(DpMaxRel);
NEW_PROP_TAG(DsMax);
NEW_PROP_TAG(PriVarOscilationThreshold);
NEW_PROP_TAG(ProjectSaturations);
SET_SCALAR_PROP(NewtonMethod, DpMaxRel, 0.3);
SET_SCALAR_PROP(NewtonMethod, DsMax, 0.2);
SET_SCALAR_PROP(NewtonMethod, PriVarOscilationThreshold, 1e-5);
SET_BOOL_PROP(NewtonMethod, ProjectSaturations, false);
END_PROPERTIES
@ -74,6 +76,7 @@ public:
priVarOscilationThreshold_ = EWOMS_GET_PARAM(TypeTag, Scalar, PriVarOscilationThreshold);
dpMaxRel_ = EWOMS_GET_PARAM(TypeTag, Scalar, DpMaxRel);
dsMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, DsMax);
projectSaturations_ = EWOMS_GET_PARAM(TypeTag, bool, ProjectSaturations);
}
/*!
@ -98,6 +101,7 @@ public:
EWOMS_REGISTER_PARAM(TypeTag, Scalar, DsMax, "Maximum absolute change of any saturation in a single iteration");
EWOMS_REGISTER_PARAM(TypeTag, Scalar, PriVarOscilationThreshold,
"The threshold value for the primary variable switching conditions after its meaning has switched to hinder oscilations");
EWOMS_REGISTER_PARAM(TypeTag,bool, ProjectSaturations, "Option for doing saturation projection");
}
/*!
@ -309,6 +313,9 @@ protected:
if (wasSwitched_[globalDofIdx])
++ numPriVarsSwitched_;
if(projectSaturations_){
nextValue.chopAndNormalizeSaturations();
}
nextValue.checkDefined();
}
@ -319,6 +326,7 @@ private:
Scalar priVarOscilationThreshold_;
Scalar dpMaxRel_;
Scalar dsMax_;
bool projectSaturations_;
// keep track of cells where the primary variable meaning has changed
// to detect and hinder oscillations

View File

@ -565,6 +565,83 @@ public:
return false;
}
bool chopAndNormalizeSaturations(){
if (primaryVarsMeaning() == OnePhase_p){
return false;
}
Scalar Sw = 0.0;
if (waterEnabled)
Sw = (*this)[Indices::waterSaturationIdx];
if (primaryVarsMeaning() == Sw_po_Sg) {
Scalar Sg = 0.0;
if (compositionSwitchEnabled)
Sg = (*this)[Indices::compositionSwitchIdx];
Scalar Ssol = 0.0;
if (enableSolvent)
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);
Scalar St = Sw + So + Sg + Ssol;
Sw = Sw/St;
Sg = Sg/St;
Ssol = Ssol/St;
assert(St>0.5);
if (waterEnabled)
(*this)[Indices::waterSaturationIdx]= Sw;
if (compositionSwitchEnabled)
(*this)[Indices::compositionSwitchIdx] = Sg;
if (enableSolvent)
(*this)[Indices::solventSaturationIdx] = Ssol;
return not(St==1);
}else if (primaryVarsMeaning() == Sw_po_Rs) {
Scalar Ssol = 0.0;
if (enableSolvent)
Ssol = (*this)[Indices::solventSaturationIdx];
Scalar So = 1.0 - Sw - Ssol;
Sw = std::min(std::max(Sw,0.0),1.0);
So = std::min(std::max(So,0.0),1.0);
Ssol = std::min(std::max(Ssol,0.0),1.0);
//Sg = 0.0;
Scalar St = Sw + So + Ssol;
assert(St>0.5);
Sw=Sw/St;
Ssol=Ssol/St;
if (waterEnabled)
(*this)[Indices::waterSaturationIdx]= Sw;
if (enableSolvent)
(*this)[Indices::solventSaturationIdx] = Ssol;
return not(St==1);
}else{
assert(primaryVarsMeaning() == Sw_pg_Rv);
assert(compositionSwitchEnabled);
Scalar Ssol=0.0;
if (enableSolvent)
Ssol = (*this)[Indices::solventSaturationIdx];
Scalar Sg = 1.0 - Sw - Ssol;
//So = 0.0;
Sw = std::min(std::max(Sw,0.0),1.0);
Sg = std::min(std::max(Sg,0.0),1.0);
Ssol = std::min(std::max(Ssol,0.0),1.0);
Scalar St = Sw + Sg + Ssol;
assert(St>0.5);
Sw=Sw/St;
Ssol=Ssol;
if (waterEnabled)
(*this)[Indices::waterSaturationIdx]= Sw;
if (enableSolvent)
(*this)[Indices::solventSaturationIdx] = Ssol;
return not(St==1);
}
assert(false);
}
BlackOilPrimaryVariables& operator=(const BlackOilPrimaryVariables& other) = default;
BlackOilPrimaryVariables& operator=(Scalar value)
{