fluid_problem.c 56.2 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
/*
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
2
 * MigFlow - Copyright (C) <2010-2018>
3
4
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
 * 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
 * List of the contributors to the development of MigFlow: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
7
8
 * Description and complete License: see LICENSE file.
 * 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
9
 * This program (MigFlow) is free software: 
10
 * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
11
12
13
14
15
16
 * Public License as published by the Free Software Foundation, either version
 * 3 of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
 * GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
18
 * 
19
20
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program (see COPYING and COPYING.LESSER files).  If not, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
21
22
23
 * see <http://www.gnu.org/licenses/>.
 */

24
#include "tools.h"
25
26
#include <stdlib.h>
#include <stdio.h>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
27
#include <float.h>
28
29
#include <math.h>
#include "fluid_problem.h"
30
#include "mesh_find.h"
31
#include "mbxml.h"
32

33
#define D DIMENSION
34

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
35
#define  deps 1e-8
36
#if D==2
37
38
39

#define N_SF 3
#define N_N 3
40
#define N_QP 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
41
#define N_LQP 2
42

43
44
static double LQP[] = {0.21132486540518713, 0.7886751345948129};
static double LQW[] = {0.5,0.5};
45
46
static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {1./6,1./6,1./6};
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
static void shape_functions(double *uvw, double *sf)
{
  sf[0] = 1-uvw[0]-uvw[1];
  sf[1] = uvw[0];
  sf[2] = uvw[1];
}
static void grad_shape_functions(const double dxidx[2][2], double gsf[3][2])
{
  for (int i = 0;  i < 2; ++i) {
    gsf[0][i] = -dxidx[0][i]-dxidx[1][i];
    gsf[1][i] = dxidx[0][i];
    gsf[2][i] = dxidx[1][i];
  }
}
static double detDD(const double m[2][2])
{
  return m[0][0]*m[1][1] - m[0][1]*m[1][0];
}
static double invDD(const double m[2][2], double inv[2][2])
{
  double d = detDD(m);
  inv[0][0] = m[1][1]/d;
  inv[0][1] = -m[0][1]/d;
  inv[1][0] = -m[1][0]/d;
  inv[1][1] = m[0][0]/d;
  return d;
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
74
#else
75
76
77

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
78
#define N_QP 5
79
#define N_LQP 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
81
82
83
84
85
86
static double QP[][3] = {
  {0.25, 0.25, 0.25},
  {0.5, 1./6, 1./6},
  {1./6, 0.5, 1./6},
  {1./6, 1./6, 0.5},
  {1./6, 1./6, 1./6}
};
87
static double LQP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
88
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
89
static double LQW[] = {1./6,1./6,1./6};
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
static void shape_functions(double *uvw, double *sf)
{
  sf[0] = 1-uvw[0]-uvw[1]-uvw[2];
  sf[1] = uvw[0];
  sf[2] = uvw[1];
  sf[3] = uvw[2];
}
static void grad_shape_functions(const double dxidx[3][3], double gsf[4][3])
{
  for (int i = 0; i < 3; ++i) {
    gsf[0][i] = -dxidx[0][i]-dxidx[1][i]-dxidx[2][i];
    gsf[1][i] = dxidx[0][i];
    gsf[2][i] = dxidx[1][i];
    gsf[3][i] = dxidx[2][i];
  }
};
static double detDD(const double m[3][3])
{
  return m[0][0]*(m[1][1]*m[2][2]-m[2][1]*m[1][2])
        -m[1][0]*(m[0][1]*m[2][2]-m[2][1]*m[0][2])
        +m[2][0]*(m[0][1]*m[1][2]-m[1][1]*m[0][2]);
}
static double invDD(const double m[3][3], double inv[3][3]) {
  double d = detDD(m);
  inv[0][0] = (m[1][1]*m[2][2]-m[2][1]*m[1][2])/d;
  inv[0][1] = -(m[0][1]*m[2][2]-m[2][1]*m[0][2])/d;
  inv[0][2] = (m[0][1]*m[1][2]-m[1][1]*m[0][2])/d;
  inv[1][0] = -(m[1][0]*m[2][2]-m[2][0]*m[1][2])/d;
  inv[1][1] = (m[0][0]*m[2][2]-m[2][0]*m[0][2])/d;
  inv[1][2] = -(m[0][0]*m[1][2]-m[1][0]*m[0][2])/d;
  inv[2][0] = (m[1][0]*m[2][1]-m[2][0]*m[1][1])/d;
  inv[2][1] = -(m[0][0]*m[2][1]-m[2][0]*m[0][1])/d;
  inv[2][2] = (m[0][0]*m[1][1]-m[1][0]*m[0][1])/d;
  return d;
};

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
126
127
#endif

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
128
129
#define P  D
#define U 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
130
#define A (D+1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
131

132
int fluid_problem_n_fields(const FluidProblem *p) {return D + p->n_fluids;}
133

134
#if D==2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
135
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
136
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
137
  double d = 2*sqrt(vol/M_PI);
138
139
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
140
  double Re_O_u = d*c/mu*rho;
141
142
143
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
  double f = pow(c,-1.8);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
144
  return f*Cd_u*rho/2*d;
145
}
146
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
147
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
148
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
149
  double d = 2*cbrt(3./4.*vol/M_PI);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
