Basilisk source code (http://basilisk.fr/src/)

root / src / all-mach.h.page

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
/**
# An "all Mach" flow solver

We wish to solve the generic momentum equation
$$
\partial_t\mathbf{q} + \nabla\cdot(\mathbf{q}\mathbf{u}) = 
- \nabla p + \nabla\cdot(\mu\nabla\mathbf{u}) + \rho\mathbf{a}
$$
with $\mathbf{q}=\rho\mathbf{u}$ the momentum, $\mathbf{u}$ the
velocity, $\rho$ the density, $\mu$ the dynamic viscosity, $p$ the
pressure and $\mathbf{a}$ an
acceleration. The pressure is defined through an equation of state and
verifies the evolution equation
$$
\partial_t p + \mathbf{u}\cdot\nabla p = -\rho c^2\nabla\cdot\mathbf{u}
$$
with $c$ the speed of sound. By default the solver sets $c=\infty$,
$\rho=1$ and the pressure equation reduces to 
$$
\nabla\cdot\mathbf{u} = 0
$$

The advection of momentum is not performed by this solver (so that
different schemes can be used) i.e. in the end, by default, we solve
the incompressible (linearised) Euler equations with a projection
method.

We build the solver using the generic time loop and the implicit
viscous solver (which includes the multigrid Poisson--Helmholtz
solver). */

#include "run.h"
#include "timestep.h"
#include "viscosity.h"

/**
The primitive variables are the momentum $\mathbf{q}$, pressure $p$,
density $\rho$, (face) specific volume $\alpha=1/\rho$, (face) dynamic
viscosity $\mu$ (which is zero by default) and (face) velocity field
$\mathbf{u}_f$. */

vector q[];
scalar p[];
face vector uf[];
(const) face vector alpha = unityf, mu = zerof;
(const) scalar rho = unity;

/**
The equation of state is defined by the pressure field *ps* and $\rho
c^2$. In the incompressible limit $\rho c^2\rightarrow\infty$. Rather
than trying to compute this limit, we set both fields to zero and
check for this particular case when computing the pressure (see
below). This means that by default the fluid is incompressible. */

scalar ps[];
(const) scalar rhoc2 = zeroc;
(const) face vector a = zerof;

/**
We store the combined pressure gradient and acceleration field in
*g*. */

vector g[];

event defaults (i = 0) {

  /**
  The default density field is set to unity (times the metric). */

  if (alpha.x.i == unityf.x.i)
    alpha = fm;
}

/**
We apply boundary conditions and define the face velocity field after
user initialisation. */

event init (i = 0) {
  boundary ({q,p,rho});

  /**
  The face velocity field is obtained by simple linear interpolation
  from the momentum field. We make sure that the specific volume
  $\alpha$ is defined by calling the "properties" event (see
  below). */
  
  event ("properties");
  foreach_face()
    uf.x[] = alpha.x[]*(q.x[] + q.x[-1])/2.;
  boundary ((scalar *){uf});
}

/**
The timestep is computed using the CFL condition on the face velocity
field. */

double dtmax;

event set_dtmax (i++,last) dtmax = DT;

event stability (i++,last) {
  dt = dtnext (timestep (uf, dtmax));
}

/**
Tracers (including momentum $\mathbf{q}$) are advected by these events. */

event vof (i++,last);
event tracer_advection (i++,last);

/**
The equation of state (i.e. fields $\alpha$, $\rho$, $\rho c^2$ and
*ps*) is defined by this event. */

event properties (i++,last)
{
  boundary ({rho, rhoc2, ps, alpha, mu});
  
  /**
  If the acceleration is not constant, we reset it to zero. */
  
  if (!is_constant(a.x)) {
    face vector af = a;
    foreach_face()
      af.x[] = 0.;
  }
}

/**
This event can be overloaded to add acceleration terms. */

event acceleration (i++, last)
{
  boundary ((scalar *){a});
}

/**
The equation for the pressure is a Poisson--Helmoltz problem which we
will solve with the [multigrid solver](poisson.h). The statistics for
the solver will be stored in *mgp* (resp. *mgu* for the viscosity
solver). */

mgstats mgp, mgu;

event pressure (i++, last)
{

  /**
  If the viscosity is not zero, we use the implicit viscosity solver
  to obtain the velocity field at time $t + \Delta t$. The pressure
  term is taken into account using the pressure gradient at time
  $t$. A provisionary momentum (without the pressure gradient) is then
  built from the velocity field. */
  
  if (constant(mu.x) != 0.) {
    foreach()
      foreach_dimension()
        q.x[] = (q.x[] + dt*g.x[])/rho[];
    boundary ((scalar *){q});
    mgu = viscosity (q, mu, rho, dt, mgu.nrelax);
    foreach()
      foreach_dimension()
        q.x[] = q.x[]*rho[] - dt*g.x[];
    boundary ((scalar *){q});
  }  
  
  /**
  We first define a temporary face velocity field $\mathbf{u}_\star$
  using simple averaging from $\mathbf{q}_{\star}$, $\alpha_{n+1}$ and
  the acceleration term. */

  foreach_face()
    uf.x[] = alpha.x[]*(q.x[] + q.x[-1])/2. + dt*fm.x[]*a.x[];
  boundary ((scalar *){uf});

  /**
  The evolution equation for the pressure is
  $$\partial_tp +\mathbf{u} \cdot \nabla p = - \rho c^2 \nabla \cdot \mathbf{u}$$
  with $\rho$ the density and $c$ the speed of sound. Following the
  classical [projection
  method](navier-stokes/centered.h#approximate-projection) for
  incompressible flows, we set
  $$
  \mathbf{u}_{n + 1} = \mathbf{u}_{\star} - \Delta t (\alpha\nabla p)_{n+1}
  $$
  The evolution equation for the pressure can then be discretised as
  $$
  \frac{p_{n + 1} - p_n}{\Delta t} +\mathbf{u}_n \cdot \nabla p_n = 
     - \rho c^2_{n + 1} \nabla \cdot \mathbf{u}_{n + 1}
  $$
  which gives, after some manipulations, the Poisson--Helmholtz equation
  $$
  \lambda_{n + 1} p_{n + 1} + \nabla \cdot \left( \alpha \nabla p
  \right)_{n + 1} = \lambda_{n + 1} p_{\star} + \frac{1}{\Delta t} \nabla \cdot
  \mathbf{u}_{\star}
  $$
  with
  $$
  p_{\star} = p_n - \Delta t\mathbf{u}_n \cdot \nabla p_n
  $$
  and
  $$
  \lambda = \frac{- 1}{\Delta t^2 \rho c^2}
  $$
  */

  scalar lambda = rhoc2, rhs = ps;
  foreach() {

    /**
    We compute $\lambda$ and the corresponding term in the
    right-hand-side of the Poisson--Helmholtz equation. */

    if (constant(lambda) == 0.)
      rhs[] = 0.;
    else {
      lambda[] = - cm[]/(sq(dt)*rhoc2[]);
      rhs[] = lambda[]*ps[];
    }
      
    /**
    We add the divergence of the velocity field to the right-hand-side. */

    double div = 0.;
    foreach_dimension()
      div += uf.x[1] - uf.x[];
    rhs[] += div/(dt*Delta);
  }
  boundary ({lambda, rhs});
  
  /**
  The Poisson--Helmholtz solver is called with a [definition of the
  tolerance](poisson.h#project) identical to that used for
  incompressible flows. */
  
  mgp = poisson (p, rhs, alpha, lambda, tolerance = TOLERANCE/sq(dt));

  /**
  The pressure gradient is applied to $\mathbf{u}_\star$ to obtain the
  face velocity field at time $n + 1$. 
  
  We also compute the combined face pressure gradient and acceleration
  field *gf*. */

  face vector gf[];
  foreach_face() {
    double dp = alpha.x[]*(p[] - p[-1])/Delta;
    uf.x[] -= dt*dp;
    gf.x[] = a.x[] - dp/fm.x[];
  }
  boundary_flux ({gf});

  /**
  And finally we apply the pressure gradient/acceleration term to the
  flux/momentum. We also store the centered, combined pressure
  gradient and acceleration field *g*. */
  
  foreach()
    foreach_dimension() {
      g.x[] = rho[]*(gf.x[] + gf.x[1])/2.;
      q.x[] += dt*g.x[];
    }
  boundary ((scalar *){q,g,uf});
}

/**
After mesh adaptation fluid properties need to be updated. */

#if TREE
event adapt (i++,last) {
  event ("properties");
}
#endif