fluid_problem.c 56.5 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"
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
31
#include "fluid_problem_fem.h"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
32

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
34
35
#define P  D
#define U 0

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
36
37
typedef void f_cb(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double rho, const double mu, const double dt, int eid, const double *v);

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
38
int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}
39

40
#if D==2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
41
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
42
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
43
  double d = 2*sqrt(vol/M_PI);
44
45
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
46
  double Re_O_u = d*c/mu*rho;
47
48
49
  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
50
  return f*Cd_u*rho/2*d;
51
}
52
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
53
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
54
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
55
  double d = 2*cbrt(3./4.*vol/M_PI);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
56
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
57
58
59
60
61
62
  //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
63
  return f*Cd_u*rho/2*(M_PI*r*r);
64
65
}
#endif
66

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
67
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
68
  for (int i = 0; i < problem->n_particles*D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
69
70
71
72
    particle_force[i] = problem->particle_force[i];
  }
}

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
73
static void particle_force_f(FluidProblem *problem, double *f0, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip, int in_derivative, int reduced_gravity) {
74
75
76
77
78
79
  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];
  }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
80
  double mu, rho;
81

82
  if (problem->n_fluids == 1){
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    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
98
  double g = problem->g;
99
100
101
102
  if (reduced_gravity){
    double rhop = mass/vol;
    g = g*(rhop-problem->rho[0])/rhop;
  }
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
103
  double gamma;
Matthieu Constant's avatar
Matthieu Constant committed
104
  
105
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
106
107
108
109
110
111
112
    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);
  }
113
114
  double fcoeff = mass/(mass+dt*gamma);
  for (int j = 0; j < D; ++j){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
115
116
    if(!in_derivative)
      problem->particle_force[ip*D+j] = fcoeff*(-gamma*Du[j]-vol*dp[j]);
117
118
    f0[U+j] = -fcoeff*gamma*(Du[j]-dt/mass*vol*dp[j])-vol*dp[j];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
119
120
  if(!in_derivative)
    problem->particle_force[ip*D+1] += fcoeff*g*mass;
121
122
123
  f0[U+1] += -fcoeff*gamma*dt*g;
}

124
static void node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity)
125
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
126

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
  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;
143
144
    for (int j = 0; j < D; ++j)
      h0[i*D+j] = 0;
145
146
147
  }
  for (int ip=0; ip < problem->n_particles; ++ip) {
    int iel = problem->particle_element_id[ip];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
    if(iel < 0){
149
      for (int d = 0; d < D; ++d)
150
151
        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
152
153
      continue;
    }
154
    double *xi = problem->particle_uvw + D*ip;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
155
156
    double phi[N_SF];
    shape_functions(xi,phi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
157
    const int *el = &mesh->elements[iel*N_N];
158
159
160
161
162
    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
163
    grad_shape_functions(dxidx, dphi);
164
165
166
167
168
169
170
    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
171
172
      }
    }
173
    double c = 0;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
174
    double a = 0;
175
176
177
178
179
    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];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
180
181
182
      if (problem->n_fluids == 2) {
        a += problem->concentration[iel*N_N+i]*phi[i];
      }
183
184
185
186
187
188
189
190
191
192
193
      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
194
    }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
195
    particle_force_f(problem,f0,s,ds,sold,c,dc,a,dt,iel,ip,0,reduced_gravity);
196
197
198
199
    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
200
    }
201
202
203
    for (int jfield = 0; jfield < n_fields; ++jfield) {
      double store = s[jfield];
      s[jfield] += deps;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
204
      particle_force_f(problem,g0,s,ds,sold,c,dc,a,dt,iel,ip,1,reduced_gravity);
205
206
207
208
      s[jfield] = store;
      for (int jd =0; jd < D; ++jd){
        store = ds[jfield*D+jd];
        ds[jfield*D+jd] += deps;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
209
        particle_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,a,dt,iel,ip,1,reduced_gravity);
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
        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
226
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
227
228
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
229
  }
