# hub.darcs.net :: basilisk -> basilisk -> files

Basilisk source code

## 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; }