- Generated by
1.9.8
|
Basilisk CFD
Adaptive Cartesian mesh PDE framework
|
Go to the source code of this file.
Macros | |
| #define | NH 1 |
Functions | |
| void | event_defaults (void) |
| Event: defaults (i = 0) | |
| void | event_viscous_term (void) |
| Event: viscous_term (i++) | |
| static void | box_matrix (Point point, scalar phi, scalar rhs, vector hf, scalar eta, double H[nl *nl], double d[nl]) |
| static trace void | relax_nh (scalar *phil, scalar *rhsl, int lev, void *data) |
| static trace double | residual_nh (scalar *phil, scalar *rhsl, scalar *resl, void *data) |
| void | event_pressure (void) |
| Event: pressure (i++) | |
| void | event_cleanup (void) |
| Event: cleanup (i = end, last) | |
Variables | |
| scalar | w |
| scalar | phi |
| The electric potential and the volume charge density are scalars while the permittivity and conductivity are face vectors. | |
| mgstats | mgp |
| double | breaking = HUGE |
| vector | hf |
| #define NH 1 |
This adds the non-hydrostatic terms of the vertically-Lagrangian multilayer solver for free-surface flows described in Popinet, 2020. The corresponding system of equations is
\[ \begin{aligned} \partial_t h_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \right)_k & = 0,\\ \partial_t \left( h \mathbf{u} \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \mathbf{u} \right)_k & = - gh_k \mathbf{{\nabla}} (\eta) {\color{blue} - \mathbf{{\nabla}} (h \phi)_k + \left[ \phi \mathbf{{\nabla}} z \right]_k},\\ {\color{blue} \partial_t (hw)_k + \mathbf{{\nabla}} \cdot \left( hw \mathbf{u} \right)_k} & {\color{blue} = - [\phi]_k,}\\ {\color{blue} \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \right)_k + \left[ w - \mathbf{u} \cdot \mathbf{{\nabla}} z \right]_k} & {\color{blue} = 0}, \end{aligned} \]
where the terms in blue are non-hydrostatic.
The additional \(w_k\) and \(\phi_k\) fields are defined. The convergence statistics of the multigrid solver are stored in mgp.
Wave breaking is parameterised usng the breaking parameter, which is turned off by default (see Section 3.6.4 in Popinet, 2020).
Note that this version differs in many ways from that presented in Popinet, 2020. The most important differences are the time-implicit integration of the barotropic free-surface evolution (explicit in the previous version) and the exact projection of the baroclinic velocities (approximate in the previous version). This results in the solution of a coupled system for the non-hydrostatic pressure \(\phi^{n+1}\) and the free-surface elevation \(\eta^{n+1}\).
See also the [general introduction](README).
|
static |
For the Keller box scheme, the linear system of equations verified by the non-hydrostatic pressure \(\phi\) is expressed as an Hessenberg matrix for each column.
The Hessenberg matrix \(\mathbf{H}\) for a column at a particular point is stored in a one-dimensional array with nl*nl elements. It encodes the coefficients of the left-hand-side of the Poisson equation as
\[ \begin{aligned} (\mathbf{H}\mathbf{\phi} - \mathbf{d})_l & = - \text{rhs}_l + h_l \nabla\cdot g_l^{n + \theta} +\\ & 4 (\phi_{l + 1 / 2} - \phi_{l - 1 / 2}) + 8 h_l \sum^{l - 1}_{k = 0} (- 1)^{l + k} \frac{\phi_{k + 1 / 2} - \phi_{k - 1 / 2}}{h_k}\\ g_l^{n + \theta} & = \nabla (h^{\star}_k \phi^{n + \theta}_{k - 1 / 2} + h^{\star}_k \phi^{n + \theta}_{k + 1 / 2}) - 2 [\phi \nabla \hat{z}]^{n + \theta}_k + 2 \theta gh^{\star}_k \nabla \eta^{n + 1} \end{aligned} \]
where \(\mathbf{\phi}\) is the vector of \(\phi_l\) for this column and \(\mathbf{d}\) is a vector dependent only on the values of \(\phi\) in the neighboring columns. Note that in contrast with Popinet, 2020, all the (metric) terms are retained.
Definition at line 116 of file nh.h.
References a, a_baro, cm, d, dimension, dry, eta, foreach_layer, gmetric, h, hf, k, l, m(), nl, phi, s, slope_limited, sq(), theta_H, vector::x, x, coord::x, and zb.
Referenced by relax_nh().
Event: defaults (i = 0)
The \(w_k\) and \(\phi_k\) scalar fields are allocated and the \(w_k\) are added to the list of advected tracers.
The default density field is set to unity (times the metric).
We reset the multigrid parameters to their default values.
If the viscosity is non-zero, we need to allocate the face-centered viscosity field.
We also initialize the list of tracers to be advected with the VOF function \(f\) (or its complementary function).
We set limiting.
On trees, we ensure that limiting is also applied to prolongation and refinement.
We add the interface and the density to the default display.
We switch to a pure minmod limiter by default for increased robustness.
With the MUSCL scheme we use the CFL depends on the dimension of the problem.
On trees we need to replace the default bilinear refinement/prolongation with linear so that reconstructed values also use slope limiting.
The restriction/refine attributes of the charge density are those of a tracer otherwise the conservation is not guaranteed.
By default the permittivity is unity and other quantities are zero.
The (velocity) CFL is limited by the unsplit advection scheme, so is dependent on the dimension. The (gravity wave) CFL is set to 1/2 (if not already set by the user).
The gradient and prolongation/restriction functions are set for all tracer fields.
We setup the default display.
Definition at line 59 of file nh.h.
References assert, hydrostatic, linearised, list_append(), mgp, nl, mgstats::nrelax, phi, reset, tracers, and w.
Event: pressure (i++)
The coupled system for \(\phi_k^{n+1}\) and \(\eta^{n+1}\) is solved using the multigrid solver.
The r.h.s. is computed as
\[ \frac{2 h_k}{\theta \Delta t} \left( 2 w^n_k + \nabla \cdot (hu)_k^{\star} - [u \cdot \nabla \hat{z}]^\star_k + 4 \sum^{k - 1}_{l = 0} (- 1)^{k + l} w^n_l \right) \]
Note that the discrete approximation below must verify Galilean invariance i.e.
\[ \begin{aligned} \nabla \cdot (h(u + U_0))_k^{\star} - [(u + U_0)\cdot \nabla \hat{z}]^\star_k & = \nabla \cdot (hu)_k^{\star} - [u \cdot \nabla \hat{z}]^\star_k \end{aligned} \]
with \(U_0\) an arbitrary constant vector. This can be simplified as
\[ \nabla \cdot (h_k^{\star}U_0) - [U_0\cdot \nabla \hat{z}]^\star_k = 0 \]
or further
\[ \nabla h_k^\star = [\nabla \hat{z}]^\star_k \]
which is obviously verified (analytically) since by definition
\[ h_k^\star = [\hat{z}]^\star_k \]
Note that it is not so obvious that this is verified numerically, as this depends on the choices made for several approximations. In particular, in the expressions below, Galilean invariance implies the relations
which depend on the detail of the calculation of hu. Note also that the slope limiter will break Galilean invariance.
The fields used by the relaxation function above need to be restricted to all levels. Note that the other fields (cm, fm, alpha_eta) were already restricted in the hydrostatic implicit solver.
We then call the multigrid solver, using the relaxation and residual functions defined above, to get both the non-hydrostatic pressure \(\phi\) and free-surface elevation \(\eta\).
The non-hydrostatic pressure gradient is added to the face-weighted acceleration and to the face fluxes.
The maximum allowed vertical velocity is computed as
\[ w_\text{max} = b \sqrt{g | H |_{\infty}} \]
with \(b\) the breaking parameter.
The vertical pressure gradient is added to the vertical velocity as
\[ w^{n + 1}_l = w^{\star}_l - \Delta t \frac{[\phi]_l}{h^{n+1}_l} \]
and the vertical velocity is limited by \(w_\text{max}\) as
\[ w^{n + 1}_l \leftarrow \text{sign} (w^{n + 1}_l) \min \left( | w^{n + 1}_l |, w_\text{max} \right) \]
The r.h.s. for \(\eta^{n+1}\) is updated. It will be used in a second pass (neglecting the non-hydrostatic terms) in the semi-implicit free-surface solver. This should be a small correction which is only necessary to limit the accumulation of divergence noise for long integration times.
Definition at line 313 of file nh.h.
References _i, alpha_eta, breaking, cm, dimension, dry, dt, dv, eta_r, fm, foreach_layer, G, h, ha, hf, hpg, hu, HUGE, scalar::i, k, mg_solve(), mgp, nl, phi, point, relax_nh(), res_eta, residual_nh(), restriction, rhs_eta, s, slope_limited, sq(), theta_H, TOLERANCE, u, w, wmax, vector::x, x, coord::x, and zb.
The updated values of \(\phi\) in a column are obtained as
\[ \mathbf{\phi} = \mathbf{H}^{-1}\mathbf{b} \]
were \(\mathbf{H}\) and \(\mathbf{b}\) are the Hessenberg matrix and vector constructed by the function above.
The value of \(\eta\) also needs to be updated since it is solved implicitly and depends on \(\phi\).
Definition at line 173 of file nh.h.
References _i, a_baro, alpha, b, box_matrix(), cm, d, data, diagonalize(), dimension, dt, eta, foreach_layer, hf, hpg, Point::i, l, level, n, nl, phi, phil, point, rhs_eta, rigid, solve_hessenberg(), sq(), theta_H, vector::x, and x.
Referenced by event_pressure().
The residual for \(\phi\) is computed as
\[ \begin{aligned} \text{res}_l = & \text{rhs}_l - h_l \nabla\cdot g_l^{n + \theta} - 4 (\phi_{l + 1 / 2} - \phi_{l - 1 / 2}) - 8 h_l \sum^{l - 1}_{k = 0} (- 1)^{l + k} \frac{\phi_{k + 1 / 2} - \phi_{k - 1 / 2}}{h_k} \end{aligned} \]
The residual for \(\eta\) is computed as
\[ \text{res}_\eta = \text{rhs}_\eta - \beta\eta + \theta_H \Delta t^2 \sum_l \nabla \cdot (- \nabla_z\phi_l + \theta_H h_l g \nabla \eta) \]
where \(\beta = 1\) for a free surface and \(\beta = 0\) for a rigid lid.
Definition at line 233 of file nh.h.
References _i, a_baro, cm, dimension, dry, dt, eta, fm, foreach_layer, g, h, hf, hpg, k, max, nl, phi, phil, point, res_eta, rhs_eta, rigid, s, slope_limited, sq(), theta_H, vector::x, x, coord::x, and zb.
Referenced by event_pressure().
Definition at line 50 of file nh.h.
Referenced by event_pressure().
| vector hf |
Definition at line 170 of file nh.h.
Referenced by box_matrix(), event_pressure(), relax_nh(), and residual_nh().
| mgstats mgp |
Definition at line 49 of file nh.h.
Referenced by event_defaults(), and event_pressure().
| scalar phi |
The electric potential and the volume charge density are scalars while the permittivity and conductivity are face vectors.
This function computes the fluxes due to ohmic conduction appearing in the Nernst–Planck equation. The species charge concentrations are then updated using the explicit scheme
\[ c^{n+1}_i = c^n_i +\Delta t \, \nabla \cdot( K_i c^n_i \nabla \phi^n) \]
where \(c_i\) is the volume density of the \(i\)-specie, \(K_i\) its volume electric conductivity and \(\phi\) the electric potential.
In mgphi we will store the statistics for the multigrid resolution of the electric poisson equation.
Definition at line 48 of file nh.h.
Referenced by box_matrix(), event_cleanup(), event_defaults(), event_pressure(), ohmic_flux(), relax_nh(), and residual_nh().
| scalar w |
Definition at line 48 of file nh.h.
Referenced by event_cleanup(), event_defaults(), event_pressure(), and event_viscous_term().
1.9.8