230
231
232
233
234
235
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
236
237
}

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
238
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 rho, const double mu, const double dt, int eid, const double *v){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
239
  double p = s[P];
240
241
242
243
244
  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];
245
  }
Matthieu Constant's avatar
Matthieu Constant committed
246

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
247
248
249
  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
250
  double un = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
251
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
252
    un += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
253
254
    divu += du[i][i];
    for (int j = 0; j < D; ++j) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
255
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
256
      utau[i] += u[j]*tau[i][j];
Matthieu Constant's avatar
Matthieu Constant committed
257
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
258
  }
Matthieu Constant's avatar
Matthieu Constant committed
259
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
260
  f[P] = un;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
261
  for (int id = 0; id < D; ++id) {
262
    f[U+id] = v[0]*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
263
264
265
  }
}

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
266
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 rho, const double mu, const double dt, int eid, const double *v){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
267
  double p = s[P];
268
269
270
271
272
  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
273
274
275
276
277
  }

  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
278
  double rho_ext;
Matthieu Constant's avatar
Matthieu Constant committed
279
  double unint = 0;
Matthieu Constant's avatar
Matthieu Constant committed
280
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
281
282
  double un = 0;
  for (int i = 0; i < D; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
283
    un += c*v[i]*n[i];
Matthieu Constant's avatar
Matthieu Constant committed
284
    unint += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
285
286
287
288
    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
289
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
290
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
291
  f[P] = un;
292
  if(problem->n_fluids == 1){
293
    rho_ext = rho;
Matthieu Constant's avatar
Matthieu Constant committed
294
295
  }
  else{
296
    rho_ext = problem->rho[0]*v[D]+problem->rho[1]*(1-v[D]);
Matthieu Constant's avatar
Matthieu Constant committed
297
298
  }

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
299
  for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
300
    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
301
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
302
303
}

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
304
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 rho, const double mu, const double dt, int eid, const double *v)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
305
306
{
  double p = s[P];
307
308
309
310
311
312
  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
313
  double un = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
314
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
315
    un += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
316
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
317
318
  double h = problem->element_size[eid];
  f[P] = (v==NULL ? 0 : -(s[P]-v[0])/h*problem->taup[eid]);
Matthieu Constant's avatar
Matthieu Constant committed
319

Matthieu Constant's avatar
Matthieu Constant committed
320
  for (int id = 0; id < D; ++id) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
321
    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
322
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
323
324
325
326
327
}

static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
  const Mesh *mesh = problem->mesh;
328

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
329
  problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
330
331
332
  problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
  problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
333
334
335
  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
336
337
  double maxtaa = 0;
  double mintaa = DBL_MAX;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
338
339
340
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    double dxdxi[DIMENSION][DIMENSION];
    double normu = {0.};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
341
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
342
    double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
343
344
    double amax = 0.5;
    double amin = 0.5;
Matthieu Constant's avatar
Matthieu Constant committed
345
    double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
346
347
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
Matthieu Constant's avatar
Matthieu Constant committed
348
      c += problem->porosity[el[i]]/N_N;
349
      if(problem->n_fluids == 2){
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
350
351
352
        a += problem->concentration[iel*N_N+i];
        amax = fmax(amax,problem->concentration[iel*N_N+i]);
        amin = fmin(amin,problem->concentration[iel*N_N+i]);
Matthieu Constant's avatar
Matthieu Constant committed
353
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
354

355
      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
356
357
358
359
360
      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
361
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
362
363
364
365
366
#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
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
367
    problem->element_size[iel] = h;
Matthieu Constant's avatar
Matthieu Constant committed
368
    double rho, mu;
369
    if(problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
370
371
372
373
      rho = problem->rho[0];
      mu = problem->mu[0];
    }else{
      a /= N_N;
374
      a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
375
376
377
      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
378
    double nu = mu/rho;
Matthieu Constant's avatar
Matthieu Constant committed
379
    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
380
    problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
381
382
  }
}
383

