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 "PorousFlowPorosity.h"
11 :
12 : registerMooseObject("PorousFlowApp", PorousFlowPorosity);
13 :
14 : InputParameters
15 10181 : PorousFlowPorosity::validParams()
16 : {
17 10181 : InputParameters params = PorousFlowPorosityExponentialBase::validParams();
18 20362 : params.addParam<bool>(
19 20362 : "mechanical", false, "If true, porosity will be a function of total volumetric strain");
20 20362 : params.addParam<bool>(
21 20362 : "fluid", false, "If true, porosity will be a function of effective porepressure");
22 20362 : params.addParam<bool>("thermal", false, "If true, porosity will be a function of temperature");
23 20362 : params.addParam<bool>("chemical", false, "If true, porosity will be a function of precipitate");
24 20362 : params.addRequiredCoupledVar("porosity_zero",
25 : "The porosity at zero volumetric strain and "
26 : "reference temperature and reference effective "
27 : "porepressure and reference chemistry. This must be a real number "
28 : "or a constant monomial variable (not a linear lagrange or other "
29 : "type of variable)");
30 20362 : params.addParam<Real>("thermal_expansion_coeff",
31 : "Volumetric thermal expansion coefficient of the drained porous solid "
32 : "skeleton (only used if thermal=true)");
33 20362 : params.addRangeCheckedParam<Real>(
34 : "biot_coefficient", 1, "biot_coefficient>=0 & biot_coefficient<=1", "Biot coefficient");
35 20362 : params.addParam<Real>("biot_coefficient_prime",
36 : "Biot coefficient that appears in the term (biot_coefficient_prime - 1) * "
37 : "(P - reference_porepressure) / solid_bulk. If not provided, this "
38 : "defaults to the standard biot_coefficient");
39 20362 : params.addParam<MooseFunctorName>(
40 : "solid_bulk", "Bulk modulus of the drained porous solid skeleton (only used if fluid=true)");
41 20362 : params.addCoupledVar(
42 : "reference_temperature", 0.0, "Reference temperature (only used if thermal=true)");
43 20362 : params.addCoupledVar(
44 : "reference_porepressure", 0.0, "Reference porepressure (only used if fluid=true)");
45 20362 : params.addCoupledVar("reference_chemistry",
46 : "Reference values of the solid mineral concentrations "
47 : "(m^3(precipitate)/m^3(porous material)), entered as "
48 : "a vector (one value per mineral). (Only used if chemical=true)");
49 20362 : params.addCoupledVar(
50 : "initial_mineral_concentrations",
51 : "Initial mineral concentrations (m^3(precipitate)/m^3(porous material)), entered as "
52 : "a vector (one value per mineral). (Only used if chemical=true)");
53 20362 : params.addParam<std::vector<Real>>("chemical_weights",
54 : "When chemical=true, porosity is a linear combination of the "
55 : "solid mineral concentrations multiplied by these weights. "
56 : "Default=1 for all minerals.");
57 10181 : params.addClassDescription("This Material calculates the porosity PorousFlow simulations");
58 10181 : return params;
59 0 : }
60 :
61 8009 : PorousFlowPorosity::PorousFlowPorosity(const InputParameters & parameters)
62 : : PorousFlowPorosityExponentialBase(parameters),
63 :
64 8009 : _mechanical(getParam<bool>("mechanical")),
65 16018 : _fluid(getParam<bool>("fluid")),
66 16018 : _thermal(getParam<bool>("thermal")),
67 16018 : _chemical(getParam<bool>("chemical")),
68 8009 : _phi0(coupledValue("porosity_zero")),
69 16018 : _biot(getParam<Real>("biot_coefficient")),
70 17948 : _exp_coeff(isParamValid("thermal_expansion_coeff") ? getParam<Real>("thermal_expansion_coeff")
71 : : 0.0),
72 21914 : _solid_bulk(isParamValid("solid_bulk") ? &(getFunctor<Real>("solid_bulk")) : nullptr),
73 24159 : _coeff(isParamValid("biot_coefficient_prime") ? (getParam<Real>("biot_coefficient_prime") - 1.0)
74 7877 : : (_biot - 1.0)),
75 :
76 8009 : _t_reference(_nodal_material ? coupledDofValues("reference_temperature")
77 12049 : : coupledValue("reference_temperature")),
78 8009 : _p_reference(_nodal_material ? coupledDofValues("reference_porepressure")
79 12049 : : coupledValue("reference_porepressure")),
80 8009 : _num_c_ref(coupledComponents("reference_chemistry")),
81 8009 : _c_reference(_num_c_ref),
82 8009 : _num_initial_c(coupledComponents("initial_mineral_concentrations")),
83 8009 : _initial_c(_num_initial_c),
84 24060 : _c_weights(isParamValid("chemical_weights") ? getParam<std::vector<Real>>("chemical_weights")
85 8009 : : std::vector<Real>(_num_c_ref, 1.0)),
86 :
87 8009 : _porosity_old(_chemical ? (_nodal_material
88 103 : ? &getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")
89 8149 : : &getMaterialPropertyOld<Real>("PorousFlow_porosity_qp"))
90 : : nullptr),
91 11091 : _vol_strain_qp(_mechanical ? &getMaterialProperty<Real>("PorousFlow_total_volumetric_strain_qp")
92 : : nullptr),
93 11091 : _dvol_strain_qp_dvar(_mechanical ? &getMaterialProperty<std::vector<RealGradient>>(
94 : "dPorousFlow_total_volumetric_strain_qp_dvar")
95 : : nullptr),
96 :
97 8009 : _pf(_fluid ? (_nodal_material
98 2950 : ? &getMaterialProperty<Real>("PorousFlow_effective_fluid_pressure_nodal")
99 10999 : : &getMaterialProperty<Real>("PorousFlow_effective_fluid_pressure_qp"))
100 : : nullptr),
101 8009 : _dpf_dvar(_fluid ? (_nodal_material ? &getMaterialProperty<std::vector<Real>>(
102 : "dPorousFlow_effective_fluid_pressure_nodal_dvar")
103 10999 : : &getMaterialProperty<std::vector<Real>>(
104 : "dPorousFlow_effective_fluid_pressure_qp_dvar"))
105 : : nullptr),
106 :
107 16018 : _temperature(_thermal
108 8009 : ? (_nodal_material ? &getMaterialProperty<Real>("PorousFlow_temperature_nodal")
109 8881 : : &getMaterialProperty<Real>("PorousFlow_temperature_qp"))
110 : : nullptr),
111 8009 : _dtemperature_dvar(
112 8009 : _thermal
113 8009 : ? (_nodal_material
114 967 : ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
115 8881 : : &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar"))
116 : : nullptr),
117 :
118 8009 : _mineral_conc_old(_chemical ? (_nodal_material ? &getMaterialPropertyOld<std::vector<Real>>(
119 : "PorousFlow_mineral_concentration_nodal")
120 8149 : : &getMaterialPropertyOld<std::vector<Real>>(
121 : "PorousFlow_mineral_concentration_qp"))
122 : : nullptr),
123 8009 : _reaction_rate(_chemical ? (_nodal_material ? &getMaterialProperty<std::vector<Real>>(
124 : "PorousFlow_mineral_reaction_rate_nodal")
125 8149 : : &getMaterialProperty<std::vector<Real>>(
126 : "PorousFlow_mineral_reaction_rate_qp"))
127 : : nullptr),
128 8009 : _dreaction_rate_dvar(_chemical ? (_nodal_material
129 103 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
130 : "dPorousFlow_mineral_reaction_rate_nodal_dvar")
131 8149 : : &getMaterialProperty<std::vector<std::vector<Real>>>(
132 : "dPorousFlow_mineral_reaction_rate_qp_dvar"))
133 : : nullptr),
134 8009 : _aq_ph(_dictator.aqueousPhaseNumber()),
135 16018 : _saturation(_chemical
136 8009 : ? (_nodal_material
137 103 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
138 8149 : : &getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_qp"))
139 : : nullptr),
140 16018 : _dsaturation_dvar(_chemical
141 8009 : ? (_nodal_material ? &getMaterialProperty<std::vector<std::vector<Real>>>(
142 : "dPorousFlow_saturation_nodal_dvar")
143 8149 : : &getMaterialProperty<std::vector<std::vector<Real>>>(
144 : "dPorousFlow_saturation_qp_dvar"))
145 8009 : : nullptr)
146 : {
147 10908 : if (_thermal && !isParamValid("thermal_expansion_coeff"))
148 2 : mooseError("PorousFlowPorosity: When thermal=true you must provide a thermal_expansion_coeff");
149 8007 : if (_fluid && !_solid_bulk)
150 2 : mooseError("PorousFlowPorosity: When fluid=true you must provide a solid_bulk");
151 8005 : if (_chemical && _num_c_ref != _dictator.numAqueousKinetic())
152 4 : mooseError("PorousFlowPorosity: When chemical=true you must provide the reference_chemistry "
153 : "values. The Dictator proclaims there should be ",
154 2 : _dictator.numAqueousKinetic(),
155 : " of these");
156 8003 : if (_chemical && _num_initial_c != _dictator.numAqueousKinetic())
157 4 : mooseError("PorousFlowPorosity: When chemical=true you must provide the "
158 : "initial_mineral_concentrations. "
159 : "The Dictator proclaims there should be ",
160 2 : _dictator.numAqueousKinetic(),
161 : " of these");
162 8001 : if (_chemical && _c_weights.size() != _dictator.numAqueousKinetic())
163 0 : mooseError(
164 : "PorousFlowPorosity: When chemical=true you must provde the correct number of "
165 : "chemical_weights (which the Dictator knows is ",
166 0 : _dictator.numAqueousKinetic(),
167 : ") or do not provide any chemical_weights and use the default value of 1 for each mineral");
168 :
169 8133 : for (unsigned i = 0; i < _num_c_ref; ++i)
170 : {
171 198 : _c_reference[i] = (_nodal_material ? &coupledDofValues("reference_chemistry", i)
172 198 : : &coupledValue("reference_chemistry", i));
173 264 : _initial_c[i] = (_nodal_material ? &coupledDofValues("initial_mineral_concentrations", i)
174 198 : : &coupledValue("initial_mineral_concentrations", i));
175 : }
176 8001 : }
177 :
178 : Real
179 9588449 : PorousFlowPorosity::atNegInfinityQp() const
180 : {
181 : /*
182 : *
183 : * Note the use of the OLD value of porosity here.
184 : * This strategy, which breaks the cyclic dependency between porosity
185 : * and mineral concentration, is used in
186 : * Kernel: PorousFlowPreDis
187 : * Material: PorousFlowPorosity
188 : * Material: PorousFlowAqueousPreDisChemistry
189 : * Material: PorousFlowAqueousPreDisMineral
190 : *
191 : */
192 9588449 : Real result = _biot;
193 9588449 : if (_chemical)
194 : {
195 652 : if (_t_step == 0 && !_app.isRestarting())
196 132 : for (unsigned i = 0; i < _num_c_ref; ++i)
197 76 : result -= _c_weights[i] * (*_initial_c[i])[_qp];
198 : else
199 666 : for (unsigned i = 0; i < _num_c_ref; ++i)
200 368 : result -= _c_weights[i] * ((*_mineral_conc_old)[_qp][i] + _dt * (*_porosity_old)[_qp] *
201 368 : (*_saturation)[_qp][_aq_ph] *
202 368 : (*_reaction_rate)[_qp][i]);
203 : }
204 9588449 : return result;
205 : }
206 :
207 : Real
208 22791972 : PorousFlowPorosity::datNegInfinityQp(unsigned pvar) const
209 : {
210 : Real result = 0.0;
211 22791972 : if (_chemical && (_t_step >= 1 || _app.isRestarting()))
212 654 : for (unsigned i = 0; i < _num_c_ref; ++i)
213 362 : result -= _c_weights[i] * _dt * (*_porosity_old)[_qp] *
214 362 : ((*_saturation)[_qp][_aq_ph] * (*_dreaction_rate_dvar)[_qp][i][pvar] +
215 362 : (*_dsaturation_dvar)[_qp][_aq_ph][pvar] * (*_reaction_rate)[_qp][i]);
216 22791972 : return result;
217 : }
218 :
219 : Real
220 9588449 : PorousFlowPorosity::atZeroQp() const
221 : {
222 : // note the [0] below: _phi0 is a constant monomial and we use [0] regardless of _nodal_material
223 9588449 : Real result = _phi0[0];
224 9588449 : if (_chemical)
225 : {
226 652 : if (_t_step == 0 && !_app.isRestarting())
227 132 : for (unsigned i = 0; i < _num_c_ref; ++i)
228 76 : result -= _c_weights[i] * ((*_initial_c[i])[_qp] - (*_c_reference[i])[_qp]);
229 : else
230 666 : for (unsigned i = 0; i < _num_c_ref; ++i)
231 368 : result -= _c_weights[i] * ((*_mineral_conc_old)[_qp][i] +
232 368 : _dt * (*_porosity_old)[_qp] * (*_saturation)[_qp][_aq_ph] *
233 368 : (*_reaction_rate)[_qp][i] -
234 368 : (*_c_reference[i])[_qp]);
235 : }
236 9588449 : return result;
237 : }
238 :
239 : Real
240 22791972 : PorousFlowPorosity::datZeroQp(unsigned pvar) const
241 : {
242 : Real result = 0.0;
243 22791972 : if (_chemical && (_t_step >= 1 || _app.isRestarting()))
244 654 : for (unsigned i = 0; i < _num_c_ref; ++i)
245 362 : result -= _c_weights[i] * _dt * (*_porosity_old)[_qp] *
246 362 : ((*_saturation)[_qp][_aq_ph] * (*_dreaction_rate_dvar)[_qp][i][pvar] +
247 362 : (*_dsaturation_dvar)[_qp][_aq_ph][pvar] * (*_reaction_rate)[_qp][i]);
248 22791972 : return result;
249 : }
250 :
251 : Real
252 9588449 : PorousFlowPorosity::decayQp() const
253 : {
254 : Real result = 0.0;
255 :
256 9588449 : if (_thermal)
257 107536 : result += _exp_coeff * ((*_temperature)[_qp] - _t_reference[_qp]);
258 :
259 9588449 : if (_fluid)
260 : {
261 : Real solid_bulk;
262 : // Using Qp 0 can leverage the functor caching
263 : // TODO: Find a way to effectively use subdomain-constant-ness
264 2895096 : unsigned int qp_used = (_constant_option == ConstantTypeEnum::NONE) ? _qp : 0;
265 2895096 : if (_nodal_material)
266 : {
267 1469216 : const std::set<SubdomainID> subdomain_set = {_current_elem->subdomain_id()};
268 1469216 : const Moose::NodeArg space_arg = {_current_elem->node_ptr(qp_used), &subdomain_set};
269 1469216 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
270 : }
271 1425880 : else if (_bnd)
272 : {
273 : const Moose::ElemSideQpArg space_arg = {
274 24924 : _current_elem, _current_side, qp_used, _qrule, _q_point[qp_used]};
275 24924 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
276 : }
277 : else
278 : {
279 1400956 : const Moose::ElemQpArg space_arg = {_current_elem, qp_used, _qrule, _q_point[qp_used]};
280 1400956 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
281 : }
282 2895096 : if (solid_bulk <= 0)
283 0 : mooseError("PorousFlowPorosity: solid_bulk must be larger than Zero");
284 2895096 : result += _coeff / solid_bulk * ((*_pf)[_qp] - _p_reference[_qp]);
285 : }
286 :
287 9588449 : if (_mechanical)
288 : {
289 : // Note that in the following _strain[_qp] is evaluated at q quadpoint
290 : // So _porosity_nodal[_qp], which should be the nodal value of porosity
291 : // actually uses the strain at a quadpoint. This
292 : // is OK for LINEAR elements, as strain is constant over the element anyway.
293 : const unsigned qp_to_use =
294 2903266 : (_nodal_material && (_bnd || _strain_at_nearest_qp) ? nearestQP(_qp) : _qp);
295 2903266 : result += -(*_vol_strain_qp)[qp_to_use];
296 : }
297 :
298 9588449 : return result;
299 : }
300 :
301 : Real
302 22791972 : PorousFlowPorosity::ddecayQp_dvar(unsigned pvar) const
303 : {
304 : Real result = 0.0;
305 :
306 22791972 : if (_thermal)
307 494072 : result += _exp_coeff * (*_dtemperature_dvar)[_qp][pvar];
308 :
309 22791972 : if (_fluid)
310 : {
311 : Real solid_bulk;
312 : // Using Qp 0 can leverage the functor caching
313 : // TODO: Find a way to effectively use subdomain-constant-ness
314 11547142 : unsigned int qp_used = (_constant_option == ConstantTypeEnum::NONE) ? _qp : 0;
315 11547142 : if (_nodal_material)
316 : {
317 5891048 : const std::set<SubdomainID> subdomain_set = {_current_elem->subdomain_id()};
318 5891048 : const Moose::NodeArg space_arg = {_current_elem->node_ptr(qp_used), &subdomain_set};
319 5891048 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
320 : }
321 5656094 : else if (_bnd)
322 : {
323 : const Moose::ElemSideQpArg space_arg = {
324 98576 : _current_elem, _current_side, qp_used, _qrule, _q_point[qp_used]};
325 98576 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
326 : }
327 : else
328 : {
329 5557518 : const Moose::ElemQpArg space_arg = {_current_elem, qp_used, _qrule, _q_point[qp_used]};
330 5557518 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
331 : }
332 11547142 : if (solid_bulk <= 0)
333 0 : mooseError("PorousFlowPorosity: solid_bulk must be larger than Zero.");
334 11547142 : result += _coeff / solid_bulk * (*_dpf_dvar)[_qp][pvar];
335 : }
336 :
337 22791972 : return result;
338 : }
339 :
340 : RealGradient
341 22791972 : PorousFlowPorosity::ddecayQp_dgradvar(unsigned pvar) const
342 : {
343 : RealGradient result(0.0, 0.0, 0.0);
344 22791972 : if (_mechanical)
345 : {
346 : const unsigned qp_to_use =
347 11590232 : (_nodal_material && (_bnd || _strain_at_nearest_qp) ? nearestQP(_qp) : _qp);
348 11590232 : result += -(*_dvol_strain_qp_dvar)[qp_to_use][pvar];
349 : }
350 22791972 : return result;
351 : }
|
