fluid_problem.c 56.1 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

32
#define D DIMENSION
33

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

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

42
43
static double LQP[] = {0.21132486540518713, 0.7886751345948129};
static double LQW[] = {0.5,0.5};
44
45
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};
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
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
73
#else
74
75
76

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
77
#define N_QP 5
78
#define N_LQP 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
79
80
81
82
83
84
85
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}
};
86
static double LQP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
87
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
88
static double LQW[] = {1./6,1./6,1./6};
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
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
125
126
#endif

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

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

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

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

166
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) {
167
168

  double p = sol[P];
169
170
171
172
173
174
  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];
  }
175
176
  double Ra, a, mu, rho;

177
  if (problem->n_fluids == 1){
178
179
180
181
182
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    a = sol[A];
183
    a = fmin(1.,fmax(0.,a));
184
185
186
187
188
189
190
191
192
193
194
    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
195
  double g = problem->g;
196
197
198
199
  if (reduced_gravity){
    double rhop = mass/vol;
    g = g*(rhop-problem->rho[0])/rhop;
  }
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
200
  double gamma;
Matthieu Constant's avatar
Matthieu Constant committed
201
  
202
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
203
204
205
206
207
208
209
    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);
  }
210
211
  double fcoeff = mass/(mass+dt*gamma);
  for (int j = 0; j < D; ++j){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
212
213
    if(!in_derivative)
      problem->particle_force[ip*D+j] = fcoeff*(-gamma*Du[j]-vol*dp[j]);
214
215
    f0[U+j] = -fcoeff*gamma*(Du[j]-dt/mass*vol*dp[j])-vol*dp[j];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
216
217
  if(!in_derivative)
    problem->particle_force[ip*D+1] += fcoeff*g*mass;
218
219
220
  f0[U+1] += -fcoeff*gamma*dt*g;
}

221
static void node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity)
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
{
  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;
239
240
    for (int j = 0; j < D; ++j)
      h0[i*D+j] = 0;
241
242
243
  }
  for (int ip=0; ip < problem->n_particles; ++ip) {
    int iel = problem->particle_element_id[ip];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
244
    if(iel < 0){
245
      for (int d = 0; d < D; ++d)
246
247
        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
248
249
      continue;
    }
250
    double *xi = problem->particle_uvw + D*ip;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
251
252
    double phi[N_SF];
    shape_functions(xi,phi);
253
254
255
256
257
258
    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
259
    grad_shape_functions(dxidx, dphi);
260
261
262
263
264
265
266
    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
267
268
      }
    }
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
    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
286
    }
287
    particle_force_f(problem,f0,s,ds,sold,c,dc,cold,dt,iel,ip,0,reduced_gravity);
288
289
290
291
    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
292
    }
293
294
295
    for (int jfield = 0; jfield < n_fields; ++jfield) {
      double store = s[jfield];
      s[jfield] += deps;
296
      particle_force_f(problem,g0,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
297
298
299
300
      s[jfield] = store;
      for (int jd =0; jd < D; ++jd){
        store = ds[jfield*D+jd];
        ds[jfield*D+jd] += deps;
301
        particle_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
        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
318
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
319
320
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
321
  }
322
323
324
325
326
327
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
328
329
}

330
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
331
  double p = s[P];
332
333
334
335
336
  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];
337
  }
Matthieu Constant's avatar
Matthieu Constant committed
338

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

370
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
371
  double p = s[P];
372
373
374
375
376
  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
377
378
379
380
381
  }

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
411
  for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
412
    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
413
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
415
}

416
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
417
418
{
  double p = s[P];
419
420
421
422
423
424
  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
425
  double un = 0;
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
426
  double a;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
427
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
428
    un += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
429
  }
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
430
431
432
433
434
435
436
437
438
  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
439
440

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

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

  for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
448
    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
449
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
450
451
452
453
454
}

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

  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
461
462
463
  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
464
465
  double maxtaa = 0;
  double mintaa = DBL_MAX;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
466
467
468
469
  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
470
    double a = 0;
Nathan Coppin's avatar
Nathan Coppin committed
471
472
    double amax = 0.5;
    double amin = 0.5;
Matthieu Constant's avatar
Matthieu Constant committed
473
    double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
474
475
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
Matthieu Constant's avatar
Matthieu Constant committed
476
      c += problem->porosity[el[i]]/N_N;
477
      if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
478
        a += problem->solution[el[i]*n_fields+A];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
479
480
        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
481
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
482

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

Matthieu Constant's avatar
Matthieu Constant committed
513
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) {
514
  double p = sol[P];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
515
516
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
517
  double taua = problem->taua[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
518
  double extraa = problem->extraa[iel];
519
  double u[D], uold[D], du[D][D], dp[D];
520
521
522
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
523
    dp[iD] = dsol[P*D+iD];
524
525
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = dsol[(U+iD)*D+jD];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
526
  double Ra, a, aold, mu, rho, rhoold, drhodt=0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
527

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

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

571
  if(problem->n_fluids == 2){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
572
    f0[A] = (a*c-aold*cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
573
574
  }
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
575
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
576
    double dudt = (u[i]-uold[i])/dt;
577
578
    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
579
    for (int j = 0; j < D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
580
      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
581
    }
Matthieu Constant's avatar
Matthieu Constant committed
582
583
    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]);
584
    if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
585
      f1[P*D+i] +=  taua*(divu+(c-cold)/dt)*u[i];
Matthieu Constant's avatar
Matthieu Constant committed
586
      f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(extraa);
Matthieu Constant's avatar
Matthieu Constant committed
587
    }
588
  }
Matthieu Constant's avatar
Matthieu Constant committed
589
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
590
    f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
Matthieu Constant's avatar
Matthieu Constant committed
591
  }
592
593
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
594
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
595
596
597
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
598
599
600
601
602
603
604
605
  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
606
  int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
607
  f_cb *f;
608
609
610
611
612
#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
613
  for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
Matthieu Constant's avatar
Matthieu Constant committed
614
615
    for (int i =0; i < mesh->n_nodes; ++i)
      bnd_node_local_id[i] = -1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
616
    if(mesh->phys_dimension[iphys] != D-1) continue;
Matthieu Constant's avatar
Matthieu Constant committed
617
618
619
620
621
622
623
624
    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
625
    f = NULL;
626
627
    int n_value;
    double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
628
    for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
Matthieu Constant's avatar
compile    
Matthieu Constant committed
629
630
631
      WeakBoundary *bnd = problem->weak_boundaries + ibnd;
      if (strcmp(mesh->phys_name[iphys],bnd->tag) == 0){
        f = bnd->cb;
632
633
634
635
636
        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
637
638
639
640
        break;
      }
    }
    if (f == NULL) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
641
642
      printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
      exit(1);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
643
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
644
645
646
647

    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];
648
649
650
      const int iebnd = bnd[1];
      int nodes[D];
      double *x[D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
651
      const unsigned int *el = &mesh->elements[iel*N_N];
652
653
654
655
      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
656
657
658
659
660
661
662
      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);

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
784

Matthieu Constant's avatar
Matthieu Constant committed
785
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
786
787
{
  const Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
788
  const double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
789
790
791
792

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