Matthieu Constant's avatar
Matthieu Constant committed
384
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, const double rho, const double *drho, const double mu, double dt, int iel, int reduced_gravity, double stab_param) {
385
  double p = sol[P];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
386
387
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
388
  double u[D], uold[D], du[D][D], dp[D];
389
390
391
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
392
    dp[iD] = dsol[P*D+iD];
393
394
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = dsol[(U+iD)*D+jD];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
395

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
396
397
398
  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
399
  double uudrho[D] = {0};
Matthieu Constant's avatar
Matthieu Constant committed
400
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
401
402
403
  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
404
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
405
      utau[i] += u[j]*tau[i][j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
406
    }
407
  }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
408
409
410
411
412
413
  double udrho = 0;
  for (int j=0;j<D;++j){
    udrho+= u[j]/c*drho[j];
  }
  for (int i = 0; i < D; ++i) {
    uudrho[i] += u[i]*udrho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
  }
Matthieu Constant's avatar
Matthieu Constant committed
415

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
418
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
419
    double dudt = (u[i]-uold[i])/dt;
420
421
    double R = dudt + u[i]*udrho/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uudrho[i]/rho + (i==1 ? -(reduced_gravity ? (rho-problem->rho[0])/rho : 1)*c*problem->g : 0.);
    f0[U+i] = rho*dudt + u[i]*udrho + (i==1 ? -(reduced_gravity ? (rho-problem->rho[0])/rho : 1)*c*rho*problem->g : 0.);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
422
    for (int j = 0; j < D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
423
      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
424
    }
Matthieu Constant's avatar
Matthieu Constant committed
425
426
    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]);
Matthieu Constant's avatar
Matthieu Constant committed
427
  }
428
429
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
430
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
431
432
433
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
434
435
436
437
438
439
440
441
  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
442
  int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
