-
Notifications
You must be signed in to change notification settings - Fork 248
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.
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
[GeoMechanicsApplication] Clean up / simplify the calculation of the permeability update factor #12372
Conversation
Use the member functions' return type rather than assigning it to a member of an output argument.
No longer pass a reference to an `ElementVariables` instance, but only a const-reference to a strain vector.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Small and clear PR which once more reduces our use of ElementVariables. Also, good catch with the std::pow
!
rVariables.PermeabilityUpdateFactor = pow(10.0, permLog10); | ||
} else { | ||
rVariables.PermeabilityUpdateFactor = 1.0; | ||
return std::pow(10.0, permLog10); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Good catch, it seems that before it found a boost version (in corecrt_math.h) and now just uses cmath. It's good to keep this in mind indeed.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi Anne,
This touches on the same lines as Richards pull request for strain computation.
Please be careful when merging.
Regards, Wijtze Pieter
rVariables.PermeabilityUpdateFactor = 1.0; | ||
if (const auto& r_prop = this->GetProperties(); r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) { | ||
const double InverseCK = r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR]; | ||
const double epsV = StressStrainUtilities::CalculateTrace(rStrainVector); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This resolves a lot of created confusion about the meaning of the volumetric strain. There was no need to mingle it with names indicating stress.
double SmallStrainUPwDiffOrderElement::CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const | ||
{ | ||
KRATOS_TRY | ||
|
||
const PropertiesType& rProp = this->GetProperties(); | ||
|
||
if (rProp[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) { | ||
const double InverseCK = rProp[PERMEABILITY_CHANGE_INVERSE_FACTOR]; | ||
const double epsV = StressStrainUtilities::CalculateTrace(rVariables.StrainVector); | ||
const double ePrevious = rProp[POROSITY] / (1.0 - rProp[POROSITY]); | ||
if (const auto& r_prop = this->GetProperties(); r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) { | ||
const double InverseCK = r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR]; | ||
const double epsV = StressStrainUtilities::CalculateTrace(rStrainVector); | ||
const double ePrevious = r_prop[POROSITY] / (1.0 - r_prop[POROSITY]); | ||
const double eCurrent = (1.0 + ePrevious) * std::exp(epsV) - 1.0; | ||
const double permLog10 = (eCurrent - ePrevious) * InverseCK; | ||
rVariables.PermeabilityUpdateFactor = pow(10.0, permLog10); | ||
} else { | ||
rVariables.PermeabilityUpdateFactor = 1.0; | ||
return std::pow(10.0, permLog10); | ||
} | ||
|
||
return 1.0; | ||
|
||
KRATOS_CATCH("") | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Repeated in several elements. I do not see element specific code, should this be a utility function. The strain dependency points towards equation of motion or stress_strain utilities. Its use points to the transport equation utilities.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I agree that we can move this to a separate utility function. I would suggest to do that in a separate PR, just to keep the history clear.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Anne, nice work in getting rid of Variables and code smells. After the changes made, it looks that both elements use almost the same function. Is it useful to put it in a common place? Thank you.
if (const auto& r_prop = this->GetProperties(); r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) { | ||
const double InverseCK = r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR]; | ||
const double epsV = StressStrainUtilities::CalculateTrace(rStrainVector); | ||
const double ePrevious = r_prop[POROSITY] / (1.0 - r_prop[POROSITY]); | ||
const double eCurrent = (1.0 + ePrevious) * std::exp(epsV) - 1.0; | ||
const double permLog10 = (eCurrent - ePrevious) * InverseCK; | ||
return std::pow(10.0, permLog10); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It looks that both elements use the same function and a difference is only in this->GetProperties. Is it useful to move the function to utilities with r_Prop as input?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, I agree that we can move this to a separate utility function. I would suggest to do that in a separate PR, just to keep the history clear. As you indicated, we'd better provide the element properties as an input parameter to that new utility function.
📝 Description
Simplified and cleaned up the implementation of the members that compute the permeability update factor.
🆕 Changelog
const
.else
keywords.std::pow
rather than C-stylepow
.