150
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
151
152
153
154
155
156
  //Reynolds/|u_p-u_f|
  double Re_O_u = d*c/mu*rho;
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
  double f = pow(c,-1.8);
  double r = d/2.;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
157
  return f*Cd_u*rho/2*(M_PI*r*r);
158
159
}
#endif
160

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
161
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
162
  for (int i = 0; i < problem->n_particles*D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
163
164
165
166
    particle_force[i] = problem->particle_force[i];
  }
}

167
static void particle_force_f(FluidProblem *problem, double *f0, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int ip, int in_derivative, int reduced_gravity) {
168
169

  double p = sol[P];
170
171
172
173
174
175
  double u[D], uold[D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    dp[iD] = dsol[P*D+iD];
  }
176
177
  double Ra, a, mu, rho;

178
  if (problem->n_fluids == 1){
179
180
181
182
183
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    a = sol[A];
184
    a = fmin(1.,fmax(0.,a));
185
186
187
188
189
190
191
192
193
194
195
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
    mu = problem->mu[0]*a + problem->mu[1]*(1-a);
  }

  double Du[D], Due[D];
  for (int j = 0; j < D; ++j) {
    Du[j] = problem->particle_velocity[ip*D+j]-u[j]/c;
    Due[j] = problem->particle_velocity[ip*D+j]-uold[j]/c;
  }
  double mass = problem->particle_mass[ip];
  double vol = problem->particle_volume[ip];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
196
  double g = problem->g;
197
198
199
200
  if (reduced_gravity){
    double rhop = mass/vol;
    g = g*(rhop-problem->rho[0])/rhop;
  }
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
201
  double gamma;
Matthieu Constant's avatar
Matthieu Constant committed
202
  
203
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
204
205
206
207
208
209
210
    double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
    double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
    gamma = gamma0*a + gamma1*(1-a);
  }
  else{
    gamma = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
  }
211
212
  double fcoeff = mass/(mass+dt*gamma);
  for (int j = 0; j < D; ++j){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
213
214
    if(!in_derivative)
      problem->particle_force[ip*D+j] = fcoeff*(-gamma*Du[j]-vol*dp[j]);
215
216
    f0[U+j] = -fcoeff*gamma*(Du[j]-dt/mass*vol*dp[j])-vol*dp[j];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
217
218
  if(!in_derivative)
    problem->particle_force[ip*D+1] += fcoeff*g*mass;
219
220
221
  f0[U+1] += -fcoeff*gamma*dt*g;
}

