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 "PorousFlowHysteresisOrder.h"
11 : #include "Conversion.h"
12 :
13 : registerMooseObject("PorousFlowApp", PorousFlowHysteresisOrder);
14 :
15 : InputParameters
16 5954 : PorousFlowHysteresisOrder::validParams()
17 : {
18 5954 : InputParameters params = PorousFlowMaterial::validParams();
19 11908 : params.addPrivateParam<std::string>("pf_material_type", "hysteresis_order");
20 11908 : params.addParam<unsigned>("liquid_phase", 0, "The phase number of the liquid phase");
21 11908 : params.addParam<unsigned>(
22 : "initial_order",
23 11908 : 0,
24 : "The initial order. 0 means the simulation will begin on the primary drying curve. 1 means "
25 : "the simulation will begin on the first-order wetting curve. 2 means the simulation will "
26 : "begin on the second-order drying curve, etc.");
27 5954 : params.addParam<std::vector<Real>>(
28 : "previous_turning_points",
29 5954 : std::vector<Real>(),
30 : "The turning points (liquid saturation values) at occured prior to the simulation. There "
31 : "must be exactly initial_order of these values. The first value is the liquid saturation at "
32 : "the transition from primary-drying to first-order-wetting; the second value is the liquid "
33 : "saturation at the transition from first-order-wetting to second-order-drying; and so on. "
34 : "You must ensure that your initial saturation field is commensurate with the "
35 : "previous_turning_points.");
36 5954 : params.addClassDescription("Computes hysteresis order for use in hysteretic capillary pressures "
37 : "and relative permeabilities");
38 5954 : return params;
39 0 : }
40 :
41 4579 : PorousFlowHysteresisOrder::PorousFlowHysteresisOrder(const InputParameters & parameters)
42 : : DerivativeMaterialInterface<PorousFlowMaterial>(parameters),
43 4579 : _liquid_ph_num(getParam<unsigned>("liquid_phase")),
44 4579 : _liquid_phase(Moose::stringify(_liquid_ph_num)),
45 9158 : _initial_order(getParam<unsigned>("initial_order")),
46 9158 : _previous_turning_points(getParam<std::vector<Real>>("previous_turning_points")),
47 9158 : _sat_old(_nodal_material
48 4579 : ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")
49 12291 : : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_qp")),
50 9158 : _sat_older(_nodal_material
51 4579 : ? getMaterialPropertyOlder<std::vector<Real>>("PorousFlow_saturation_nodal")
52 12291 : : getMaterialPropertyOlder<std::vector<Real>>("PorousFlow_saturation_qp")),
53 4579 : _hys_order(_nodal_material ? declareProperty<unsigned>("PorousFlow_hysteresis_order_nodal")
54 8435 : : declareProperty<unsigned>("PorousFlow_hysteresis_order_qp")),
55 9158 : _hys_order_old(_nodal_material
56 4579 : ? getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_nodal")
57 12291 : : getMaterialPropertyOld<unsigned>("PorousFlow_hysteresis_order_qp")),
58 9158 : _hys_sat_tps(_nodal_material
59 4579 : ? declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
60 : "PorousFlow_hysteresis_saturation_tps_nodal")
61 8435 : : declareProperty<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
62 : "PorousFlow_hysteresis_saturation_tps_qp")),
63 4579 : _hys_sat_tps_old(
64 4579 : _nodal_material
65 4579 : ? getMaterialPropertyOld<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
66 : "PorousFlow_hysteresis_saturation_tps_nodal")
67 12291 : : getMaterialPropertyOld<std::array<Real, PorousFlowConstants::MAX_HYSTERESIS_ORDER>>(
68 4579 : "PorousFlow_hysteresis_saturation_tps_qp"))
69 : {
70 4579 : if (_liquid_ph_num >= _dictator.numPhases())
71 2 : paramError("liquid_phase",
72 : "The Dictator proclaims that the number of fluid phases is ",
73 2 : _dictator.numPhases(),
74 : " while you have foolishly entered ",
75 2 : _liquid_ph_num,
76 : " for your liquid phase number. Remember that indexing starts at 0. Be aware that "
77 : "the Dictator does not tolerate mistakes.");
78 4577 : if (_initial_order > PorousFlowConstants::MAX_HYSTERESIS_ORDER)
79 2 : paramError("initial_order",
80 : "The initial order must be less than the max order, which is hard-coded to ",
81 : PorousFlowConstants::MAX_HYSTERESIS_ORDER);
82 4575 : if (_previous_turning_points.size() != _initial_order)
83 2 : paramError("previous_turning_points",
84 : "initial_order is ",
85 : _initial_order,
86 : " so there must be this many previous_turning_points");
87 9977 : for (const auto & tp : _previous_turning_points)
88 5408 : if (tp < 0.0 || tp > 1.0)
89 4 : paramError("previous_turning_points",
90 : "Each previous_turning_point must lie in the range [0.0, 1.0]");
91 : // check that the previous_turning_points are ordered correctly
92 : Real previous_tp = 1.0;
93 : Real previous_jump = std::numeric_limits<Real>::max();
94 : bool drying = true;
95 9967 : for (unsigned i = 0; i < _initial_order; ++i)
96 : {
97 5404 : const Real this_jump = _previous_turning_points[i] - previous_tp;
98 5404 : if (drying)
99 : {
100 3628 : if (this_jump >= 0)
101 2 : paramError("previous_turning_points",
102 : "The previous_turning_points vector, {a, b, c, ...}, must be obey the following "
103 : "order: a < 1; b > a; a < c < b, etc");
104 : }
105 : else
106 : {
107 1776 : if (this_jump <= 0)
108 2 : paramError("previous_turning_points",
109 : "The previous_turning_points vector, {a, b, c, ...}, must be obey the following "
110 : "order: a < 1; b > a; a < c < b, etc");
111 : }
112 5400 : if (std::abs(this_jump) >= previous_jump)
113 2 : paramError("previous_turning_points",
114 : "The previous_turning_points vector, {a, b, c, ...}, must be obey the following "
115 : "order: a < 1; b > a; a < c < b, etc");
116 5398 : drying = !drying;
117 : previous_tp = _previous_turning_points[i];
118 : previous_jump = std::abs(this_jump);
119 : }
120 4563 : }
121 :
122 : void
123 213836 : PorousFlowHysteresisOrder::initQpStatefulProperties()
124 : {
125 213836 : _hys_order[_qp] = _initial_order;
126 482224 : for (unsigned i = 0; i < _initial_order; ++i)
127 268388 : _hys_sat_tps[_qp].at(i) = _previous_turning_points[i];
128 213836 : }
129 :
130 : void
131 663882 : PorousFlowHysteresisOrder::computeQpProperties()
132 : {
133 : // NOTE to readers: the computed properties depend on old and older quantities only. Hence, it
134 : // would be nice if we could compute them only at the start of each timestep instead of every
135 : // nonlinear iteration
136 :
137 : // For finite elements involved in boundary-condition calculations, number of nodes > number of
138 : // quadpoints. Material values on these elements are not computed during initial setup, which
139 : // means that for timestep = 1, _sat_older.size() = number of quadpoints. Since
140 : // computeQpProperties() is computed with 0 <= _qp < number of nodes (for nodal materials), this
141 : // results in an out-of-bounds error on _sat_older. PorousFlow's implementation of hysteresis
142 : // lags one timestep behind, and assumes that _sat_older = _sat_old at timestep = 1 anyway, so
143 : // just use it:
144 : const MaterialProperty<std::vector<Real>> & sat_older =
145 663882 : (_nodal_material && _sat_old.size() != _sat_older.size()) ? _sat_old : _sat_older;
146 :
147 : // whether saturation has been reducing ("drying" or "draining"):
148 663882 : const bool drying = (_hys_order_old[_qp] % 2 == 0);
149 : // whether have been drying but are now increasing saturation ("wetting" or "imbibing"):
150 663882 : const bool dry2wet = drying && (_sat_old[_qp][_liquid_ph_num] > sat_older[_qp][_liquid_ph_num]);
151 : // whether have been wetting but are now decreasing saturation:
152 663882 : const bool wet2dry = !drying && (_sat_old[_qp][_liquid_ph_num] < sat_older[_qp][_liquid_ph_num]);
153 663882 : if ((dry2wet || wet2dry) && _hys_order_old[_qp] < PorousFlowConstants::MAX_HYSTERESIS_ORDER)
154 : {
155 16452 : _hys_order[_qp] = _hys_order_old[_qp] + 1;
156 16452 : _hys_sat_tps[_qp] = _hys_sat_tps_old[_qp];
157 16452 : _hys_sat_tps[_qp].at(_hys_order_old[_qp]) = sat_older[_qp][_liquid_ph_num];
158 : }
159 : else
160 : {
161 : // no change in wetting/drying direction
162 647430 : _hys_order[_qp] = _hys_order_old[_qp];
163 647430 : _hys_sat_tps[_qp] = _hys_sat_tps_old[_qp];
164 : }
165 663882 : if (_sat_old[_qp][_liquid_ph_num] >= 1.0)
166 14588 : _hys_order[_qp] = 0; // assume the system gets back to the original drying curve
167 :
168 : // reduce order by 1 if _hys_order = max and saturation exceeds bounds of last turning-point
169 663882 : if (_hys_order[_qp] == PorousFlowConstants::MAX_HYSTERESIS_ORDER)
170 : {
171 : // the maxiumum-order curve is a drying curve, and the saturation exceeds the previous turning
172 : // point:
173 : const bool drying_and_high_sat =
174 : (_hys_order[_qp] % 2 == 0) &&
175 : (_sat_old[_qp][_liquid_ph_num] >= _hys_sat_tps[_qp].at(_hys_order[_qp] - 1));
176 : // the maxiumum-order curve is a wetting curve, and the saturation is less than the previous
177 : // turning point:
178 : const bool wetting_and_low_sat =
179 79828 : (_hys_order[_qp] % 2 == 1) &&
180 79828 : (_sat_old[_qp][_liquid_ph_num] <= _hys_sat_tps[_qp].at(_hys_order[_qp] - 1));
181 79828 : if (drying_and_high_sat || wetting_and_low_sat)
182 1536 : _hys_order[_qp] -= 1;
183 : }
184 :
185 : // reduce order by 2 if saturation exceeds the bounds of the second-to-last turning-point
186 : // encountered
187 663882 : bool can_reduce_order = (_hys_order[_qp] > 1);
188 669590 : while (can_reduce_order)
189 : {
190 : can_reduce_order = false;
191 : // are drying and the saturation is <= a previously-encountered turning point:
192 : const bool drying_and_low_sat =
193 205300 : (_hys_order[_qp] % 2 == 0) &&
194 127008 : (_sat_old[_qp][_liquid_ph_num] <= _hys_sat_tps[_qp].at(_hys_order[_qp] - 2));
195 : // are wetting and the saturation is >= a previously-encountered turning point:
196 : const bool wetting_and_high_sat =
197 205300 : (_hys_order[_qp] % 2 == 1) &&
198 78292 : (_sat_old[_qp][_liquid_ph_num] >= _hys_sat_tps[_qp].at(_hys_order[_qp] - 2));
199 205300 : if (drying_and_low_sat || wetting_and_high_sat)
200 : {
201 5708 : _hys_order[_qp] -= 2;
202 5708 : can_reduce_order = (_hys_order[_qp] > 1);
203 : }
204 : }
205 663882 : }
|
