OGS 6.2.0-97-g4a610c866
TESReactionAdaptor.cpp
Go to the documentation of this file.
1 
10 #include <cassert>
11 
12 #include <logog/include/logog.hpp>
13 
16 
20 
21 #include "TESLocalAssemblerInner.h"
22 #include "TESReactionAdaptor.h"
23 
24 namespace ProcessLib
25 {
26 namespace TES
27 {
28 std::unique_ptr<TESFEMReactionAdaptor> TESFEMReactionAdaptor::newInstance(
29  TESLocalAssemblerData const& data)
30 {
31  auto const* ads = data.ap.react_sys.get();
32  if (dynamic_cast<Adsorption::AdsorptionReaction const*>(ads) != nullptr)
33  {
34  return std::make_unique<TESFEMReactionAdaptorAdsorption>(data);
35  }
36  if (dynamic_cast<Adsorption::ReactionInert const*>(ads) != nullptr)
37  {
38  return std::make_unique<TESFEMReactionAdaptorInert>(data);
39  }
40  if (dynamic_cast<Adsorption::ReactionSinusoidal const*>(ads) != nullptr)
41  {
42  return std::make_unique<TESFEMReactionAdaptorSinusoidal>(data);
43  }
44  if (dynamic_cast<Adsorption::ReactionCaOH2 const*>(ads) != nullptr)
45  {
46  return std::make_unique<TESFEMReactionAdaptorCaOH2>(data);
47  }
48 
49  OGS_FATAL("No suitable TESFEMReactionAdaptor found. Aborting.");
50  return nullptr;
51 }
52 
54  TESLocalAssemblerData const& data)
55  // caution fragile: this relies in this constructor b eing called __after__
56  // data.solid_density has been properly set up!
57  : _bounds_violation(data.solid_density.size(), false),
58  _d(data)
59 {
60  assert(dynamic_cast<Adsorption::AdsorptionReaction const*>(
61  data.ap.react_sys.get()) != nullptr &&
62  "Reactive system has wrong type.");
63  assert(!_bounds_violation.empty());
64 }
65 
68  const unsigned int_pt)
69 {
70  assert(_d.ap.number_of_try_of_iteration <= 20);
71 
72  const double loading = Adsorption::AdsorptionReaction::getLoading(
74 
75  double react_rate_R =
76  _d.ap.react_sys->getReactionRate(_d.p_V, _d.T, _d.ap.M_react, loading) *
78 
79  // set reaction rate based on current damping factor
80  react_rate_R = (_reaction_damping_factor > 1e-3)
81  ? _reaction_damping_factor * react_rate_R
82  : 0.0;
83 
84  if (_d.p_V <
86  _d.T) &&
87  react_rate_R > 0.0)
88  {
89  react_rate_R = 0.0;
90  }
91  else if (_d.p_V < 100.0 ||
93  getEquilibriumVapourPressure(_d.T))
94  {
95  // use equilibrium reaction for dry regime
96 
97  // in the case of zeroth try in zeroth iteration: _p_V and loading are
98  // the values
99  // at the end of the previous timestep
100 
101  const double pV_eq = estimateAdsorptionEquilibrium(_d.p_V, loading);
102  // TODO [CL]: it would be more correct to subtract pV from the previous
103  // timestep here
104  const double delta_pV = pV_eq - _d.p_V;
105  const double delta_rhoV =
106  delta_pV * _d.ap.M_react /
108  _d.ap.poro;
109  const double delta_rhoSR = delta_rhoV / (_d.ap.poro - 1.0);
110  double react_rate_R2 = delta_rhoSR / _d.ap.delta_t;
111 
112  if (_bounds_violation[int_pt])
113  {
114  react_rate_R2 *= 0.5;
115  }
116 
117  // 0th try: make sure reaction is not slower than allowed by local
118  // estimation
119  // nth try: make sure reaction is not made faster by local estimation
120  if ((_d.ap.number_of_try_of_iteration == 1 &&
121  std::abs(react_rate_R2) > std::abs(react_rate_R)) ||
123  std::abs(react_rate_R2) < std::abs(react_rate_R)))
124  {
125  react_rate_R = react_rate_R2;
126  }
127  }
128 
129  // smooth out readjustment of reaction rate
131  {
133  {
134  // update reaction rate for for five iterations
135  const auto N = _d.ap.iteration_in_current_timestep - 4;
136 
137  // take average s.t. does not oscillate so much
138  react_rate_R = 1.0 / (1.0 + N) *
139  (N * _d.reaction_rate[int_pt] + 1.0 * react_rate_R);
140  }
141  else
142  {
143  // afterwards no update anymore
144  react_rate_R = _d.reaction_rate[int_pt];
145  }
146  }
147 
149  {
150  // assert that within tries reaction does not get faster
151  // (e.g. due to switch equilibrium reaction <--> kinetic reaction)
152 
153  // factor of 0.9*N: in fact, even slow down reaction over tries
154  const double r = std::pow(0.9, _d.ap.number_of_try_of_iteration - 1) *
155  _d.reaction_rate[int_pt];
156  if (std::abs(react_rate_R) > std::abs(r))
157  {
158  react_rate_R = r;
159  }
160  }
161 
162  return {react_rate_R,
163  _d.solid_density_prev_ts[int_pt] + react_rate_R * _d.ap.delta_t};
164 }
165 
167  const double p_V0, const double C0) const
168 {
169  auto f = [this, p_V0, C0](double pV) -> double {
170  // pV0 := _p_V
171  const double C_eq =
172  _d.ap.react_sys->getEquilibriumLoading(pV, _d.T, _d.ap.M_react);
173  return (pV - p_V0) * _d.ap.M_react /
175  _d.ap.poro +
176  (1.0 - _d.ap.poro) * (C_eq - C0) * _d.ap.rho_SR_dry;
177  };
178 
179  // range where to search for roots of f
180  const double C_eq0 =
181  _d.ap.react_sys->getEquilibriumLoading(p_V0, _d.T, _d.ap.M_react);
182  const double limit =
183  (C_eq0 > C0)
184  ? 1e-8
186  _d.T);
187 
188  // search for roots
189  auto rf = MathLib::Nonlinear::makeRegulaFalsi<MathLib::Nonlinear::Pegasus>(
190  f, p_V0, limit);
191  rf.step(3);
192 
193  // set vapour pressure
194  return rf.getResult();
195 }
196 
198  std::vector<double> const& local_x,
199  std::vector<double> const& local_x_prev_ts)
200 {
201  double alpha = 1.0;
202 
203  const double min_xmV = 1e-6;
204  const std::size_t nnodes = local_x.size() / NODAL_DOF;
205  const std::size_t xmV_offset = COMPONENT_ID_MASS_FRACTION * nnodes;
206 
207  for (std::size_t i = 0; i < nnodes; ++i)
208  {
209  auto const xnew = local_x[xmV_offset + i];
210  auto const xold = local_x_prev_ts[xmV_offset + i];
211 
212  if (xnew < min_xmV)
213  {
214  const auto a = xold / (xold - xnew);
215  alpha = std::min(alpha, a);
216  _bounds_violation[i] = true;
217  }
218  else if (xnew > 1.0)
219  {
220  const auto a = xold / (xnew - xold);
221  alpha = std::min(alpha, a);
222  _bounds_violation[i] = true;
223  }
224  else
225  {
226  _bounds_violation[i] = false;
227  }
228  }
229 
230  assert(alpha > 0.0);
231 
232  if (alpha != 1.0)
233  {
234  if (alpha > 0.5)
235  {
236  alpha = 0.5;
237  }
238  if (alpha < 0.05)
239  {
240  alpha = 0.05;
241  }
242 
244  {
245  _reaction_damping_factor *= sqrt(alpha);
246  }
247  else
248  {
249  _reaction_damping_factor *= alpha;
250  }
251  }
252 
253  return alpha == 1.0;
254 }
255 
257 {
258  if (_reaction_damping_factor < 1e-3)
259  {
261  }
262 
264  10.0 * _reaction_damping_factor);
265 }
266 
268  TESLocalAssemblerData const& data)
269  : _d(data)
270 {
271 }
272 
274 {
275  return {0.0, _d.solid_density_prev_ts[int_pt]};
276  // _d._qR = 0.0;
277  // _d._reaction_rate[int_pt] = 0.0;
278 }
279 
281  TESLocalAssemblerData const& data)
282  : _d(data)
283 {
284  assert(dynamic_cast<Adsorption::ReactionSinusoidal const*>(
285  data.ap.react_sys.get()) != nullptr &&
286  "Reactive system has wrong type.");
287 }
288 
290  const unsigned /*int_pt*/)
291 {
292  const double t = _d.ap.current_time;
293 
294  // Cf. OGS5
295  const double rhoSR0 = 1.0;
296  const double rhoTil = 0.1;
297  const double omega = 2.0 * 3.1416;
298  const double poro = _d.ap.poro;
299 
300  return {rhoTil * omega * cos(omega * t) / (1.0 - poro),
301  rhoSR0 + rhoTil * std::sin(omega * t) / (1.0 - poro)};
302 }
303 
305  TESLocalAssemblerData const& data)
306  : _d(data),
307  _react(dynamic_cast<Adsorption::ReactionCaOH2&>(*data.ap.react_sys))
308 {
309  _ode_solver = MathLib::ODE::createODESolver<1>(_react.getOdeSolverConfig());
310  // TODO invalidate config
311 
312  _ode_solver->setTolerance(1e-10, 1e-10);
313 
314  auto f = [this](const double /*t*/,
317  ydot) -> bool {
318  ydot[0] = _react.getReactionRate(y[0]);
319  return true;
320  };
321 
322  _ode_solver->setFunction(f, nullptr);
323 }
324 
326 {
327  // TODO if the first holds, the second also has to hold
330  {
331  return {_d.reaction_rate[int_pt], _d.solid_density[int_pt]};
332  }
333 
334  // TODO: double check!
335  // const double xv_NR = SolidProp->non_reactive_solid_volume_fraction;
336  // const double rho_NR = SolidProp->non_reactive_solid_density;
337  const double xv_NR = 0.0;
338  const double rho_NR = 0.0;
339 
340  const double t0 = 0.0;
341  const double y0 =
342  (_d.solid_density_prev_ts[int_pt] - xv_NR * rho_NR) / (1.0 - xv_NR);
343 
344  const double t_end = _d.ap.delta_t;
345 
347  _d.solid_density_prev_ts[int_pt]);
348 
349  _ode_solver->setIC(t0, {y0});
350  _ode_solver->preSolve();
351  _ode_solver->solve(t_end);
352 
353  const double time_reached = _ode_solver->getTime();
354  (void)time_reached;
355  assert(std::abs(t_end - time_reached) <
356  std::numeric_limits<double>::epsilon());
357 
358  auto const& y_new = _ode_solver->getSolution();
359  auto const& y_dot_new = _ode_solver->getYDot(t_end, y_new);
360 
361  double rho_react;
362 
363  // cut off when limits are reached
364  if (y_new[0] < _react.rho_low)
365  {
366  rho_react = _react.rho_low;
367  }
368  else if (y_new[0] > _react.rho_up)
369  {
370  rho_react = _react.rho_up;
371  }
372  else
373  {
374  rho_react = y_new[0];
375  }
376 
377  return {y_dot_new[0] * (1.0 - xv_NR),
378  (1.0 - xv_NR) * rho_react + xv_NR * rho_NR};
379 }
380 
381 } // namespace TES
382 } // namespace ProcessLib
TESFEMReactionAdaptorCaOH2(TESLocalAssemblerData const &data)
static const double r
MappedConstMatrix< N, 1 > MappedConstVector
void updateParam(double T_solid, double _p_gas, double _x_react, double rho_s_initial)
static double getEquilibriumVapourPressure(const double T_Ads)
Definition: Adsorption.cpp:28
double estimateAdsorptionEquilibrium(const double p_V0, const double C0) const
static MATERIALLIB_EXPORT const double rho_up
lower density limit
Definition: ReactionCaOH2.h:84
ReactionRate initReaction(const unsigned) override
bool checkBounds(std::vector< double > const &local_x, std::vector< double > const &local_x_prev_ts) override
TESFEMReactionAdaptorSinusoidal(TESLocalAssemblerData const &data)
ReactionRate initReaction_slowDownUndershootStrategy(const unsigned int_pt)
TESFEMReactionAdaptorAdsorption(TESLocalAssemblerData const &data)
const BaseLib::ConfigTree & getOdeSolverConfig() const
Definition: ReactionCaOH2.h:41
ReactionRate initReaction(const unsigned) override
static double getLoading(const double rho_curr, const double rho_dry)
Definition: Adsorption.cpp:152
const unsigned COMPONENT_ID_MASS_FRACTION
TESFEMReactionAdaptorInert(TESLocalAssemblerData const &)
std::unique_ptr< MathLib::ODE::ODESolver< 1 > > _ode_solver
MappedMatrix< N, 1 > MappedVector
std::unique_ptr< Adsorption::Reaction > react_sys
const unsigned NODAL_DOF
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
double getReactionRate(const double, const double, const double, const double) const override
static std::unique_ptr< TESFEMReactionAdaptor > newInstance(TESLocalAssemblerData const &data)
ReactionRate initReaction(const unsigned int_pt) override
static MATERIALLIB_EXPORT const double rho_low
Definition: ReactionCaOH2.h:83