- Generated by
1.9.8
|
Basilisk CFD
Adaptive Cartesian mesh PDE framework
|
Go to the source code of this file.
Macros | |
| #define | dx(s) ((s[1,0] - s[-1,0])/(2.*Delta)) |
| We first define some useful macros, following the notations in Bonneton et al, 2011. | |
| #define | dy(s) ((s[0,1] - s[0,-1])/(2.*Delta)) |
| #define | d2x(s) ((s[1,0] + s[-1,0] - 2.*s[])/sq(Delta)) |
| #define | d2y(s) ((s[0,1] + s[0,-1] - 2.*s[])/sq(Delta)) |
| #define | d2xy(s) ((s[1,1] - s[1,-1] - s[-1,1] + s[-1,-1])/sq(2.*Delta)) |
| #define | R1(h, zb, w) (-h[]*(h[]/3.*dx(w) + w[]*(dx(h) + dx(zb)/2.))) |
| The definitions of the \(\mathcal{R}_1\) and \(\mathcal{R}_2\) operators. | |
| #define | R2(h, zb, w) (h[]/2.*dx(w) + w[]*(dx(zb) + dx(h))) |
Functions | |
| static double | update_green_naghdi (scalar *current, scalar *updates, double dtmax) |
| void | event_defaults (void) update |
| Event: defaults (i = 0) | |
| static double | residual_GN (scalar *a, scalar *r, scalar *resl, void *data) |
| static void | relax_GN (scalar *a, scalar *r, int l, void *data) |
| The relaxation function is built by copying and pasting the residual implementation above and inverting manually for \(D_x\). | |
Variables | |
| vector | D [] |
| The linear system can be inverted with the multigrid Poisson solver. | |
| mgstats | mgD |
| double | breaking = 1. |
| double | alpha_d = 1.153 |
Definition at line 84 of file green-naghdi.h.
Definition at line 86 of file green-naghdi.h.
Definition at line 85 of file green-naghdi.h.
We first define some useful macros, following the notations in Bonneton et al, 2011.
Simple centered-difference approximations of the first and second derivatives of a field.
Definition at line 82 of file green-naghdi.h.
The definitions of the \(\mathcal{R}_1\) and \(\mathcal{R}_2\) operators.
\[ \begin{array}{lll} \mathcal{R}_1 \left[ h, z_b \right] w & = & - \frac{1}{3 h} \nabla \left( h^3 w \right) - \frac{h}{2} w \nabla z_b\\ & = & - h \left[ \frac{h^{}}{3} \nabla w + w \left( \nabla h + \frac{1}{2} \nabla z_b \right)\right]\\ \mathcal{R}_2 \left[ h, z_b \right] w & = & \frac{1}{2 h} \nabla \left( h^2 w \right) + w \nabla z_b\\ & = & \frac{h}{2} \nabla w + w \nabla \left( z_b + h \right) \end{array} \]
Definition at line 102 of file green-naghdi.h.
Definition at line 103 of file green-naghdi.h.
Event: defaults (i = 0)
Boundary conditions for VOF-advected tracers usually depend on boundary conditions for the VOF field.
Event: defaults (i = 0)
Event: defaults (i = 0)
We set the default values.
The EHD force density, \(\mathbf{f}_e\), can be computed as the divergence of the Maxwell stress tensor \(\mathbf{M}\),
\[ M_{ij} = \varepsilon (E_i E_j - \frac{E^2}{2}\delta_{ij}) \]
where \(E_i\) is the \(i\)-component of the electric field, \(\mathbf{E}=-\nabla \phi\) and \(\delta_{ij}\) is the Kronecker delta.
We need to add the corresponding acceleration to the Navier–Stokes solver.
If the acceleration vector a (defined by the Navier–Stokes solver) is constant, we make it variable.
Event: defaults (i = 0)
Event: defaults (i = 0 )
On trees we need to ensure conservation of the tracer when refining/coarsening.
Event: defaults (i = 0)
it is an acceleration. If necessary, we allocate a new vector field to store it.
Event: defaults (i = 0)
Event: defaults (i = 0)
\[ \begin{aligned} 0 & = - \sum_k \nabla \cdot [\theta_H (hu)_k^{n + 1} + (1 - \theta_H) (hu)^n_k] \\ \frac{(hu)^{n + 1}_k - (hu)_k^n}{\Delta t} & = - \Delta tgh^{n + 1 / 2}_k (\theta_H \nabla \eta_r^{n + 1} + (1 - \theta_H) \nabla \eta_r^n) \end{aligned} \]
where \(\eta_r\) is the equivalent free-surface height (i.e. pressure) applied on the rigid lid.
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.
Event: defaults (i = 0)
Event: defaults (i = 0)
Event: defaults (i = 0)
The default convention in Basilisk is no-flow through the boundaries of the domain, i.e. they are a streamline i.e. \(\psi=\)constant on the boundary. We set the default value for the CFL (the default in utils.h is 0.5). This is done once at the beginning of the simulation.
Event: defaults (i = 0)
Event: defaults (i = 0)
We will need to add the acceleration term \(w^2/y\) in the evolution equation for \(u_y\). If the acceleration field is not allocated yet, we do so.
Event: defaults (i = 0)
We use the main time loop (in the predictor-corrector scheme) to setup the initial defaults.
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.
By default we set a zero Neumann boundary condition for all the components except if the bottom is an axis of symmetry.
We use (strict) minmod slope limiting for all components.
We reset the multigrid parameters to their default values.
The pressures are never dumped.
The default density field is set to unity (times the metric and the solid factors).
On trees, refinement of the face-centered velocity field needs to preserve the divergence-free condition.
When using embedded boundaries, the restriction and prolongation operators need to take the boundary into account.
We set the dimensions of the velocity field.
On trees, the refinement and restriction functions above rely on the volume fraction field f being refined/restricted before the components of velocity. To ensure this, we move f to the front of the field list (all).
We then set the refinement and restriction functions for the components of the velocity field. The boundary conditions on \(\mathbf{u}\) now depend on those on \(f\).
The default display.
We set limiting for q, h.
As well as default values.
On trees, we ensure that limiting is also applied to prolongation.
We setup the default display.
By default, all the layers have the same relative thickness.
We overload the default 'advance' and 'update' functions of the predictor-corrector scheme and setup the prolongation and restriction methods on trees.
On trees we make sure that slope-limiting is also used for refinement and prolongation. The prolongation/restriction functions for \(\eta\) are set and they depend on boundary conditions on \(z_b\) and \(h\).
We setup the default display.
If the viscosity is non-zero, we need to allocate the face-centered viscosity field.
We add the interface to the default display.
Definition at line 40 of file advection.h.
References _i, a, advance, advance_saint_venant(), all, alpha, alphav, assert, av, beta, c, calloc(), CFL, CFL_H, cm, cs, dimension, display(), epsilon, eta, eta_r, evolving, f, f_r, f_s, fclone(), fE1, fE2, fenep(), fm, fraction_refine(), frho1, frho2, G, gammap, gotm_buoy_method, gotm_cnpar, gradient, h, hydrostatic, scalar::i, i, interfaces, K, kappa, kappa1, kappa2, kpp_init_kpp(), l, lambdav, layer, left, linearised, list_add(), list_append(), list_concat(), list_copy(), list_len(), meanflow_depth, meanflow_gravity, meanflow_h, meanflow_post_init_meanflow(), meanflow_rho_0, mgp, mgpf, mgu, minmod2(), momentum_refine(), momentum_restriction(), mtridiagonal_init_tridiagonal(), mu, mu1, mu2, muv, nl, mgstats::nrelax, number_of_interfaces, p, periodic_bc(), pf, phi, phil, q, qmalloc, qzl, refine_embed_linear(), refine_eta(), refine_face(), refine_face_solenoidal(), refine_linear(), reset, restriction_embed_linear(), restriction_eta(), restriction_volume_average(), rho, rhoc2, rhoc2v, rhoe, rhov, rigid, s, scalars, set_prolongation(), set_restriction(), stokes, swap, t, T, tau_p, theta, trA, tracers, turbulence_post_init_turbulence(), turbulence_turb_method, u, uf, ul, unity, unityf, update, update_conservation(), update_saint_venant(), vectors, vectors_append(), vof_concentration_refine(), w, wl, vector::x, tensor::x, x, vector::y, and zb.
The relaxation function is built by copying and pasting the residual implementation above and inverting manually for \(D_x\).
Definition at line 200 of file green-naghdi.h.
References _i, a, alpha_d, b, cube(), D, d2x, d2xy, data, dimension, dry, dx, dy, h, Point::i, l, level, list, point, sq(), vector, vector::x, x, vector::y, and zb.
Referenced by update_green_naghdi().
To invert the linear system which defines \(D\), we need to write discretised versions of the residual and relaxation operators. The functions take generic multilevel parameters and a user-defined structure which contains solution-specific parameters, in our case a list of the \(h\), \(zb\) and wet fields.
We first recover all the parameters from the generic pointers and rename them according to our notations.
The general form for \(\mathcal{T}\) is
\[ \mathcal{T} \left[ h, z_b \right] w = \mathcal{R}_1 \left[ h, z_b \right] \left( \nabla \cdot w \right) + \mathcal{R}_2 \left[ h, z_b \right] \left( \nabla z_b \cdot w \right) \]
which gives the linear problem
\[ \alpha_d h \mathcal{T} \left( D \right) + hD = b \]
\[ \alpha_d h \mathcal{R}_1 \left( \nabla \cdot D \right) + \alpha_d h \mathcal{R}_2 \left( \nabla z_b \cdot D \right) + hD = b \]
\[ - \alpha_d h^2 \left[ \frac{h}{3} \nabla \left( \nabla \cdot D \right) + \left( \nabla \cdot D \right) \left( \nabla h + \frac{1}{2} \nabla z_b \right)\right] + \alpha_d h \left[ \frac{h}{2} \nabla \left( \nabla z_b \cdot D \right) + \left( \nabla z_b \cdot D \right) \nabla \eta \right] + hD = b \]
Expanding the operators for the \(x\)-component gives
\[ \begin{aligned} - \frac{\alpha_d}{3} \partial_x \left( h^3 \partial_x D_x \right) + h \left[ \alpha_d \left( \partial_x \eta \partial_x z_b + \frac{h}{2} \partial^2_x z_b \right) + 1 \right] D_x + & \\ \alpha_d h \left[ \left( \frac{h}{2} \partial^2_{xy} z_b + \partial_x \eta \partial_y z_b \right) D_y + \frac{h}{2} \partial_y z_b \partial_x D_y - \frac{h^2}{3} \partial^2_{xy} D_y - h \partial_y D_y \left( \partial_x h + \frac{1}{2} \partial_x z_b \right) \right] & = b_x \end{aligned} \]
The \(y\)-component is obtained by rotation of the indices.
Finally we translate the formula above to its discrete version.
The function also need to return the maximum residual.
The maximum residual is normalised by gravity i.e. the tolerance is the relative acceleration times the depth.
Definition at line 114 of file green-naghdi.h.
References a, alpha_d, b, cube(), D, d2x, d2xy, data, dimension, dx, dy, G, h, list, max, sq(), vector, vector::x, x, vector::y, and zb.
Referenced by update_green_naghdi().
The Green-Naghdi equations (also known as the Serre or fully nonlinear Boussinesq equations) can be seen as an extension of the Saint-Venant equations to the next order \(O(\mu^2)\) in term of the shallowness parameter*
\[ \mu = \frac{h_0^2}{L^2} \]
with \(h_0\) the typical depth and \(L\) the typical horizontal scale. In contrast to the Saint-Venant equations the Green-Naghdi equations have dispersive* wave solutions.
A more detailed description of the context and numerical scheme implemented here is given in Popinet, 2015.
The solver is built by adding a source term to the momentum equation of the Saint-Venant solver. Following Bonneton et al, 2011, this source term can be written
\[ \partial_t \left( hu \right) + \ldots = h \left( \frac{g}{\alpha_d} \nabla \eta - D \right) \]
where \(D\) verifies
\[ \alpha_d h\mathcal{T} \left( D \right) + hD = b \]
and
\[ b = \left[ \frac{g}{\alpha_d} \nabla \eta +\mathcal{Q}_1 \left( u \right) \right] \]
With \(\mathcal{T}\) a linear operator to be defined below, as well as \(\mathcal{Q}_1 \left( u \right)\).
Before including the Saint-Venant solver, we need to overload the default update function of the predictor-corrector scheme in order to add our source term.
To add the source term to the Saint-Venant system we overload the default update function with this one. The function takes a list of the current evolving scalar fields and fills the corresponding updates with the source terms.
We first compute the right-hand-side \(b\). The general form for the \(\mathcal{Q}_1\) operator is (eq. (9) of Bonneton et al, 2011).
\[ \mathcal{Q}_1 \left[ h, z_b \right] \left( u \right) = - 2 \mathcal{R}_1 \left( c \right) + \mathcal{R}_2 \left( d \right) \]
with
\[ \begin{aligned} c & = \partial_1 u \cdot \partial_2 u^{\perp} + \left( \nabla \cdot u \right)^2\\ & = - \partial_x u_x \partial_y u_y + \partial_x u_y \partial_y u_x + \left( \partial_x u_x + \partial_y u_y \right)^2\\ d & = u \cdot \left( u \cdot \nabla \right) \nabla z_b\\ & = u_x \left( u_x \partial_x^2 z_b + u_y \partial_{xy}^2 z_b \right) + u_y \left( u_y \partial_y^2 z_b + u_x \partial_{xy}^2 z_b \right)\\ & = u^2_x \partial_x^2 z_b + u^2_y \partial_y^2 z_b + 2 u_x u_y \partial_{xy}^2 z_b \end{aligned} \]
Note that we declare field c and d in a new scope, so that the corresponding memory is freed after we have computed \(b\).
We can now compute \(b\)
\[ b = \left[ \frac{g}{\alpha_d} \nabla \eta +\mathcal{Q}_1 \left( u \right) \right] \]
We declare a new field which will track whether cells are completely wet. This is necessary to turn off dispersive terms close to the shore so that lake-at-rest balance is maintained.
Finally we solve the linear problem with the multigrid solver using the relaxation and residual functions defined above. We need to restrict \(h\) and \(z_b\) as their values will be required for relaxation on coarse grids.
We can then compute the updates for \(hu\).
We only apply the Green-Naghdi source term when the slope of the free surface is smaller than breaking. The idea is to turn off dispersion in areas where the wave is "breaking" (i.e. in hydraulic jumps). We also turn off dispersive terms close to shore so that lake-at-rest balance is maintained.
Definition at line 239 of file green-naghdi.h.
References _i, alpha_d, b, breaking, c, d, D, d2x, d2xy, d2y, dimension, dry, dt, dtmax, dx, dy, eta, G, h, list, mg_solve(), mgD, mgstats::nrelax, R1, R2, relax_GN(), residual_GN(), restriction, sq(), u, update_saint_venant(), vector, vector::x, x, and zb.
| double alpha_d = 1.153 |
Definition at line 73 of file green-naghdi.h.
Referenced by relax_GN(), residual_GN(), and update_green_naghdi().
| double breaking = 1. |
Definition at line 73 of file green-naghdi.h.
Referenced by update_green_naghdi().
| vector D[] |
The linear system can be inverted with the multigrid Poisson solver.
We declare D as a global variable so that it can be re-used as initial guess for the Poisson solution. The solver statistics will be stored in mgD. The breaking parameter defines the slope above which dispersive terms are turned off. The \(\alpha_d\) parameter controls the optimisation of the dispersion relation (see Bonneton et al, 2011).
Definition at line 71 of file green-naghdi.h.
Referenced by diagonalization_2D(), dphidt(), event_tracer_diffusion(), h_relax(), h_residual(), horizontal_diffusion(), implicit_horizontal_diffusion(), relax_GN(), residual_GN(), update_green_naghdi(), and vertical_diffusion().
| mgstats mgD |
Definition at line 72 of file green-naghdi.h.
Referenced by update_green_naghdi().
1.9.8