RichardsMechanics; Update mass lumping procedure. #2921
Conversation
endJunction
commented
Apr 26, 2020
endJunction
commented
- Feature description was added to the changelog
- Tests covering your feature were added? updated.
- Any new feature or behavior change was documented?
This allows for mass lumping. Part of the pressure-pressure jacobian block is now also diagnalized if mass lumping is used.
|
looks good |
It replaces S' with S/Delta p N p'_nodes, in order to turn a nodal vector into the product of a "mass" matrix (still don't like this term) and a nodal vector. This way, we can apply mass lumping. In contrast to a \partial S/\partial p formulation we maintain the good mass conservation properties of the S' formulation. |
| storage_p_a_S.noalias() += N_p.transpose() * rho_LR * | ||
| specific_storage_a_S * (S_L - S_L_prev) / | ||
| dt * w; | ||
| if (p_cap_dot_ip != 0) // prevent division by zero. |
There was a problem hiding this comment.
If p_cap_dot_ip is very small, this condition would cause a problem.
There was a problem hiding this comment.
If
p_cap_dot_ipis very small, this condition would cause a problem.
Yeah. In the saturated range p_cap_dot_ip = -p_L_dot_ip, and the potentially small pressure rate goes along with S_L = S_L_prev. So this case should be fine, I guess.
In the unsaturated case, the ratio of both depends on the shape of the SWCC, but in many cases it should work.
@endJunction if none of the benchmarks nor our additional tests run into any problems of this kind, then I'd be fine with this for now and only address it if an actual problem occurs.
There was a problem hiding this comment.
So far I didn't come across a case where it would produce anything strange.
Thanks for the explanation. |
|
OpenGeoSys development has been moved to GitLab. |