443
444
445
446
447
448
449
450
451
452
453
454
455
456
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    double *v = malloc(wbnd->n_value*bnd->n_nodes*sizeof(double));
    for (int i = 0; i < mesh->n_nodes; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
457
      bnd_node_local_id[i] = -1;
Matthieu Constant's avatar
Matthieu Constant committed
458
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
459
460
461
462
463
464
465
    if (wbnd->n_value != 0) {
      double *x = (double*)malloc(bnd->n_nodes*sizeof(double)*D);
      for (int i = 0; i < bnd->n_nodes; ++i){
        int j = bnd->nodes[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
466
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
467
468
      wbnd->apply(bnd->n_nodes, x, v);
      free(x);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
469
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
470
471
472
473
474
475
476
477
478
    f_cb *bndcb;
    switch(wbnd->type) {
      case BND_WALL : bndcb = f_wall; break;
      case BND_PRESSURE : bndcb = f_pressure; break;
      case BND_VELOCITY : bndcb = f_velocity; break;
      default :
        printf("Unknown boundary type : %i\n", wbnd->type);
        exit(1);
    };
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
479
480
481
482
483
    double *vbnd = malloc(sizeof(double)*wbnd->n_value);
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
484
485
      int nodes[D];
      double *x[D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
486
      const int *el = &mesh->elements[iel*N_N];
487
      for (int i = 0; i < D; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
488
        nodes[i] = el[elbnd[icl][i]];
489
490
        x[i] = &mesh->x[nodes[i]*3];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
491
492
493
494
495
      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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
496
497
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
498
499
500
      grad_shape_functions(dxidx, dphi);

      for (int iqp = 0; iqp < N_LQP; ++iqp) {
501
        double xi[D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
502
        {
503
#if DIMENSION == 2 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
504
          double a = LQP[iqp][0];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
505
506
507
508
509
510
          switch(icl) {
            case 0 : xi[0] = a;   xi[1] = 0;   break;
            case 1 : xi[0] = 1-a; xi[1] = a;   break;
            case 2 : xi[0] = 0;   xi[1] = 1-a; break;
            default:;
          }
511
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
512
513
514
515
516
517
518
519
520
          double a = LQP[iqp][0];
          double b = LQP[iqp][1];
          switch(icl) {
            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:;
          }
521
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
522
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
523
524
525
526
527
528
529
530
531
532
533
534
535
        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;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
536
        for (int i = 0; i < wbnd->n_value; ++i){
537
538
          vbnd[i] = 0;
        }
Matthieu Constant's avatar
Matthieu Constant committed
539
        double dc[D] = {0};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
540
        double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
541
        for (int i = 0; i < N_SF; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
542
          c += problem->porosity[el[i]]*phi[i];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
543
544
545
          if (problem->n_fluids == 2) {
            a += problem->concentration[iel*N_N+i]*phi[i];
          }
Matthieu Constant's avatar
Matthieu Constant committed
546
547
          int node_local_id = bnd_node_local_id[el[i]];
          if (node_local_id != -1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
548
549
            for (int j = 0; j < wbnd->n_value; j++){
              vbnd[j] += v[node_local_id*wbnd->n_value+j]*phi[i];
550
            }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
551
552
553
554
555
556
          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
557
558
559
          for (int k=0; k<D; ++k){
            dc[k] += problem->porosity[el[i]]*dphi[i][k];
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
560
        }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
561
562
563
564
565
        double rho = problem->rho[0], mu = problem->mu[0];
        if (problem->n_fluids == 2) {
          rho = problem->rho[0]*a + problem->rho[1]*(1-a);
          mu = problem->mu[0]*a + problem->mu[1]*(1-a);
        }
566
        const double jw = LQW[iqp]*detbnd;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
567
        bndcb(problem,n,f0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0 ? NULL :vbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
568
569
570
571
572
573
574
575
        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;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
576
          bndcb(problem,n,g0,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value == 0? NULL :vbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
577
578
579
580
          s[jfield] = store;
          for (int jd =0; jd < D; ++jd){
            store = ds[jfield*D+jd];
            ds[jfield*D+jd] += deps;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
581
            bndcb(problem,n,h0+jd*n_fields,s,ds,c,dc,rho,mu,dt,iel,wbnd->n_value ==0? NULL : vbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
            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;
              }
            }
          }
        }
      }
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
602
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
603
    free(vbnd);
Matthieu Constant's avatar
Matthieu Constant committed
604
    free(v);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
605
  }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
606
  free(bnd_node_local_id);
607
608
609
610
611
612
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
613
}
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
614

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
615

Matthieu Constant's avatar
Matthieu Constant committed
616
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
617
618
{
  const Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
619
  const double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
620
621
622
623

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
624
625
626
627
628
629
630
631
632
  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
633
  for (int iel=0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
634
    const int *el = &mesh->elements[iel*N_N];
635
636
637
    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)
638
639
640
        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
641
642
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
643
644
645
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
646
647
648
649
650
651
652
      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
broken    
Jonathan Lambrechts committed
653
654
655
      double c = 0, a = 0;
      double cold = 0, aold = 0;
      double dc[D] = {0}, da[D]={0};
656
      for (int i = 0; i < N_SF; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
657
658
        c += problem->porosity[el[i]]*phi[i];
        cold += problem->oldporosity[el[i]]*phi[i];
659
660
661
662
663
664
665
666
        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
667
668
669
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
        if (problem->n_fluids==2) {
          a += problem->concentration[iel*N_N+i]*phi[i];
          for (int k=0; k<D; ++k){
            da[k] += problem->concentration[iel*N_N+i]*dphi[i][k];
          }
        }
      }
      double mu = problem->mu[0];
      double rho = problem->rho[0];
      double drho[D] = {0};
      if(problem->n_fluids == 2) {
        mu = a*problem->mu[0] + (1-a)*problem->mu[1];
        rho = a*problem->rho[0] + (1-a)*problem->rho[1];
        for (int k=0; k<D; ++k){
          drho[k] = da[k]*(problem->rho[0]-problem->rho[1]);
        }
686
687
      }
      const double jw = QW[iqp]*detj;
Matthieu Constant's avatar
Matthieu Constant committed
688
      fluid_problem_f(problem,f0,f1,s,ds,sold,c,dc,cold,rho,drho,mu,dt,iel, reduced_gravity, stab_param);
689
690
691
692
693
      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
694
695
          }
        }
696
697
698
699
      }
      for (int jfield = 0; jfield < n_fields; ++jfield) {
        double store = s[jfield];
        s[jfield] += deps;
Matthieu Constant's avatar
Matthieu Constant committed
700
        fluid_problem_f(problem,g0,g1,s,ds,sold,c,dc,cold,rho,drho,mu,dt,iel,reduced_gravity,stab_param);
701
702
703
704
        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
705
          fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dc,cold,rho,drho,mu,dt,iel,reduced_gravity,stab_param);
706
          ds[jfield*D+jd] = store;
707
        }
708
709
710
711
712
713
714
715
716
717
718
719
720
        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;
721
722
                }
              }
723
724
725
              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
726
              }
727
728
              int IDX = (iphi*N2)*n_fields + jphi*n_fields;
              LOCAL_MATRIX2(ifield,jfield) += m;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
729
            }
Matthieu Constant's avatar
Matthieu Constant committed
730
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
731
732
733
        }
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
734
735
736
737
738
739
740
741
742
743
744
745
746
  }

  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(f1);
  free(g0);
  free(g1);
  free(h0);
  free(h1);
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
747
748
749
750
751
752
753
754
755
756
757
758
static void fluid_problem_surface_tension(FluidProblem *problem, double *all_local_vector){
  if (problem->n_fluids == 1) return;
  const Mesh *mesh = problem->mesh;
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxdxi[D][D],dxidx[D][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);
  }

