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

root / src / conservation.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
274
275
276
277
278
279
280
281
282
283
/**
# A generic solver for systems of conservation laws

Using the ideas of [Kurganov and Tadmor,
2000](references.bib#kurganov2000) it is possible to write a generic
solver for systems of conservation laws of the form
$$
\partial_t\left(\begin{array}{c}
    s_i\\
    \mathbf{v}_j\\
 \end{array}\right) + \nabla\cdot\left(\begin{array}{c}
    \mathbf{F}_i\\
    \mathbf{T}_j\\
 \end{array}\right) = 0
$$
where $s_i$ is a list of scalar fields, $\mathbf{v}_j$ a list of
vector fields and $\mathbf{F}_i$, $\mathbf{T}_j$ are the corresponding
vector (resp. tensor) fluxes. 

Note that the [Saint-Venant solver](saint-venant.h) is a particular
case of this generic algorithm.

The user must provide the lists of conserved scalar and vector fields
*/

extern scalar * scalars;
extern vector * vectors;

/**
as well as a function which, given the state of each quantity,
returns the fluxes and the minimum/maximum eigenvalues
(i.e. `eigenvalue[0]`/`eigenvalue[1]`) of the corresponding Riemann
problem. */

void flux (const double * state, double * flux, double * eigenvalue);

/**
## Time-integration

### Setup

Time integration will be done with a generic
[predictor-corrector](predictor-corrector.h) scheme. */

#include "predictor-corrector.h"

/**
The generic time-integration scheme in `predictor-corrector.h` needs
to know which fields are updated i.e. all the scalars + the components
of all the vector fields. It also needs a function to compute the
time-derivatives of the evolving variables. */

scalar * evolving;
double update_conservation (scalar * conserved, scalar * updates, double dtmax);

event defaults (i = 0)
{
  evolving = list_concat (scalars, (scalar *) vectors);
  update = update_conservation;

  /**
  We switch to a pure minmod limiter by default for increased
  robustness. */
  
  theta = 1.;

  /**
  With the MUSCL scheme we use the CFL depends on the dimension of the
  problem. */

  if (CFL > 1./dimension)
    CFL = 1./dimension;
  
  /**
  On trees we need to replace the default bilinear
  refinement/prolongation with linear so that reconstructed values
  also use slope limiting. */

  #if TREE
  for (scalar s in evolving) {
    s.refine = s.prolongation = refine_linear;
    s.restriction = restriction_volume_average;
  }
  #endif
}

/**
At the end of the run we need to free the list (to avoid a memory
leak). */

event cleanup (i = end) free (evolving);

/**
We force boundary conditions on all fields after initialisation.
*/

event init (i = 0)
{
  boundary (all);
}

/**
### Computing fluxes

The core of the central-upwind scheme (see e.g. section 3.1 of
[Kurganov & Levy, 2002](references.bib#kurganov2002)) is the
approximate solution of the Riemann problem given by the left and
right states to get the fluxes `f`. */

static double riemann (const double * right, const double * left,
		       double Delta, double * f, int len, 
		       double dtmax)
{
  double fr[len], fl[len], er[2], el[2];
  flux (right, fr, er);
  flux (left,  fl, el);
  double ap = max(er[1], el[1]); ap = max(ap, 0.);
  double am = min(er[0], el[0]); am = min(am, 0.);
  double a = max(ap, -am); 

  if (a > 0.) {
    for (int i = 0; i < len; i++)
      f[i] = (ap*fl[i] - am*fr[i] + ap*am*(right[i] - left[i]))/(ap - am);
    double dt = CFL*Delta/a;
    if (dt < dtmax)
      dtmax = dt;
  }
  else
    for (int i = 0; i < len; i++)
      f[i] = 0.;
  return dtmax;
}

