Water/steam model for PorousFlow #13163
Conversation
|
Job Documentation on 5329425 wanted to post the following: View the site here This comment will be updated on new commits. |
|
Hi @lindsayad I ran into an issue when I had something like I verified that doing something like works as expected (outputs 1). The stack trace looks like: |
|
What's the easiest way for me to reproduce? |
|
Hi @WilkAndy, I added some write up of some tests that I have added - see https://mooseframework.inl.gov/docs/PRs/13163/site/modules/porous_flow/verification_problems/index.html Is this the sort of documentation that you were thinking of for these tests that have analytical or benchmark solutions? |
|
Yea, this looks excellent (great agreement with external sources too!). My current setup is: for all tests in porous_flow/test/tests/some_dir, i write a markdown document porous_flow/doc/content/modules/porous_flow/tests/some_dir/some_dir_tests.md that describes all the tests. In that "doc" directory i also keep all the png files and python files needed to create those png files for the some_dir tests. By the way, is it possible to compare the CPU time in MOOSE with TOUGH for these tests, or is that annoying/impossible? |
|
Shouldn't be too hard to do |
|
@lindsayad - I haven't been able to come up with a simple example yet. I'll try again tomorrow |
|
Hi @lindsayad I can't replicate this fpe, so I guess that it is happening somewhere else. In case you wanted to take a look, you can grab this branch and apply the following patch: Then compile in dbg. The test |
cpgr
force-pushed
the
pf_watervapor
branch
3 times, most recently
from
c932036
to
d9ed1c4
Compare
| design = 'PorousFlowFluidPropertyIC.md' | ||
| issues = '#13108' | ||
| [../] | ||
| [./fluidpropic_celcius] |
| [../] | ||
| [./fluidpropic_celcius] | ||
| type = 'CSVDiff' | ||
| input = 'fluidpropic_celcius.i' |
| # Test the correct calculation of fluid properties using PorousFlwoFluidPropertyIC | ||
| # | ||
| # Variables: | ||
| # Pressure: 1 Mpa |
There was a problem hiding this comment.
Mpa -> MPa, also in the other file in the patch
| @@ -0,0 +1,114 @@ | |||
| # Test the correct calculation of fluid properties using PorousFlwoFluidPropertyIC | |||
There was a problem hiding this comment.
PorousFlwoFluidPropertyIC -> PorousFlowFluidPropertyIC
| params.addRequiredParam<UserObjectName>("fp", "The name of the user object for the fluid"); | ||
| MooseEnum property_enum("enthalpy internal_energy density"); | ||
| params.addRequiredParam<MooseEnum>( | ||
| "property", property_enum, "The fluid property that this auxillary kernel is to calculate"); |
There was a problem hiding this comment.
You mean "initial condition" not "auxiliary kernel", right?
|
|
||
| /** | ||
| * Calculates all required thermophysical properties and derivatives for each phase | ||
| * and fluid component. Must override in all derived classes. |
There was a problem hiding this comment.
You say this must be overridden, but you provide an implementation. That can confuse your users...
There was a problem hiding this comment.
Good catch. This was a copy and paste as well, but the other one should be fixed as well. This method doesn't need to be overridden in most cases
| @@ -277,7 +277,7 @@ class SinglePhaseFluidProperties : public FluidProperties | |||
| * @return internal energy (J/kg) | |||
| */ | |||
| virtual Real e_from_p_T(Real p, Real T) const; | |||
| virtual DualReal e_from_p_T(DualReal p, DualReal T) const; | |||
| DualReal e_from_p_T(const DualReal & p, const DualReal & T) const; | |||
There was a problem hiding this comment.
Are these not virtual, because the way you compute it is always the same and you actually rely on calling the virtual versions inside the methods? If so, we need this to be mentioned in the doco for this class.
There was a problem hiding this comment.
I did this to be consistent with the definitions in SinglePhaseFluidProperties.h in the propfunc macro for the other methods.
I'm not sure that this is the most flexible design though. For the simple fluids like the ideal gas, it is trivial to do the derivatives, so this idea makes sense. For more complicated fluids (like the ones I use), the derivatives are a pain to code, so it might be nice to code the DualReal version instead.
There was a problem hiding this comment.
Let me demonstrate the idea behind propfunc:
- let's say we need to compute
e_from_p_Tand our primary variables arerho,rhou,rhoE. - The
DualRealforphas this structure:(value of p, dp_drho, dp_drhou, dp_rhoE), similarly forT. - The value we need back for
eis:(value of e, de_drho, de_drhou, de_drhoE). - So we need to compute:
(e_from_p_T, de_dp * dp_drho + de_dT * dT_drho, de_dp * dp_drhou + de_dT * dT_drhou, de_dp * dp_drhoE + de_dT * dT_drhoE).
The d{p|T}_d{rho|rhou|rhoE} is coming from the AD system and de_dp and de_dT must be computed.
If de_dp and de_dT are difficult to compute (as you say), I am not sure that computing the DualReal version would be easier.
Also, if we are solving in terms of p and T, then the dp_dp and dT_dT from the AD system should be 1 and the chain rule still holds.
Keep in mind that this API has to be independent from the primary variables used. It looks a little bit that you are trying to invert the design and use the AD system to compute the non-AD value, like in T_from_p_h. That would create an inconsistency and this API would be usable only with p, h as primary variables. We do not want that.
May be, we need to look at what you are trying to do and how to fit it into this class.
There was a problem hiding this comment.
What I added here is basically a copy and paste from the definition in propfunc but with p1 and p2 renamed p and T. Are you just saying that it should be declared as
DualReal e_from_p_T(const DualReal & p1, const DualReal & p2) const
etc?
There was a problem hiding this comment.
I have just had a look at what I did here and I think that it implements exactly what you describe, and what is implemented in propfunc for all of the other methods. The input DualReals p and T can be functions of whatever with all of the derivatives carried through, and they are multiplied by the derivatives de_dp and de_dT computed in the e_from_p_T(p, T, e, de_dp, de_dT) method.
e_from_p_T and T_from_p_h aren't currently declared using propfunc as they had default implementations, so currently there is no AD version of these methods (I guess that they aren't declared using propfunc(e,p,T) as you can't declare the methods twice?)
I do compute the (Real) derivatives in WaterFluidProperties::T_from_p_h(p, h, T, dT_dp, dT_dh) using a seperate AD method to avoid coding all the derivatives for all of the subregions, but the public AD method T_from_p_h(p, h) is unchanged, so should work with p and h defined in terms of other variables as you described above.
There was a problem hiding this comment.
I think that it implements exactly what you describe, and what is implemented in
propfuncfor all of the other methods.
So, why did not you use the propfunc macro? When we add the derivative version, then the macro should be used.
There was a problem hiding this comment.
Because these methods are implemented in SinglePhaseFluidProperties (with a default implementation), they can't be declared using propfunc because they would be defined twice in the same class.
I could add a SinglePhaseFluidPropertiesBase that has all of the methods declared using propfunc and then override the ones with default implementations in SinglePhaseFluidProperties to get around this I guess.
There was a problem hiding this comment.
I finally understand why you did it this way. I forgot that the default implementation of e_from_p_T is actually a useful one. It does not throw the error of not being implemented...
Sorry it took so long...
|
@cpgr - David is on vacation for the remainder of this week. If you need a review on this sooner. Please let me know. |
|
No rush @permcody - it is my fault as I was away all of last week. |
Did you also do this verification in // if_else is necessary here to handle cases where a is negative but b
// is 0; we should have a contribution of 0 from those, not NaN.
DualNumber_std_binary(pow,
funcval * (b.value() * a.derivatives() / a.value() +
MetaPhysicL::if_else(b.derivatives(), b.derivatives() * std::log(a.value()), b.derivatives())))So the |
|
Wow, isn't that a sneaky one. I hadn't thought to test MetaPhysicL with FPEs enabled, since I usually just use tests for NaN myself. And even if I had turned FPEs on, I'm not certain we have internal test coverage for The reason for MetaPhysicL::if_else in general isn't to avoid manual Without expression-template-based types (which would be too hard to make happen on every class we want to use) there's no way to avoid making that log call there. And all the other options I can see are lousy. We could have a compile-time option to temporarily disable FPEs for that call, but that would mean multiple system calls trashing performance for scalar types (although if you're only enabling FPEs in debug modes maybe that's acceptable??). We could truncate the range of a before calling the log to get effectively a log_without_FPEs, then examine a in an outer if_else call and put quiet_NaNs back in place where needed, but that would mean having to fall back on NaN infectiousness to avoid falsely trusting your output, losing that FPE signal and its help debugging... |
In preparation for a water/vapor EOS, enable AD calculation of some other fluid properties Refs idaholab#13108
A single component fluid state does not need to have a compositional flash, so rearrange the base classes to make the inheritance more suitable for both single and multicomponent fluid states. Includes a new multicomponent base class that new multicomponent fluid states should derive from. Refs idaholab#13108
Takes the formulation in a single component fluid state, like PorousFlowWaterVapor, and translate the fluid properties into material properties for use in the existing kernels. Refs idaholab#13108
Given a pressure and temperature, set an IC for a fluid property (such as enthalpy) using a FluidProperties UserObject. Refs idaholab#13108
Simplified tests of verification problems Refs idaholab#13108
An fpe was thrown in dbg when std::pow(-0.x, 0) was encountered. Factoring out the negative seems to work around this. Refs idaholab#13108
Calculate the AD version in the SinglePhaseFluidProperties base class as with the other properties Refs idaholab#13108
Implements a water/steam model based on pressure and enthalpy using AD to compute all of the tricky derivatives.
A lot of this is some machinery to make it easier to implement this model - the actual water/steam model is only a few hundred lines - as well as documentation and tests
This model will be useful for geothermal, as it allows phase changes to occur (from water to vapour and vice-versa).
Closes #13108