Matthieu Constant's avatar
Matthieu Constant committed
759
  double sigma = 22.27*1e-3;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
  int n_fields = fluid_problem_n_fields(problem);
  double *a_cg = malloc(sizeof(double)*mesh->n_nodes);
  for (int i = 0; i< mesh->n_nodes; ++i) {
    a_cg[i] = 0;
  }
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxdxi[D][D],dxidx[D][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];
    double det = invDD(dxdxi, dxidx);
#if DIMENSION == 2
    double vol = det/6;
#else
    double vol = det/24;
#endif
    for (int i = 0; i< N_N; ++i) {
      a_cg[el[i]] += vol * problem->concentration[iel*N_N+i];
    }
  }
  for (int i = 0; i< mesh->n_nodes; ++i) {
    a_cg[i] /= problem->node_volume[i];
  }
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxdxi[D][D],dxidx[D][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];
Matthieu Constant's avatar
Matthieu Constant committed
790
791
792
793
794
795
    double det = invDD(dxdxi, dxidx);
#if DIMENSION == 2
    double vol = det/2;
#else
    double vol = det/6;
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
796
797
    double dphi[N_N][D];
    grad_shape_functions(dxidx, dphi);
798
    double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
799
800
    double da[D] = {0};
    double nda = 0;
801
802
803
804
    for (int i = 0; i < N_SF; ++i) {
      a += a_cg[el[i]];
    }
    a /= N_N;
Matthieu Constant's avatar
Matthieu Constant committed
805
806
807
808
    double c = 0;
    for (int i = 0; i < N_N; ++i){
      c += problem->porosity[el[i]]/N_N;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
809
810
811
812
813
814
815
    for (int k=0; k<D; ++k){
      for (int i = 0; i < N_SF; ++i) {
        da[k] += a_cg[el[i]]*dphi[i][k];
      }
      nda += da[k]*da[k];
    }
    nda = fmax(sqrt(nda),1e-8);
816
    double *local_vector = &all_local_vector[iel*n_fields*N_N];
817
818
819
    for (int iphi = 0; iphi < N_N; ++iphi) {
      for (int id = 0; id < D; ++id) {
        for (int jd = 0; jd < D; ++jd) {
Matthieu Constant's avatar
Matthieu Constant committed
820
          local_vector[(U+id)*N_N+iphi] -= c*vol*da[id]*da[jd]*dphi[iphi][jd]/nda*sigma;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
821
        }
Matthieu Constant's avatar
Matthieu Constant committed
822
        local_vector[(U+id)*N_N+iphi] += c*vol*nda*dphi[iphi][id]*sigma;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
823
824
825
826
827
828
      }
    }
  }
  free(a_cg);
}