double update_conservation (scalar * conserved, scalar * updates, double dtmax)
{
  /**
  The gradients of each quantity are stored in a list of dynamically-allocated
  fields. First-order reconstruction is used for the gradient fields. */

  vector * slopes = NULL;
  for (scalar s in conserved) {
    vector slope = new vector;
    foreach_dimension() {
      slope.x.gradient = zero;
      #if TREE
      slope.x.prolongation = refine_linear;
      #endif
    }
    slopes = vectors_append (slopes, slope);
  }
  gradients (conserved, slopes);

  /**
  We allocated fields for storing fluxes for each scalar and vector
  quantity. */

  vector * lflux = NULL;
  int len = list_len (conserved);
  for (scalar s in conserved) {
    vector f1 = new vector;
    lflux = vectors_append (lflux, f1);
  }

  /**
  The predictor-corrector scheme treats all fields as scalars (stored in
  the `conserved` list). We need to recover vector and tensor quantities
  from these lists. To do so, knowing the number of scalar fields, we
  split the scalar list into a list of scalars and a list of vectors. */
  
  int scalars_len = list_len (scalars);

  scalar * scalars = list_copy (conserved);
  if (scalars) scalars[scalars_len].i = -1;
  vector * vectors = vectors_from_scalars (&conserved[scalars_len]);
  
  /**
  We then do the same for the gradients i.e. split the list of vectors
  into a list of vectors and a list of tensors. */

  vector * scalar_slopes = vectors_copy (slopes);
  if (scalar_slopes) scalar_slopes[scalars_len] = (vector){{-1}};
  tensor * vector_slopes = tensors_from_vectors (&slopes[scalars_len]);

  /**
  And again for the fluxes. */
  
  vector * scalar_fluxes = vectors_copy (lflux);
  if (scalar_fluxes) scalar_fluxes[scalars_len] = (vector){{-1}};
  tensor * vector_fluxes = tensors_from_vectors (&lflux[scalars_len]);

  /**
  We are ready to compute the fluxes through each face of the domain. */
  
  foreach_face (reduction (min:dtmax)) {

    /**
    #### Left/right state reconstruction 
    
    We use the central values of each scalar/vector quantity and the
    pre-computed gradients to compute the left and right states. */
    
    double r[len], l[len]; // right/left Riemann states
    double f[len];         // fluxes for each conserved quantity
    double dx = Delta/2.;
    int i = 0;
    scalar s;
    vector g;
    for (s,g in scalars,scalar_slopes) {
      r[i] = s[] - dx*g.x[];
      l[i++] = s[-1] + dx*g.x[-1];
    }
    vector v;
    tensor t;
    for (v,t in vectors,vector_slopes) {      
      r[i] = v.x[] - dx*t.x.x[];
      l[i++] = v.x[-1] + dx*t.x.x[-1];
      #if dimension > 1
        r[i] = v.y[] - dx*t.y.x[];
	l[i++] = v.y[-1] + dx*t.y.x[-1];
      #endif
      #if dimension > 2
        r[i] = v.z[] - dx*t.z.x[];
	l[i++] = v.z[-1] + dx*t.z.x[-1];
      #endif
    }

    /**
    #### Riemann problem
    
    We then call the generic Riemann solver and store the resulting fluxes
    in the pre-allocated fields. */
    
    dtmax = riemann (r, l, Delta*cm[]/fm.x[], f, len, dtmax);
    i = 0;
    for (vector fs in scalar_fluxes)
      fs.x[] = fm.x[]*f[i++];
    for (tensor fv in vector_fluxes) {
      fv.x.x[] = fm.x[]*f[i++];
      #if dimension > 1
        fv.y.x[] = fm.x[]*f[i++];
      #endif
      #if dimension > 2
        fv.z.x[] = fm.x[]*f[i++];
      #endif
    }
  }

  boundary_flux (lflux);

  /**
  #### Update
  
  The update for each scalar quantity is the divergence of the fluxes. */
  
  foreach() {
    scalar ds;
    vector f;
    for (ds,f in updates,lflux) {
      ds[] = 0.;
      foreach_dimension()
	ds[] += (f.x[] - f.x[1])/(cm[]*Delta);
    }
  }

  /**
  #### Cleanup
  
  We finally deallocate the memory used to store lists and gradient
  fields. */
  
  free (scalars);
  free (vectors);
  free (scalar_slopes);
  free (vector_slopes);
  free (scalar_fluxes);
  free (vector_fluxes);
  delete ((scalar *) slopes);
  free (slopes);
  delete ((scalar *) lflux);
  free (lflux);
  
  return dtmax;
}