Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PorousFlowFluidMass.h"
11 :
12 : #include "MooseVariable.h"
13 :
14 : #include "libmesh/quadrature.h"
15 :
16 : registerMooseObject("PorousFlowApp", PorousFlowFluidMass);
17 : registerMooseObject("PorousFlowApp", ADPorousFlowFluidMass);
18 :
19 : template <bool is_ad>
20 : InputParameters
21 5235 : PorousFlowFluidMassTempl<is_ad>::validParams()
22 : {
23 5235 : InputParameters params = ElementIntegralPostprocessor::validParams();
24 10470 : params.addParam<unsigned int>(
25 : "fluid_component",
26 10470 : 0,
27 : "The index corresponding to the fluid component that this Postprocessor acts on");
28 10470 : params.addRequiredParam<UserObjectName>(
29 : "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
30 10470 : params.addParam<std::vector<unsigned int>>("phase",
31 : {},
32 : "The index of the fluid phase that this "
33 : "Postprocessor is restricted to. Multiple "
34 : "indices can be entered");
35 15705 : params.addRangeCheckedParam<Real>("saturation_threshold",
36 10470 : 1.0,
37 : "saturation_threshold >= 0 & saturation_threshold <= 1",
38 : "The saturation threshold below which the mass is calculated "
39 : "for a specific phase. Default is 1.0. Note: only one "
40 : "phase_index can be entered");
41 10470 : params.addParam<unsigned int>("kernel_variable_number",
42 10470 : 0,
43 : "The PorousFlow variable number (according to the dictator) of "
44 : "the fluid-mass kernel. This is required only in the unusual "
45 : "situation where a variety of different finite-element "
46 : "interpolation schemes are employed in the simulation");
47 10470 : params.addParam<std::string>(
48 : "base_name",
49 : "For non-mechanically-coupled systems with no TensorMechanics strain calculators, base_name "
50 : "need not be set. For mechanically-coupled systems, base_name should be the same base_name "
51 : "as given to the TensorMechanics object that computes strain, so that this Postprocessor can "
52 : "correctly account for changes in mesh volume. For non-mechanically-coupled systems, "
53 : "base_name should not be the base_name of any TensorMechanics strain calculators.");
54 5235 : params.set<bool>("use_displaced_mesh") = false;
55 5235 : params.suppressParameter<bool>("use_displaced_mesh");
56 5235 : params.addClassDescription("Calculates the mass of a fluid component in a region");
57 5235 : return params;
58 0 : }
59 :
60 : template <bool is_ad>
61 2829 : PorousFlowFluidMassTempl<is_ad>::PorousFlowFluidMassTempl(const InputParameters & parameters)
62 : : ElementIntegralPostprocessor(parameters),
63 :
64 2829 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
65 5658 : _fluid_component(getParam<unsigned int>("fluid_component")),
66 5658 : _phase_index(getParam<std::vector<unsigned int>>("phase")),
67 5662 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
68 2829 : _has_total_strain(hasMaterialProperty<RankTwoTensor>(_base_name + "total_strain")),
69 5658 : _total_strain(_has_total_strain
70 3159 : ? &getMaterialProperty<RankTwoTensor>(_base_name + "total_strain")
71 : : nullptr),
72 5658 : _porosity(getGenericMaterialProperty<Real, is_ad>("PorousFlow_porosity_nodal")),
73 5658 : _fluid_density(getGenericMaterialProperty<std::vector<Real>, is_ad>(
74 : "PorousFlow_fluid_phase_density_nodal")),
75 2829 : _fluid_saturation(
76 2829 : getGenericMaterialProperty<std::vector<Real>, is_ad>("PorousFlow_saturation_nodal")),
77 5658 : _mass_fraction(getGenericMaterialProperty<std::vector<std::vector<Real>>, is_ad>(
78 : "PorousFlow_mass_frac_nodal")),
79 5658 : _saturation_threshold(getParam<Real>("saturation_threshold")),
80 5658 : _var(getParam<unsigned>("kernel_variable_number") < _dictator.numVariables()
81 5656 : ? &_fe_problem.getStandardVariable(
82 2827 : _tid,
83 2827 : _dictator
84 8483 : .getCoupledStandardMooseVars()[getParam<unsigned>("kernel_variable_number")]
85 2827 : ->name())
86 2829 : : nullptr)
87 : {
88 2829 : const unsigned int num_phases = _dictator.numPhases();
89 2829 : const unsigned int num_components = _dictator.numComponents();
90 :
91 : // Check that the number of components entered is not greater than the total number of components
92 2829 : if (_fluid_component >= num_components)
93 0 : paramError(
94 : "fluid_component",
95 : "The Dictator proclaims that the number of components in this simulation is ",
96 : num_components,
97 : " whereas you have used a component index of ",
98 : _fluid_component,
99 : ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
100 :
101 : // Check that the number of phases entered is not more than the total possible phases
102 2829 : if (_phase_index.size() > num_phases)
103 2 : paramError("phase",
104 : "The Dictator decrees that the number of phases in this simulation is ",
105 : num_phases,
106 : " but you have entered ",
107 : _phase_index.size(),
108 : " phases.");
109 :
110 : // Check that kernel_variable_number is OK
111 8481 : if (getParam<unsigned>("kernel_variable_number") >= _dictator.numVariables())
112 2 : paramError("kernel_variable_number",
113 : "The Dictator pronounces that the number of PorousFlow variables is ",
114 2 : _dictator.numVariables(),
115 : ", however you have used ",
116 2 : getParam<unsigned>("kernel_variable_number"),
117 : ". This is an error");
118 :
119 : // Now that we know kernel_variable_number is OK, _var must be OK,
120 : // so ensure that reinit is called on _var:
121 5650 : addMooseVariableDependency(_var);
122 :
123 : // Also check that the phase indices entered are not greater than the number of phases
124 : // to avoid a segfault. Note that the input parser takes care of negative inputs so we
125 : // don't need to guard against them
126 2825 : if (!_phase_index.empty())
127 : {
128 854 : unsigned int max_phase_num = *std::max_element(_phase_index.begin(), _phase_index.end());
129 854 : if (max_phase_num > num_phases - 1)
130 2 : paramError("phase",
131 : "The Dictator proclaims that the phase index ",
132 : max_phase_num,
133 : " is greater than the largest phase index possible, which is ",
134 : num_phases - 1);
135 : }
136 :
137 : // Using saturation_threshold only makes sense for a specific phase_index
138 2823 : if (_saturation_threshold < 1.0 && _phase_index.size() != 1)
139 2 : paramError("saturation_threshold",
140 : "A single phase_index must be entered when prescribing a saturation_threshold");
141 :
142 : // If _phase_index is empty, create vector of all phase numbers to calculate mass over all phases
143 2821 : if (_phase_index.empty())
144 4511 : for (unsigned int i = 0; i < num_phases; ++i)
145 2540 : _phase_index.push_back(i);
146 :
147 : // Error if a strain base_name is provided but doesn't exist
148 5642 : if (parameters.isParamSetByUser("base_name") && !_has_total_strain)
149 2 : paramError("base_name", "A strain base_name ", _base_name, " does not exist");
150 2819 : }
151 :
152 : template <bool is_ad>
153 : Real
154 119973 : PorousFlowFluidMassTempl<is_ad>::computeIntegral()
155 : {
156 : Real sum = 0;
157 :
158 : // The use of _test in the loops below mean that the
159 : // integral is exactly the same as the one computed
160 : // by the PorousFlowMassTimeDerivative Kernel. Because that
161 : // Kernel is lumped, this Postprocessor also needs to
162 : // be lumped. Hence the use of the "nodal" Material
163 : // Properties
164 119973 : const VariableTestValue & test = _var->phi();
165 :
166 727569 : for (unsigned node = 0; node < test.size(); ++node)
167 : {
168 : Real nodal_volume = 0.0;
169 4538052 : for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
170 : {
171 3930456 : const Real n_v = _JxW[_qp] * _coord[_qp] * test[node][_qp];
172 3930456 : if (_has_total_strain)
173 51456 : nodal_volume += n_v * (1.0 + (*_total_strain)[_qp].trace());
174 : else
175 3879000 : nodal_volume += n_v;
176 : }
177 :
178 : Real mass = 0.0;
179 1310536 : for (auto ph : _phase_index)
180 : {
181 702940 : if (MetaPhysicL::raw_value(_fluid_saturation[node][ph]) <= _saturation_threshold)
182 702870 : mass += MetaPhysicL::raw_value(_fluid_density[node][ph] * _fluid_saturation[node][ph] *
183 702870 : _mass_fraction[node][ph][_fluid_component]);
184 : }
185 607596 : sum += nodal_volume * MetaPhysicL::raw_value(_porosity[node]) * mass;
186 : }
187 :
188 119973 : return sum;
189 : }
190 :
191 : template <bool is_ad>
192 : Real
193 0 : PorousFlowFluidMassTempl<is_ad>::computeQpIntegral()
194 : {
195 0 : return 0.0;
196 : }
197 :
198 : template class PorousFlowFluidMassTempl<false>;
199 : template class PorousFlowFluidMassTempl<true>;
|