Matthieu Constant's avatar
Matthieu Constant committed
829
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt, int reduced_gravity, double stab_param)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
{
  HXTLinearSystem *lsys = problem->linear_system;
  const Mesh *mesh = problem->mesh;

  int n_fields = fluid_problem_n_fields(problem);
  
  size_t local_size = N_SF*n_fields;
  size_t N_E = mesh->n_elements;
  double *all_local_vector = (double*) malloc(N_E*local_size*sizeof(double));
  double *all_local_matrix = (double*) malloc(N_E*local_size*local_size*sizeof(double));
  for (size_t i = 0; i < N_E*local_size; ++i)
    all_local_vector[i] = 0;
  for (size_t i = 0; i < N_E*local_size*local_size; ++i)
    all_local_matrix[i] = 0;

  char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
  for (int i = 0; i < n_fields*mesh->n_nodes; ++i)
    forced_row[i] = -1;
  for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
    const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
    int iphys;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
851
852
    for (iphys = 0; iphys < mesh->n_boundaries; ++iphys) {
      if (strcmp(bnd->tag, mesh->boundary_names[iphys]) == 0)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
853
854
        break;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
855
    if (iphys == mesh->n_boundaries){
856
      printf("Boundary tag \"%s\" not found.\n", bnd->tag);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
857
858
859
860
      exit(1);
    }
    for (int i = 0; i < problem->boundaries[iphys].n_nodes; ++i){
      forced_row[problem->boundaries[iphys].nodes[i]*n_fields+bnd->row] = bnd->field;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
861
862
    }
  }
863
  node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity);
Matthieu Constant's avatar
Matthieu Constant committed
864
  compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
Matthieu Constant's avatar
Matthieu Constant committed
865
  fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity,stab_param);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
866
  fluid_problem_surface_tension(problem, all_local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
867
868
869
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
870
871
    for (int inode = 0; inode < N_SF; ++inode) {
      for(int ifield = 0; ifield < n_fields; ++ifield) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
872
        const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
873
874
875
876
877
878
879
880
881
882
883
884
        char forced_field = forced_row[el[inode]*n_fields+ifield];
        if (forced_field == -1) continue;
        int i = inode*n_fields + ifield;
        for (int jnode = 0; jnode< N_SF; ++jnode){
          for (int jfield= 0; jfield< n_fields; ++jfield){
            local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] =
              (inode==jnode && jfield == forced_field) ? 1. : 0.;
          }
        }
        local_vector[inode+ifield*N_SF] = 0;
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
885
886
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
887
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
888
  free(forced_row);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
889
890
  free(all_local_matrix);
  free(all_local_vector);
891
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
892

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
893
static void apply_boundary_conditions(FluidProblem *problem)
894
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
895
896
897
898
  const Mesh *mesh = problem->mesh;
  int n_fields = fluid_problem_n_fields(problem);
  for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
    StrongBoundary *bnd = problem->strong_boundaries + ibnd;
899
    int iphys;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
900
901
    for (iphys = 0; iphys < mesh->n_boundaries; ++iphys) {
      if (strcmp(bnd->tag, mesh->boundary_names[iphys]) == 0)
902
903
        break;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
904
    if (iphys == mesh->n_boundaries){
905
      printf("Boundary tag \"%s\" not found.\n", bnd->tag);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
906
      exit(1);
907
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
908
909
910
911
912
    MeshBoundary *mbnd = &problem->boundaries[iphys];
    double *x = (double*)malloc(mbnd->n_nodes*sizeof(double)*D);
    double *v = (double*)malloc(mbnd->n_nodes*sizeof(double));
    for (int i = 0; i < mbnd->n_nodes; ++i){
      int j = mbnd->nodes[i];
913
914
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
915
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
916
917
918
    bnd->apply(mbnd->n_nodes, x, v);
    for (int i = 0; i < mbnd->n_nodes; ++i)
      problem->solution[mbnd->nodes[i]*n_fields+