222
static void node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity)
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
{
  const Mesh *mesh = problem->mesh;
  const double *porosity = problem->porosity;

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *g0 = malloc(sizeof(double)*n_fields);
  double *h0 = malloc(sizeof(double)*n_fields*D);

  for (int i=0;i<n_fields;++i){
    f0[i] = 0;
    g0[i] = 0;
240
241
    for (int j = 0; j < D; ++j)
      h0[i*D+j] = 0;
242
243
244
  }
  for (int ip=0; ip < problem->n_particles; ++ip) {
    int iel = problem->particle_element_id[ip];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
245
    if(iel < 0){
246
      for (int d = 0; d < D; ++d)
247
248
        problem->particle_force[ip*D+d] = 0;
      problem->particle_force[ip*D+1] = problem->g*problem->particle_mass[ip];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
249
250
      continue;
    }
251
    double *xi = problem->particle_uvw + D*ip;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
252
253
    double phi[N_SF];
    shape_functions(xi,phi);
254
255
256
257
258
259
    const unsigned int *el = &mesh->elements[iel*N_N];
    double dxdxi[D][D], dxidx[D][D], dphi[N_SF][D];
    for (int i = 0; i < D; ++i)
      for (int j = 0; j < D; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
    const double detj = invDD(dxdxi, dxidx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
260
    grad_shape_functions(dxidx, dphi);
261
262
263
264
265
266
267
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
    for (int i = 0; i < n_fields; ++i) {
      s[i] = 0;
      sold[i] = 0;
      for (int j = 0; j < D; ++j) {
        ds[i*D+j] = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
268
269
      }
    }
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    double c = 0;
    double cold = 0;
    double dc[D] = {0};
    for (int i = 0; i < N_SF; ++i) {
      c += problem->porosity[el[i]]*phi[i];
      cold += problem->oldporosity[el[i]]*phi[i];
      for (int j = 0; j < n_fields; ++j) {
        double dof = solution[el[i]*n_fields+j];
        double dofold = solution_old[el[i]*n_fields+j];
        s[j] += phi[i]*dof;
        sold[j] += phi[i]*dofold;
        for (int k = 0; k < D; ++k)
          ds[j*D+k] += dphi[i][k]*dof;
      }
      for (int k=0; k<D; ++k){
        dc[k] += problem->porosity[el[i]]*dphi[i][k];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
287
    }
288
    particle_force_f(problem,f0,s,ds,sold,c,dc,cold,dt,iel,ip,0,reduced_gravity);
289
290
291
292
    for (int ifield = 0; ifield < n_fields; ++ifield) {
      for (int iphi = 0; iphi < N_SF; ++iphi){
        local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
293
    }
294
295
296
    for (int jfield = 0; jfield < n_fields; ++jfield) {
      double store = s[jfield];
      s[jfield] += deps;
297
      particle_force_f(problem,g0,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
298
299
300
301
      s[jfield] = store;
      for (int jd =0; jd < D; ++jd){
        store = ds[jfield*D+jd];
        ds[jfield*D+jd] += deps;
302
        particle_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
        ds[jfield*D+jd] = store;
      }
      int N2 = n_fields*N_SF;
#define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)]
      for (int ifield = 0; ifield < n_fields; ++ifield){
        for (int iphi = 0; iphi < N_SF; ++iphi){
          for (int jphi = 0; jphi < N_SF; ++jphi){
            double df00 = (g0[ifield]-f0[ifield])/deps;
            double m = phi[iphi]*phi[jphi]*df00;
            for (int jd = 0; jd < D; ++jd) {
              double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
              m += phi[iphi]*dphi[jphi][jd]*df01;
            }
            int IDX = (iphi*N2)*n_fields + jphi*n_fields;
            LOCAL_MATRIX2(ifield,jfield) += m;
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
319
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
320
321
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
322
  }
323
324
325
326
327
328
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
329
330
}

331
static void f_pressure(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
332
  double p = s[P];
333
334
335
336
337
  double u[D], du[D][D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
338
  }
Matthieu Constant's avatar
Matthieu Constant committed
339

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
340
341
342
  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
343
  double un = 0;
344
  double rho, mu;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
345
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
346
    un += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
347
348
    divu += du[i][i];
    for (int j = 0; j < D; ++j) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
349
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
350
      utau[i] += u[j]*tau[i][j];
Matthieu Constant's avatar
Matthieu Constant committed
351
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
352
  }
Matthieu Constant's avatar
Matthieu Constant committed
353
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
354
  f[P] = un;
355
  if(problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
356
357
358
359
360
    mu = problem->mu[0];
    rho = problem->rho[0];
  }
  else{
    double a = s[A];
361
    a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
362
363
    rho = problem->rho[0]*a+problem->rho[1]*(1-a);
    mu = problem->mu[0]*a+problem->mu[1]*(1-a);
364
    f[A] = un*(un > 0 ? a : v[1]);
Matthieu Constant's avatar
Matthieu Constant committed
365
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
366
  for (int id = 0; id < D; ++id) {
367
    f[U+id] = v[0]*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
368
369
370
  }
}

371
static void f_velocity(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
372
  double p = s[P];
373
374
375
376
377
  double u[D], du[D][D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
378
379
380
381
382
  }

  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
383
  double rho, mu, rho_ext;
Matthieu Constant's avatar
Matthieu Constant committed
384
  double unint = 0;
Matthieu Constant's avatar
Matthieu Constant committed
385
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
386
387
  double un = 0;
  for (int i = 0; i < D; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
388
    un += c*v[i]*n[i];
Matthieu Constant's avatar
Matthieu Constant committed
389
    unint += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
390
391
392
393
    divu += du[i][i];
    for (int j = 0; j < D; ++j) {
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
      utau[i] += u[j]*tau[i][j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
394
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
395
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
396
  f[P] = un;
397
  if(problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
398
399
    mu = problem->mu[0];
    rho = problem->rho[0];
400
    rho_ext = rho;
Matthieu Constant's avatar
Matthieu Constant committed
401
402
403
  }
  else{
    double a = s[A];
404
    a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
405
    rho = problem->rho[0]*a+problem->rho[1]*(1-a);
406
    rho_ext = problem->rho[0]*v[D]+problem->rho[1]*(1-v[D]);
Matthieu Constant's avatar
Matthieu Constant committed
407
    mu = problem->mu[0]*a+problem->mu[1]*(1-a);
408
409
    f[A] = un*(un > 0 ? a : v[D]);
    //printf("%g\n",v[D]);
Matthieu Constant's avatar
Matthieu Constant committed
410
411
  }

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
412
  for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
413
    f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*unint : rho_ext*un*un*n[id]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
415
416
}

417
static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc, const double dt, int eid, const double *v)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
418
419
{
  double p = s[P];
420
421
422
423
424
425
  double u[D], du[D][D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
426
  double un = 0;
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
427
  double a;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
428
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
429
    un += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
430
  }
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
431
432
433
434
435
436
437
438
439
  double rho;
  if (problem->n_fluids == 1){
    rho = problem->rho[0];
  }
  else{
    a = s[A];
    a = fmin(1.,fmax(0.,a));
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
440
441

  f[P] = 0;
Matthieu Constant's avatar
Matthieu Constant committed
442

443
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
444
    double a = s[A];
Matthieu Constant's avatar
Matthieu Constant committed
445
    f[A] = 0;
Matthieu Constant's avatar
Matthieu Constant committed
446
  }
Matthieu Constant's avatar
Matthieu Constant committed
447
448

  for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
449
    f[U+id] = (v==NULL ? p : v[0])*n[id]+ (un> 0 ? rho*u[id]*un/c : 0);
Matthieu Constant's avatar
Matthieu Constant committed
450
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
451
452
453
454
455
}

static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
  const Mesh *mesh = problem->mesh;
456
457
458
459
460
461

  problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
  problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);
  problem->taua = realloc(problem->taua,sizeof(double)*mesh->n_elements);
  problem->extraa = realloc(problem->extraa,sizeof(double)*mesh->n_elements);

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
462
463
464
  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
465
466
  double maxtaa = 0;
  double mintaa = DBL_MAX;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
467
468
469
470
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    double dxdxi[DIMENSION][DIMENSION];
    double normu = {0.};
    const uint32_t *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
471
    double a = 0;
472
473
    double amax = 0;
    double amin = 0;
Matthieu Constant's avatar
Matthieu Constant committed
474
    double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
475
476
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
Matthieu Constant's avatar
Matthieu Constant committed
477
      c += problem->porosity[el[i]]/N_N;
478
      if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
479
        a += problem->solution[el[i]*n_fields+A];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
480
481
        amax = fmax(amax,problem->solution[el[i]*n_fields+A]);
        amin = fmin(amin,problem->solution[el[i]*n_fields+A]);
Matthieu Constant's avatar
Matthieu Constant committed
482
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
483

484
      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
485
486
487
488
489
      normu = fmax(normu,sqrt(normup));
    }
    for (int i=0; i<DIMENSION; ++i){
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
490
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
491
492
493
494
495
#if DIMENSION == 2
    const double h = 2*sqrt(detDD(dxdxi)/(2*M_PI));
#else
    const double h = 2*cbrt(detDD(dxdxi)/(8*M_PI));
#endif
Matthieu Constant's avatar
Matthieu Constant committed
496
    double rho, mu;
497
    if(problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
498
499
500
501
      rho = problem->rho[0];
      mu = problem->mu[0];
    }else{
      a /= N_N;
502
      a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
503
504
505
      rho = problem->rho[0]*a + problem->rho[1]*(1-a);
      mu = problem->mu[0]*a + problem->mu[1]*(1-a);
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
506
    double nu = mu/rho;
Matthieu Constant's avatar
Matthieu Constant committed
507
    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
508
    problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
Matthieu Constant's avatar
Matthieu Constant committed
509
    problem->taua[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h));
510
    problem->extraa[iel] = fmax(0.,fmax(-amin+0.01,amax-1+0.01))*problem->coeffStab;
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
511
512
  }
}
513

Matthieu Constant's avatar
Matthieu Constant committed
514
static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int reduced_gravity, double stab_param) {
515
  double p = sol[P];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
516
517
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
518
  double taua = problem->taua[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
519
  double extraa = problem->extraa[iel];
520
  double u[D], uold[D], du[D][D], dp[D];
521
522
523
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
524
    dp[iD] = dsol[P*D+iD];
525
526
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = dsol[(U+iD)*D+jD];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
527
  double Ra, a, aold, mu, rho, rhoold, drhodt=0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
528

529
  if (problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
530
531
532
533
534
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    a = sol[A];
535
    a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
536
    aold = solold[A];
537
    aold = fmin(1.,fmax(0.,aold));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
538
    Ra = (a*c-aold*cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
539
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
Matthieu Constant's avatar
Matthieu Constant committed
540
    rhoold = problem->rho[0]*aold + problem->rho[1]*(1-aold);
Matthieu Constant's avatar
Matthieu Constant committed
541
542
    mu = problem->mu[0]*a + problem->mu[1]*(1-a);
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
543
544
545
  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
Matthieu Constant's avatar
Matthieu Constant committed
546
547
  double da[D];
  double uugradrho[D] = {0};
Matthieu Constant's avatar
Matthieu Constant committed
548
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
549
550
551
  for (int i = 0; i < D; ++i) {
    divu += du[i][i];
    for (int j = 0; j < D; ++j) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
552
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
553
      utau[i] += u[j]*tau[i][j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
554
    }
555
    if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
556
      da[i] = dsol[A*D+i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
557
      Ra += u[i]*da[i];
558
559
560
561
562
563
    }
  }
  if(problem->n_fluids == 2){
    Ra += a*divu;
    drhodt = (rho-rhoold)/dt;
    for (int i = 0; i < D; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
564
565
566
      for (int j=0;j<D;++j){
        uugradrho[i] += u[i]*u[j]/c *(problem->rho[0]*da[j] - problem->rho[1]*da[j]);
      }
Matthieu Constant's avatar
Matthieu Constant committed
567
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
568
  }
Matthieu Constant's avatar
Matthieu Constant committed
569

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
570
  f0[P] = (c-cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
571

572
  if(problem->n_fluids == 2){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
573
    f0[A] = (a*c-aold*cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
574
575
  }
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
576
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
577
    double dudt = (u[i]-uold[i])/dt;
578
579
    double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==1 ? -(reduced_gravity ? 0 : 1)*c*problem->g : 0.);
    f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -(reduced_gravity ? 0 : 1)*c*rho*problem->g : 0.);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
580
    for (int j = 0; j < D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
581
      f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + (stab_param==0 ? 1 : 0)*R*u[j]*rho*taup;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
582
    }
Matthieu Constant's avatar
Matthieu Constant committed
583
584
    f1[(U+i)*D+i] += (stab_param==0 ? 1 : 0)*divu*tauc*rho;
    f1[P*D+i] = -u[i] + (stab_param==0 ? taup*R : stab_param*dp[i]);
585
    if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
586
      f1[P*D+i] +=  taua*(divu+(c-cold)/dt)*u[i];
Matthieu Constant's avatar
Matthieu Constant committed
587
      f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(extraa);
Matthieu Constant's avatar
Matthieu Constant committed
588
    }
589
  }
Matthieu Constant's avatar
Matthieu Constant committed
590
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
591
    f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
Matthieu Constant's avatar
Matthieu Constant committed
592
  }
593
594
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
595
static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
596
597
598
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
599
600
601
602
603
604
605
606
  const int n_fields = fluid_problem_n_fields(problem);
  const size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *g0 = malloc(sizeof(double)*n_fields);
  double *h0 = malloc(sizeof(double)*n_fields*D);
Matthieu Constant's avatar
Matthieu Constant committed
607
  int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
608
  f_cb *f;
609
610
611
612
613
#if DIMENSION == 2
  int elbnd[3][2] = {{0,1},{1,2},{2,0}};
#else
  int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
614
  for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
Matthieu Constant's avatar
Matthieu Constant committed
615
616
    for (int i =0; i < mesh->n_nodes; ++i)
      bnd_node_local_id[i] = -1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
617
    if(mesh->phys_dimension[iphys] != D-1) continue;
Matthieu Constant's avatar
Matthieu Constant committed
618
619
620
621
622
623
624
625
    double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
      bnd_node_local_id[j] = i;
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
    }
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
626
    f = NULL;
627
628
    int n_value;
    double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
629
    for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
Matthieu Constant's avatar
compile    
Matthieu Constant committed
630
631
632
      WeakBoundary *bnd = problem->weak_boundaries + ibnd;
      if (strcmp(mesh->phys_name[iphys],bnd->tag) == 0){
        f = bnd->cb;
633
634
635
636
637
        n_value = bnd->n_value;
        if (n_value > 0){
          v = realloc(v,mesh->phys_n_nodes[iphys]*sizeof(double)*n_value);
          bnd->apply(mesh->phys_n_nodes[iphys], x, v);
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
638
639
640
641
        break;
      }
    }
    if (f == NULL) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
642
643
      printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
      exit(1);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
644
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
645
646
647
648

    for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
      const int *bnd = &mesh->phys_boundaries[iphys][ibnd*2];
      const int iel = bnd[0];
649
650
651
      const int iebnd = bnd[1];
      int nodes[D];
      double *x[D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
652
      const unsigned int *el = &mesh->elements[iel*N_N];
653
654
655
656
      for (int i = 0; i < D; ++i){
        nodes[i] = el[elbnd[iebnd][i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
657
658
659
660
661
662
663
      double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
      for (int i = 0; i < D; ++i)
        for (int j = 0; j < D; ++j)
          dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
      invDD(dxdxi, dxidx);
      grad_shape_functions(dxidx, dphi);

664
665
666
667
668
669
670
671
672
673
674
#if DIMENSION == 2 
      const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
      const double detbnd = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
      const double n[2] = {dx[1]/detbnd,-dx[0]/detbnd};
#else
      const double dx[2][3] = {{x[1][0]-x[0][0], x[1][1]-x[0][1], x[1][2]-x[0][2]},{x[2][0]-x[0][0], x[2][1]-x[0][1], x[2][2]-x[0][2]}};
      const double nn[3] = {
        dx[0][1]*dx[1][2] - dx[0][2]*dx[1][1],
        dx[0][2]*dx[1][0] - dx[0][0]*dx[1][2],
        dx[0][0]*dx[1][1] - dx[0][1]*dx[1][0]};
      const double detbnd = sqrt(nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
675
      const double n[3] = {-nn[0]/detbnd,-nn[1]/detbnd,-nn[2]/detbnd};
676
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
677
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
678
679
680
681
682
        double xi[D];
#if DIMENSION == 2 
        double a = LQP[iqp];
        switch(iebnd) {
          case 0 : xi[0] = a;   xi[1] = 0;   break;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
683
684
685
686
          case 1 : xi[0] = 1-a; xi[1] = a;   break;
          case 2 : xi[0] = 0;   xi[1] = 1-a; break;
          default:;
        }
687
688
689
690
691
692
693
694
695
696
697
#else
        double a = LQP[iqp][0];
        double b = LQP[iqp][1];
        switch(iebnd) {
          case 0 : xi[0] = a;     xi[1] = b; xi[2] = 0; break;
          case 1 : xi[0] = 0;     xi[1] = a; xi[2] = b; break;
          case 2 : xi[0] = b;     xi[1] = 0; xi[2] = a; break;
          case 3 : xi[0] = 1-a-b; xi[1] = b; xi[2] = a; break;
          default:;
        }
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
698
699
700
701
702
703
704
705
706
707
708
709
710
        double dp[D]={0};
        double *local_vector = &all_local_vector[local_size*iel];
        double *local_matrix = &all_local_matrix[local_size*local_size*iel];
        double phi[N_SF];
        shape_functions(xi,phi);
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
          sold[i] = 0;
          for (int j = 0; j < D; ++j) {
            ds[i*D+j] = 0;
          }
        }
        double c = 0;
Matthieu Constant's avatar
Matthieu Constant committed
711
712
713
714
        double *vbnd = NULL;
        if (n_value>0){
          vbnd = malloc(sizeof(double)*n_value);
        }
715
716
717
        for (int i = 0; i < n_value; ++i){
          vbnd[i] = 0;
        }
Matthieu Constant's avatar
Matthieu Constant committed
718
        double dc[D] = {0};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
719
        for (int i = 0; i < N_SF; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
720
          c += problem->porosity[el[i]]*phi[i];
Matthieu Constant's avatar
Matthieu Constant committed
721
722
          int node_local_id = bnd_node_local_id[el[i]];
          if (node_local_id != -1)
723
724
725
            for (int j = 0; j < n_value; j++){
              vbnd[j] += v[node_local_id*n_value+j]*phi[i];
            }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
726
727
728
729
730
731
          for (int j = 0; j < n_fields; ++j) {
            double dof = solution[el[i]*n_fields+j];
            s[j] += phi[i]*dof;
            for (int k = 0; k < D; ++k)
              ds[j*D+k] += dphi[i][k]*dof;
          }
Matthieu Constant's avatar
Matthieu Constant committed
732
733
734
          for (int k=0; k<D; ++k){
            dc[k] += problem->porosity[el[i]]*dphi[i][k];
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
735
        }
736
        const double jw = LQW[iqp]*detbnd;
Matthieu Constant's avatar
Matthieu Constant committed
737
        f(problem,n,f0,s,ds,c,dc,dt,iel,vbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
738
739
740
741
742
743
744
745
        for (int ifield = 0; ifield < n_fields; ++ifield) {
          for (int iphi = 0; iphi < N_SF; ++iphi){
            local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
          }
        }
        for (int jfield = 0; jfield < n_fields; ++jfield) {
          double store = s[jfield];
          s[jfield] += deps;
Matthieu Constant's avatar
Matthieu Constant committed
746
          f(problem,n,g0,s,ds,c,dc,dt,iel,vbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
747
748
749
750
          s[jfield] = store;
          for (int jd =0; jd < D; ++jd){
            store = ds[jfield*D+jd];
            ds[jfield*D+jd] += deps;
Matthieu Constant's avatar
Matthieu Constant committed
751
            f(problem,n,h0+jd*n_fields,s,ds,c,dc,dt,iel,vbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
            ds[jfield*D+jd] = store;
          }
          int N2 = n_fields*N_SF;
#define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)]
          for (int ifield = 0; ifield < n_fields; ++ifield){
            for (int iphi = 0; iphi < N_SF; ++iphi){
              for (int jphi = 0; jphi < N_SF; ++jphi){
                double df00 = (g0[ifield]-f0[ifield])/deps;
                double m = jw*phi[iphi]*phi[jphi]*df00;
                for (int jd = 0; jd < D; ++jd) {
                  double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
                  m += jw*phi[iphi]*dphi[jphi][jd]*df01;
                }
                int IDX = (iphi*N2)*n_fields + jphi*n_fields;
                LOCAL_MATRIX2(ifield,jfield) += m;
              }
            }
          }
        }
771
        free(vbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
772
      }
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
773
    }
Matthieu Constant's avatar
Matthieu Constant committed
774
775
    free(x);
    free(v);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
776
  }
777
778
779
780
781
782
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
783
}
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
784

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
785

Matthieu Constant's avatar
Matthieu Constant committed
786
static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity, double stab_param)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
787
788
{
  const Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
789
  const double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
790
791
792
793

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
794
795
796
797
798
799
800
801
802
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *f1 = malloc(sizeof(double)*n_fields*D);
  double *g0 = malloc(sizeof(double)*n_fields);
  double *g1 = malloc(sizeof(double)*n_fields*D);
  double *h0 = malloc(sizeof(double)*n_fields*D);
  double *h1 = malloc(sizeof(double)*n_fields*D*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
803
  for (int iel=0; iel < mesh->n_elements; ++iel) {
804
    const unsigned int *el = &mesh->elements[iel*N_N];
805
806
807
    double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
    for (int i = 0; i < D; ++i)
      for (int j = 0; j < D; ++j)
808
809
810
        dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
    const double detj = invDD(dxdxi, dxidx);
    grad_shape_functions(dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
811
812
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
813
814
815
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
816
817
818
819
820
821
822
      for (int i = 0; i < n_fields; ++i) {
        s[i] = 0;
        sold[i] = 0;
        for (int j = 0; j < D; ++j) {
          ds[i*D+j] = 0;
        }
      }
823
      double c = 0;
Matthieu Constant's avatar
Matthieu Constant committed
824
825
      double cold = 0;
      double dc[D] = {0};
826
      for (int i = 0; i < N_SF; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
827
828
        c += problem->porosity[el[i]]*phi[i];
        cold += problem->oldporosity[el[i]]*phi[i];
829
830
831
832
833
834
835
836
        for (int j = 0; j < n_fields; ++j) {
          double dof = solution[el[i]*n_fields+j];
          double dofold = solution_old[el[i]*n_fields+j];
          s[j] += phi[i]*dof;
          sold[j] += phi[i]*dofold;
          for (int k = 0; k < D; ++k)
            ds[j*D+k] += dphi[i][k]*dof;
        }
Matthieu Constant's avatar
Matthieu Constant committed
837
838
839
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
840
841
      }
      const double jw = QW[iqp]*detj;
Matthieu Constant's avatar
Matthieu Constant committed
842
      fluid_problem_f(problem,f0,f1,s,ds,sold,c,dc,cold,dt,iel, reduced_gravity, stab_param);
843
844
845
846
847
      for (int ifield = 0; ifield < n_fields; ++ifield) {
        for (int iphi = 0; iphi < N_SF; ++iphi){
          local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
          for (int id = 0; id < D; ++id) {
            local_vector[iphi+N_SF*ifield] += dphi[iphi][id]*f1[ifield*D+id]*jw;
Matthieu Constant's avatar
Matthieu Constant committed
848
849
          }
        }
850
851
852
853
      }
      for (int jfield = 0; jfield < n_fields; ++jfield) {
        double store = s[jfield];
        s[jfield] += deps;
Matthieu Constant's avatar
Matthieu Constant committed
854
        fluid_problem_f(problem,g0,g1,s,ds,sold,c,dc,cold,dt,iel,reduced_gravity,stab_param);
855
856
857
858
        s[jfield] = store;
        for (int jd =0; jd < D; ++jd){
          store = ds[jfield*D+jd];
          ds[jfield*D+jd] += deps;
Matthieu Constant's avatar
Matthieu Constant committed
859
          fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dc,cold,dt,iel,reduced_gravity,stab_param);
860
          ds[jfield*D+jd] = store;
861
        }
862
863
864
865
866
867
868
869
870
871
872
873
874
        int N2 = n_fields*N_SF;
#define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)]
        for (int ifield = 0; ifield < n_fields; ++ifield){
          for (int iphi = 0; iphi < N_SF; ++iphi){
            for (int jphi = 0; jphi < N_SF; ++jphi){
              double df00 = (g0[ifield]-f0[ifield])/deps;
              double m = jw*phi[iphi]*phi[jphi]*df00;
              for (int id = 0; id < D; ++id) {
                double df10 = (g1[ifield*D+id]-f1[ifield*D+id])/deps;
                m += jw*dphi[iphi][id]*phi[jphi]*df10;
                for (int jd = 0; jd < D; ++jd) {
                  double df11 = (h1[jd*n_fields*D+ifield*D+id]-f1[ifield*D+id])/deps;
                  m += jw*dphi[iphi][id]*dphi[jphi][jd]*df11;
875
876
                }
              }
877
878
879
              for (int jd = 0; jd < D; ++jd) {
                double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
                m += jw*phi[iphi]*dphi[jphi][jd]*df01;
Matthieu Constant's avatar
Matthieu Constant committed
880
              }
881
882
              int IDX = (iphi*N2)*n_fields + jphi*n_fields;
              LOCAL_MATRIX2(ifield,jfield) += m;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
883
            }
Matthieu Constant's avatar
Matthieu Constant committed
884
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
885
886