fluid_problem.c 48.8 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
/*
 * Marblesbag - Copyright (C) <2010-2018>
3
4
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
 * 	
6
 * List of the contributors to the development of Marblesbag: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
7
8
9
 * Description and complete License: see LICENSE file.
 * 	
 * This program (Marblesbag) 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
Jonathan Lambrechts committed
31
#include "mbxml.h"
32

Matthieu Constant's avatar
Matthieu Constant committed
33

34

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
35
#define newton_max_it 10
Matthieu Constant's avatar
Matthieu Constant committed
36
37
#define newton_atol 5e-6
#define newton_rtol 1e-4
38
#define D DIMENSION
39

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
40
#define  deps 1e-8
41
#if D==2
42
43
44

#define N_SF 3
#define N_N 3
45
#define N_QP 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
46
#define N_LQP 2
47

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
48
49
static double LQP[] = {-5.773502691896257e-01, 5.773502691896257e-01};
static double LQW[] = {1,1};
50
51
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};
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
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
79
#else
80
81
82

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
83
84
85
86
87
88
89
90
91
#define N_QP 5
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}
};
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
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
128
129
#endif

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

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

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

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

169
170
171
172
173
174
175
176
177
static void node_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) {

  double p = sol[P];
  double dp[D] = {dsol[P*D+0],dsol[P*D+1]};
  double u[D] = {sol[U],sol[U+1]};
  double uold[D] = {solold[U],solold[U+1]};
  double Ra, a, mu, rho;
  //printf("c=%e, dc[0]=%e, dc[1]=%e, cold=%e, dt=%e\n",c,dc[0],dc[1],cold,dt);

178
  if (problem->n_fluids == 1){
179
180
181
182
183
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    a = sol[A];
184
    a = fmin(1.,fmax(0.,a));
185
186
187
188
189
190
191
192
193
194
195
196
    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];
  double g = problem->g;
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
197
  double gamma;
198
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
199
200
201
202
203
204
205
    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);
  }
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
  double fcoeff = mass/(mass+dt*gamma);

  for (int j = 0; j < D; ++j){
    problem->particle_force[ip*D+j] = fcoeff*(-gamma*Du[j]-vol*dp[j]);
    f0[U+j] = -fcoeff*gamma*(Du[j]-dt/mass*vol*dp[j])-vol*dp[j];
  }
  problem->particle_force[ip*D+1] += fcoeff*g*mass;
  f0[U+1] += -fcoeff*gamma*dt*g;
}

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
326
static void f_outflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
327
  double p = s[P];
328
329
330
331
  double dp[D];
  for (int id = 0; id < D; ++id) {
    dp[id] = ds[P*D+id];
  }
Matthieu Constant's avatar
Matthieu Constant committed
332
  double mu, rho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
333
334
  double u[D] = {s[U],s[U+1]};
  double du[D][D] = {{ds[U*D+0],ds[U*D+1]},{ds[U*D+D+0],ds[U*D+D+1]}};
Matthieu Constant's avatar
Matthieu Constant committed
335

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

static void f_inflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double *dc,const double dt, int eid){
  double p = s[P];
  double dp[D];
  for (int id = 0; id < D; ++id) {
    dp[id] = ds[P*D+id];
  }
  double u[D] = {s[U],s[U+1]};
  double du[D][D] = {{ds[U*D+0],ds[U*D+1]},{ds[U*D+D+0],ds[U*D+D+1]}};

  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
Matthieu Constant's avatar
Matthieu Constant committed
377
378
  double rho, mu;
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
379
380
381
382
383
384
385
  double un = 0;
  for (int i = 0; i < D; ++i) {
    un += u[i]*n[i];
    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
386
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
387
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
388
  f[P] = un;
Matthieu Constant's avatar
Matthieu Constant committed
389

390
  if(problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
391
392
393
394
395
    mu = problem->mu[0];
    rho = problem->rho[0];
  }
  else{
    double a = s[A];
396
    a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
397
398
    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
399
    f[A] = un*(un > 0 ? a : 1.);
Matthieu Constant's avatar
Matthieu Constant committed
400
401
  }

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
402
  for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
403
    f[U+id] = p*n[id]+ (un> 0 ? rho*u[id]*un/c : 1);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
404
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
405
406
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
407
static void f_wall_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)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
408
409
410
411
412
413
{
  double dp[D];
  for (int id = 0; id < D; ++id) {
    dp[id] = ds[P*D+id];
  }
  double p = s[P];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
415
416
  double mu = problem->mu[0];
  double u[D] = {s[U],s[U+1]};
  double du[D][D] = {{ds[U*D+0],ds[U*D+1]},{ds[U*D+D+0],ds[U*D+D+1]}};
Matthieu Constant's avatar
Matthieu Constant committed
417

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
418
  double un = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
419
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
420
    un += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
421
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
422
423

  f[P] = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
424
  for (int id = 0; id < D; ++id) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
425
426
    f[U+id] = p*n[id];
  }
Matthieu Constant's avatar
Matthieu Constant committed
427

428
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
429
    double a = s[A];
Matthieu Constant's avatar
Matthieu Constant committed
430
    f[A] = 0;
Matthieu Constant's avatar
Matthieu Constant committed
431
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
432
433
434
435
436
}

static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
  const Mesh *mesh = problem->mesh;
437
438
439
440
441
442

  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
443
444
445
446
447
448
449
  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
  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
450
    double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
451
452
    double amax = 0;
    double amin = 0;
Matthieu Constant's avatar
Matthieu Constant committed
453
    double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
454
455
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
Matthieu Constant's avatar
Matthieu Constant committed
456
      c += problem->porosity[el[i]]/N_N;
457
      if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
458
        a += problem->solution[el[i]*n_fields+A];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
459
460
        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
461
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
462

463
      for(int j = 0; j < DIMENSION; ++j){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
464
        normup += pow2(problem->solution[el[i]*n_fields+j]);
465
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
466
467
468
469
470
      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
471
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
472
473
474
475
476
#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
477
    double rho, mu;
478
    if(problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
479
480
481
482
      rho = problem->rho[0];
      mu = problem->mu[0];
    }else{
      a /= N_N;
483
      a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
484
485
486
      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
487
    double nu = mu/rho;
Matthieu Constant's avatar
Matthieu Constant committed
488
    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
489
    problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
Matthieu Constant's avatar
Matthieu Constant committed
490
491
    problem->taua[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h));
    problem->extraa[iel] = fmax(0.,fmax(-amin+0.01,amax-1+0.01))*0.01;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
492
493
    maxtaup = fmax(maxtaup, problem->taup[iel]);
    mintaup = fmin(mintaup, problem->taup[iel]);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
494
495
  }
}
496

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
497
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) {
498
499
  double p = sol[P];
  double dp[D] = {dsol[P*D+0],dsol[P*D+1]};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
500
501
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
502
  double taua = problem->taua[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
503
  double extraa = problem->extraa[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
504
505
506
  double u[D] = {sol[U],sol[U+1]};
  double uold[D] = {solold[U],solold[U+1]};
  double du[D][D] = {{dsol[U*D+0],dsol[U*D+1]},{dsol[U*D+D+0],dsol[U*D+D+1]}};
Matthieu Constant's avatar
Matthieu Constant committed
507
  double Ra, a, aold, mu, rho, rhoold, drhodt;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
508

509
  if (problem->n_fluids == 1){
Matthieu Constant's avatar
Matthieu Constant committed
510
511
512
513
514
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    a = sol[A];
515
    a = fmin(1.,fmax(0.,a));
Matthieu Constant's avatar
Matthieu Constant committed
516
    aold = solold[A];
517
    aold = fmin(1.,fmax(0.,aold));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
518
    Ra = (a*c-aold*cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
519
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
Matthieu Constant's avatar
Matthieu Constant committed
520
    rhoold = problem->rho[0]*aold + problem->rho[1]*(1-aold);
Matthieu Constant's avatar
Matthieu Constant committed
521
522
    mu = problem->mu[0]*a + problem->mu[1]*(1-a);
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
523
524
525
  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
Matthieu Constant's avatar
Matthieu Constant committed
526
527
  double da[D];
  double uugradrho[D] = {0};
Matthieu Constant's avatar
Matthieu Constant committed
528
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
529
530
531
  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
532
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
533
      utau[i] += u[j]*tau[i][j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
534
    }
535
    if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
536
      da[i] = dsol[A*D+i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
537
      Ra += u[i]*da[i];
538
539
540
541
542
543
    }
  }
  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
544
545
546
      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
547
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
548
  }
Matthieu Constant's avatar
Matthieu Constant committed
549

Matthieu Constant's avatar
Matthieu Constant committed
550
  f0[P] = (c-cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
551

552
  if(problem->n_fluids == 2){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
553
    f0[A] = (a*c-aold*cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
554
555
  }
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
556
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
557
    double dudt = (u[i]-uold[i])/dt;
Matthieu Constant's avatar
Matthieu Constant committed
558
559
    double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==(D-1) ? -c*problem->g : 0);
    f0[U+i] = rho*dudt + u[i]*drhodt + (i==(D-1) ? -c*rho*problem->g : 0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
560
    for (int j = 0; j < D; ++j) {
561
      f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + R*u[j]*rho*taup;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
562
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
563
    f1[(U+i)*D+i] += divu*tauc*rho;
Matthieu Constant's avatar
Matthieu Constant committed
564
    f1[P*D+i] = -u[i] + taup*R;
565
    if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
566
567
      f1[P*D+i] +=  taua*(divu+(c-cold)/dt)*u[i];
      f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(1e-7+extraa);
Matthieu Constant's avatar
Matthieu Constant committed
568
    }
569
  }
570
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
571
572
    f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
  }
573
574
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
575
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
576
577
578
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
579
580
581
582
583
584
585
586
587
588
  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);
  f_cb *f;
  for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
589
590
591
592
593
594
595
596
597
598
599
600
601
    if(mesh->phys_dimension[iphys] != D-1) continue;
    f = NULL;
    for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
      WeakBoundary bnd = problem->weak_boundaries[ibnd];
      if (strcmp(mesh->phys_name[iphys],bnd.tag) == 0){
        f = bnd.cb;
        break;
      }
    }
    if (f == NULL) {
      printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
      exit(1);
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620

    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];
      const int i0 = bnd[1];
      const int i1 = (bnd[1]+1) % 3;
      const unsigned int *el = &mesh->elements[iel*N_N];

      const int nodes[2] = {el[i0],el[i1]};
      const double *x[2] = {&mesh->x[nodes[0]*3], &mesh->x[nodes[1]*3]};
      const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};

      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);
      const double l = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
621
      const double n[2] = {dx[1]/l,-dx[0]/l};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644

      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double xi[2];
        double a = (1+LQP[iqp])/2;
        switch(i0) {
          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:;
        }
        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
645
        double dc[D] = {0};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
646
        for (int i = 0; i < N_SF; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
647
          c += problem->porosity[el[i]]*phi[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
648
649
650
651
652
653
          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
654
655
656
          for (int k=0; k<D; ++k){
            dc[k] += problem->porosity[el[i]]*dphi[i][k];
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
657
658
        }
        const double jw = LQW[iqp]*l/2;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
659
        f(problem,n,f0,s,ds,c,dc,dt,iel);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
660
661
662
663
664
665
666
667
        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
668
          f(problem,n,g0,s,ds,c,dc,dt,iel);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
669
670
671
672
          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
673
            f(problem,n,h0+jd*n_fields,s,ds,c,dc,dt,iel);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
            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
694
695
    }
  }
696
697
698
699
700
701
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
702
}
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
703

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
704
705
706
707

static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
  const Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
708
  const double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
709
710
711
712

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
713
714
715
716
717
718
719
720
721
  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
722
  for (int iel=0; iel < mesh->n_elements; ++iel) {
723
    const unsigned int *el = &mesh->elements[iel*N_N];
724
725
726
    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)
727
728
729
        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
730
731
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
732
733
734
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
735
736
737
738
739
740
741
      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;
        }
      }
742
      double c = 0;
Matthieu Constant's avatar
Matthieu Constant committed
743
744
      double cold = 0;
      double dc[D] = {0};
745
      for (int i = 0; i < N_SF; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
746
747
        c += problem->porosity[el[i]]*phi[i];
        cold += problem->oldporosity[el[i]]*phi[i];
748
749
750
751
752
753
754
755
        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
756
757
758
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
759
760
      }
      const double jw = QW[iqp]*detj;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
761
      fluid_problem_f(problem,f0,f1,s,ds,sold,c,dc,cold,dt,iel);
762
763
764
765
766
      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
767
768
          }
        }
769
770
771
772
      }
      for (int jfield = 0; jfield < n_fields; ++jfield) {
        double store = s[jfield];
        s[jfield] += deps;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
773
        fluid_problem_f(problem,g0,g1,s,ds,sold,c,dc,cold,dt,iel);
774
775
776
777
        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
778
          fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dc,cold,dt,iel);
779
          ds[jfield*D+jd] = store;
780
        }
781
782
783
784
785
786
787
788
789
790
791
792
793
        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;
794
795
                }
              }
796
797
798
              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
799
              }
800
801
              int IDX = (iphi*N2)*n_fields + jphi*n_fields;
              LOCAL_MATRIX2(ifield,jfield) += m;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
802
            }
Matthieu Constant's avatar
Matthieu Constant committed
803
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
804
805
806
        }
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
  }

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

static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
  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;
    for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
      if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_phys)
847
      printf("Boundary tag \"%s\" not found.\n", bnd->tag);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
848
849
850
851
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
    }
  }
852
  node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
853
854
855
856
857
  compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
  fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
  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
858
859
    for (int inode = 0; inode < N_SF; ++inode) {
      for(int ifield = 0; ifield < n_fields; ++ifield) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
860
        const unsigned int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
861
862
863
864
865
866
867
868
869
870
871
872
        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
873
874
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
875
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
876
  free(forced_row);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
877
878
  free(all_local_matrix);
  free(all_local_vector);
879
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
880

881
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int n_fields, int nBnd, StrongBoundary *bnds)
882
{
883
884
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
885
886
887
888
889
890
    int iphys;
    for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
      if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_phys){
891
      printf("Boundary tag \"%s\" not found.\n", bnd->tag);
892
    }
893
894
    double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
    double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
895
896
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
897
898
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
899
900
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
901
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
902
903
904
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
905
906
907
  }
}

908
int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
909
{
910
911
  static double totaltime = 0;
  struct timespec tic, toc;
912
  const Mesh *mesh = problem->mesh;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
913
  int n_fields = fluid_problem_n_fields(problem);
914
915
916
  apply_boundary_conditions(mesh, problem->solution, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
  double *solution_old = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
  double *rhs = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
917
918
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
919
  problem->solution_explicit = solution_old;
920
  double *dx = (double*)malloc(sizeof(double)*mesh->n_nodes*n_fields);
921
  double norm0 = 0;
922
  int i;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
923
  compute_stab_parameters(problem,dt);
924
  for (i=0; i<newton_max_it; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
925
926
    double norm = 0;
    hxtLinearSystemZeroMatrix(problem->linear_system);
927
928

    for (int i=0; i<mesh->n_nodes*n_fields; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
929
      rhs[i] = 0;
930
931
932
      //printf("solution=%g\n",problem->solution[i]);
    }
    //break;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
933
934
    fluid_problem_assemble_system(problem, rhs, solution_old, dt);
    hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
935
936
937
938
939
940
941
    printf("iter %i norm = %g\n", i, norm);
    if (norm < newton_atol)
      break;
    if (i == 0)
      norm0 = norm;
    else if(norm/norm0 < newton_rtol)
      break;
942
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
943
  hxtLinearSystemSolve(problem->linear_system, rhs, dx);
944
  clock_get(&toc);
945
  //break;
946
  totaltime += clock_diff(&tic, &toc);
947
948
949
950
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
951
  printf("total solve time : %g\n", totaltime);
952
  free(dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
953
  free(rhs);
954
  free(solution_old);
955
  return i;
956
957
}

Matthieu Constant's avatar
Matthieu Constant committed
958
959
960
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
  free(problem->node_volume);
  Mesh *mesh = problem->mesh;
961
  problem->node_volume = (double*)malloc(sizeof(double)*mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
962
963
964
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->node_volume[i] = 0;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed