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 5159 : PorousFlowPorosity::validParams()
16 : {
17 5159 : InputParameters params = PorousFlowPorosityExponentialBase::validParams();
18 10318 : params.addParam<bool>(
19 10318 : "mechanical", false, "If true, porosity will be a function of total volumetric strain");
20 10318 : params.addParam<bool>(
21 10318 : "fluid", false, "If true, porosity will be a function of effective porepressure");
22 10318 : params.addParam<bool>("thermal", false, "If true, porosity will be a function of temperature");
23 10318 : params.addParam<bool>("chemical", false, "If true, porosity will be a function of precipitate");
24 10318 : 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 10318 : params.addParam<Real>("thermal_expansion_coeff",
31 : "Volumetric thermal expansion coefficient of the drained porous solid "
32 : "skeleton (only used if thermal=true)");
33 10318 : params.addRangeCheckedParam<Real>(
34 : "biot_coefficient", 1, "biot_coefficient>=0 & biot_coefficient<=1", "Biot coefficient");
35 10318 : 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 10318 : params.addParam<MooseFunctorName>(
40 : "solid_bulk", "Bulk modulus of the drained porous solid skeleton (only used if fluid=true)");
41 10318 : params.addCoupledVar(
42 : "reference_temperature", 0.0, "Reference temperature (only used if thermal=true)");
43 10318 : params.addCoupledVar(
44 : "reference_porepressure", 0.0, "Reference porepressure (only used if fluid=true)");
45 10318 : 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 10318 : 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 10318 : 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 5159 : params.addClassDescription("This Material calculates the porosity PorousFlow simulations");
58 5159 : return params;
59 0 : }
60 :
61 4007 : PorousFlowPorosity::PorousFlowPorosity(const InputParameters & parameters)
62 : : PorousFlowPorosityExponentialBase(parameters),
63 :
64 4007 : _mechanical(getParam<bool>("mechanical")),
65 8014 : _fluid(getParam<bool>("fluid")),
66 8014 : _thermal(getParam<bool>("thermal")),
67 8014 : _chemical(getParam<bool>("chemical")),
68 4007 : _phi0(coupledValue("porosity_zero")),
69 8014 : _biot(getParam<Real>("biot_coefficient")),
70 8984 : _exp_coeff(isParamValid("thermal_expansion_coeff") ? getParam<Real>("thermal_expansion_coeff")
71 : : 0.0),
72 10934 : _solid_bulk(isParamValid("solid_bulk") ? &(getFunctor<Real>("solid_bulk")) : nullptr),
73 12081 : _coeff(isParamValid("biot_coefficient_prime") ? (getParam<Real>("biot_coefficient_prime") - 1.0)
74 3947 : : (_biot - 1.0)),
75 :
76 4007 : _t_reference(_nodal_material ? coupledDofValues("reference_temperature")
77 5974 : : coupledValue("reference_temperature")),
78 4007 : _p_reference(_nodal_material ? coupledDofValues("reference_porepressure")
79 5974 : : coupledValue("reference_porepressure")),
80 4007 : _num_c_ref(coupledComponents("reference_chemistry")),
81 4007 : _c_reference(_num_c_ref),
82 4007 : _num_initial_c(coupledComponents("initial_mineral_concentrations")),
83 4007 : _initial_c(_num_initial_c),
84 12042 : _c_weights(isParamValid("chemical_weights") ? getParam<std::vector<Real>>("chemical_weights")
85 4007 : : std::vector<Real>(_num_c_ref, 1.0)),
86 :
87 4007 : _porosity_old(_chemical ? (_nodal_material
88 55 : ? &getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")
89 4075 : : &getMaterialPropertyOld<Real>("PorousFlow_porosity_qp"))
90 : : nullptr),
91 5529 : _vol_strain_qp(_mechanical ? &getMaterialProperty<Real>("PorousFlow_total_volumetric_strain_qp")
92 : : nullptr),
93 5529 : _dvol_strain_qp_dvar(_mechanical ? &getMaterialProperty<std::vector<RealGradient>>(
94 : "dPorousFlow_total_volumetric_strain_qp_dvar")
95 : : nullptr),
96 :
97 4007 : _pf(_fluid ? (_nodal_material
98 1462 : ? &getMaterialProperty<Real>("PorousFlow_effective_fluid_pressure_nodal")
99 5461 : : &getMaterialProperty<Real>("PorousFlow_effective_fluid_pressure_qp"))
100 : : nullptr),
101 4007 : _dpf_dvar(_fluid ? (_nodal_material ? &getMaterialProperty<std::vector<Real>>(
102 : "dPorousFlow_effective_fluid_pressure_nodal_dvar")
103 5461 : : &getMaterialProperty<std::vector<Real>>(
104 : "dPorousFlow_effective_fluid_pressure_qp_dvar"))
105 : : nullptr),
106 :
107 8014 : _temperature(_thermal
108 4007 : ? (_nodal_material ? &getMaterialProperty<Real>("PorousFlow_temperature_nodal")
109 4423 : : &getMaterialProperty<Real>("PorousFlow_temperature_qp"))
110 : : nullptr),
111 4007 : _dtemperature_dvar(
112 4007 : _thermal
113 4007 : ? (_nodal_material
114 487 : ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
115 4423 : : &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar"))
116 : : nullptr),
117 :
118 4007 : _mineral_conc_old(_chemical ? (_nodal_material ? &getMaterialPropertyOld<std::vector<Real>>(
119 : "PorousFlow_mineral_concentration_nodal")
120 4075 : : &getMaterialPropertyOld<std::vector<Real>>(
121 : "PorousFlow_mineral_concentration_qp"))
122 : : nullptr),
123 4007 : _reaction_rate(_chemical ? (_nodal_material ? &getMaterialProperty<std::vector<Real>>(
124 : "PorousFlow_mineral_reaction_rate_nodal")
125 4075 : : &getMaterialProperty<std::vector<Real>>(
126 : "PorousFlow_mineral_reaction_rate_qp"))
127 : : nullptr),
128 4007 : _dreaction_rate_dvar(_chemical ? (_nodal_material
129 55 : ? &getMaterialProperty<std::vector<std::vector<Real>>>(
130 : "dPorousFlow_mineral_reaction_rate_nodal_dvar")
131 4075 : : &getMaterialProperty<std::vector<std::vector<Real>>>(
132 : "dPorousFlow_mineral_reaction_rate_qp_dvar"))
133 : : nullptr),
134 4007 : _aq_ph(_dictator.aqueousPhaseNumber()),
135 8014 : _saturation(_chemical
136 4007 : ? (_nodal_material
137 55 : ? &getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
138 4075 : : &getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_qp"))
139 : : nullptr),
140 8014 : _dsaturation_dvar(_chemical
141 4007 : ? (_nodal_material ? &getMaterialProperty<std::vector<std::vector<Real>>>(
142 : "dPorousFlow_saturation_nodal_dvar")
143 4075 : : &getMaterialProperty<std::vector<std::vector<Real>>>(
144 : "dPorousFlow_saturation_qp_dvar"))
145 4007 : : nullptr)
146 : {
147 5466 : if (_thermal && !isParamValid("thermal_expansion_coeff"))
148 2 : mooseError("PorousFlowPorosity: When thermal=true you must provide a thermal_expansion_coeff");
149 4005 : if (_fluid && !_solid_bulk)
150 2 : mooseError("PorousFlowPorosity: When fluid=true you must provide a solid_bulk");
151 4003 : 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 4001 : 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 3999 : 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 4071 : for (unsigned i = 0; i < _num_c_ref; ++i)
170 : {
171 102 : _c_reference[i] = (_nodal_material ? &coupledDofValues("reference_chemistry", i)
172 102 : : &coupledValue("reference_chemistry", i));
173 144 : _initial_c[i] = (_nodal_material ? &coupledDofValues("initial_mineral_concentrations", i)
174 102 : : &coupledValue("initial_mineral_concentrations", i));
175 : }
176 3999 : }
177 :
178 : Real
179 7254964 : 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 7254964 : Real result = _biot;
193 7254964 : if (_chemical)
194 : {
195 456 : if (_t_step == 0 && !_app.isRestarting())
196 96 : for (unsigned i = 0; i < _num_c_ref; ++i)
197 56 : result -= _c_weights[i] * (*_initial_c[i])[_qp];
198 : else
199 472 : for (unsigned i = 0; i < _num_c_ref; ++i)
200 264 : result -= _c_weights[i] * ((*_mineral_conc_old)[_qp][i] + _dt * (*_porosity_old)[_qp] *
201 264 : (*_saturation)[_qp][_aq_ph] *
202 264 : (*_reaction_rate)[_qp][i]);
203 : }
204 7254964 : return result;
205 : }
206 :
207 : Real
208 16132469 : PorousFlowPorosity::datNegInfinityQp(unsigned pvar) const
209 : {
210 : Real result = 0.0;
211 16132469 : if (_chemical && (_t_step >= 1 || _app.isRestarting()))
212 468 : for (unsigned i = 0; i < _num_c_ref; ++i)
213 262 : result -= _c_weights[i] * _dt * (*_porosity_old)[_qp] *
214 262 : ((*_saturation)[_qp][_aq_ph] * (*_dreaction_rate_dvar)[_qp][i][pvar] +
215 262 : (*_dsaturation_dvar)[_qp][_aq_ph][pvar] * (*_reaction_rate)[_qp][i]);
216 16132469 : return result;
217 : }
218 :
219 : Real
220 7254964 : PorousFlowPorosity::atZeroQp() const
221 : {
222 : // note the [0] below: _phi0 is a constant monomial and we use [0] regardless of _nodal_material
223 7254964 : Real result = _phi0[0];
224 7254964 : if (_chemical)
225 : {
226 456 : if (_t_step == 0 && !_app.isRestarting())
227 96 : for (unsigned i = 0; i < _num_c_ref; ++i)
228 56 : result -= _c_weights[i] * ((*_initial_c[i])[_qp] - (*_c_reference[i])[_qp]);
229 : else
230 472 : for (unsigned i = 0; i < _num_c_ref; ++i)
231 264 : result -= _c_weights[i] * ((*_mineral_conc_old)[_qp][i] +
232 264 : _dt * (*_porosity_old)[_qp] * (*_saturation)[_qp][_aq_ph] *
233 264 : (*_reaction_rate)[_qp][i] -
234 264 : (*_c_reference[i])[_qp]);
235 : }
236 7254964 : return result;
237 : }
238 :
239 : Real
240 16132469 : PorousFlowPorosity::datZeroQp(unsigned pvar) const
241 : {
242 : Real result = 0.0;
243 16132469 : if (_chemical && (_t_step >= 1 || _app.isRestarting()))
244 468 : for (unsigned i = 0; i < _num_c_ref; ++i)
245 262 : result -= _c_weights[i] * _dt * (*_porosity_old)[_qp] *
246 262 : ((*_saturation)[_qp][_aq_ph] * (*_dreaction_rate_dvar)[_qp][i][pvar] +
247 262 : (*_dsaturation_dvar)[_qp][_aq_ph][pvar] * (*_reaction_rate)[_qp][i]);
248 16132469 : return result;
249 : }
250 :
251 : Real
252 7254964 : PorousFlowPorosity::decayQp() const
253 : {
254 : Real result = 0.0;
255 :
256 7254964 : if (_thermal)
257 80248 : result += _exp_coeff * ((*_temperature)[_qp] - _t_reference[_qp]);
258 :
259 7254964 : 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 1892566 : unsigned int qp_used = (_constant_option == ConstantTypeEnum::NONE) ? _qp : 0;
265 1892566 : if (_nodal_material)
266 : {
267 954224 : const std::set<SubdomainID> subdomain_set = {_current_elem->subdomain_id()};
268 954224 : const Moose::NodeArg space_arg = {_current_elem->node_ptr(qp_used), &subdomain_set};
269 954224 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
270 : }
271 938342 : else if (_bnd)
272 : {
273 : const Moose::ElemSideQpArg space_arg = {
274 536 : _current_elem, _current_side, qp_used, _qrule, _q_point[qp_used]};
275 536 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
276 : }
277 : else
278 : {
279 937806 : const Moose::ElemQpArg space_arg = {_current_elem, qp_used, _qrule, _q_point[qp_used]};
280 937806 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
281 : }
282 1892566 : if (solid_bulk <= 0)
283 0 : mooseError("PorousFlowPorosity: solid_bulk must be larger than Zero");
284 1892566 : result += _coeff / solid_bulk * ((*_pf)[_qp] - _p_reference[_qp]);
285 : }
286 :
287 7254964 : 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 1899036 : (_nodal_material && (_bnd || _strain_at_nearest_qp) ? nearestQP(_qp) : _qp);
295 1899036 : result += -(*_vol_strain_qp)[qp_to_use];
296 : }
297 :
298 7254964 : return result;
299 : }
300 :
301 : Real
302 16132469 : PorousFlowPorosity::ddecayQp_dvar(unsigned pvar) const
303 : {
304 : Real result = 0.0;
305 :
306 16132469 : if (_thermal)
307 367368 : result += _exp_coeff * (*_dtemperature_dvar)[_qp][pvar];
308 :
309 16132469 : 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 7555464 : unsigned int qp_used = (_constant_option == ConstantTypeEnum::NONE) ? _qp : 0;
315 7555464 : if (_nodal_material)
316 : {
317 3833392 : const std::set<SubdomainID> subdomain_set = {_current_elem->subdomain_id()};
318 3833392 : const Moose::NodeArg space_arg = {_current_elem->node_ptr(qp_used), &subdomain_set};
319 3833392 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
320 : }
321 3722072 : else if (_bnd)
322 : {
323 : const Moose::ElemSideQpArg space_arg = {
324 1792 : _current_elem, _current_side, qp_used, _qrule, _q_point[qp_used]};
325 1792 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
326 : }
327 : else
328 : {
329 3720280 : const Moose::ElemQpArg space_arg = {_current_elem, qp_used, _qrule, _q_point[qp_used]};
330 3720280 : solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
331 : }
332 7555464 : if (solid_bulk <= 0)
333 0 : mooseError("PorousFlowPorosity: solid_bulk must be larger than Zero.");
334 7555464 : result += _coeff / solid_bulk * (*_dpf_dvar)[_qp][pvar];
335 : }
336 :
337 16132469 : return result;
338 : }
339 :
340 : RealGradient
341 16132469 : PorousFlowPorosity::ddecayQp_dgradvar(unsigned pvar) const
342 : {
343 : RealGradient result(0.0, 0.0, 0.0);
344 16132469 : if (_mechanical)
345 : {
346 : const unsigned qp_to_use =
347 7589884 : (_nodal_material && (_bnd || _strain_at_nearest_qp) ? nearestQP(_qp) : _qp);
348 7589884 : result += -(*_dvol_strain_qp_dvar)[qp_to_use][pvar];
349 : }
350 16132469 : return result;
351 : }
|
