fluid_problem.c 52 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
Jonathan Lambrechts committed
31
#include "mbxml.h"
32

33
#define D DIMENSION
34

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

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

43
44
static double LQP[] = {0.21132486540518713, 0.7886751345948129};
static double LQW[] = {0.5,0.5};
45
46
static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {1./6,1./6,1./6};
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
static void shape_functions(double *uvw, double *sf)
{
  sf[0] = 1-uvw[0]-uvw[1];
  sf[1] = uvw[0];
  sf[2] = uvw[1];
}
static void grad_shape_functions(const double dxidx[2][2], double gsf[3][2])
{
  for (int i = 0;  i < 2; ++i) {
    gsf[0][i] = -dxidx[0][i]-dxidx[1][i];
    gsf[1][i] = dxidx[0][i];
    gsf[2][i] = dxidx[1][i];
  }
}
static double detDD(const double m[2][2])
{
  return m[0][0]*m[1][1] - m[0][1]*m[1][0];
}
static double invDD(const double m[2][2], double inv[2][2])
{
  double d = detDD(m);
  inv[0][0] = m[1][1]/d;
  inv[0][1] = -m[0][1]/d;
  inv[1][0] = -m[1][0]/d;
  inv[1][1] = m[0][0]/d;
  return d;
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
74
#else
75
76
77

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
78
#define N_QP 5
79
#define N_LQP 3
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
81
82
83
84
85
86
static double QP[][3] = {
  {0.25, 0.25, 0.25},
  {0.5, 1./6, 1./6},
  {1./6, 0.5, 1./6},
  {1./6, 1./6, 0.5},
  {1./6, 1./6, 1./6}
};
87
static double LQP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
88
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
89
static double LQW[] = {1./6,1./6,1./6};
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
static void shape_functions(double *uvw, double *sf)
{
  sf[0] = 1-uvw[0]-uvw[1]-uvw[2];
  sf[1] = uvw[0];
  sf[2] = uvw[1];
  sf[3] = uvw[2];
}
static void grad_shape_functions(const double dxidx[3][3], double gsf[4][3])
{
  for (int i = 0; i < 3; ++i) {
    gsf[0][i] = -dxidx[0][i]-dxidx[1][i]-dxidx[2][i];
    gsf[1][i] = dxidx[0][i];
    gsf[2][i] = dxidx[1][i];
    gsf[3][i] = dxidx[2][i];
  }
};
static double detDD(const double m[3][3])
{
  return m[0][0]*(m[1][1]*m[2][2]-m[2][1]*m[1][2])
        -m[1][0]*(m[0][1]*m[2][2]-m[2][1]*m[0][2])
        +m[2][0]*(m[0][1]*m[1][2]-m[1][1]*m[0][2]);
}
static double invDD(const double m[3][3], double inv[3][3]) {
  double d = detDD(m);
  inv[0][0] = (m[1][1]*m[2][2]-m[2][1]*m[1][2])/d;
  inv[0][1] = -(m[0][1]*m[2][2]-m[2][1]*m[0][2])/d;
  inv[0][2] = (m[0][1]*m[1][2]-m[1][1]*m[0][2])/d;
  inv[1][0] = -(m[1][0]*m[2][2]-m[2][0]*m[1][2])/d;
  inv[1][1] = (m[0][0]*m[2][2]-m[2][0]*m[0][2])/d;
  inv[1][2] = -(m[0][0]*m[1][2]-m[1][0]*m[0][2])/d;
  inv[2][0] = (m[1][0]*m[2][1]-m[2][0]*m[1][1])/d;
  inv[2][1] = -(m[0][0]*m[2][1]-m[2][0]*m[0][1])/d;
  inv[2][2] = (m[0][0]*m[1][1]-m[1][0]*m[0][1])/d;
  return d;
};

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
407
  for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
408
    f[U+id] = p*n[id]+ (un> 0 ? rho*u[id]*un/c : 1);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
409
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
410
411
}

Matthieu Constant's avatar
Matthieu Constant committed
412
static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc,const double dt,int eid,double v)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
413
414
{
  double p = s[P];
415
416
417
418
419
420
  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
421
  double un = 0;
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
422
  double a;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
423
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
424
    un += u[i]*n[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
425
  }
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
426
427
428
429
430
431
432
433
434
  double rho;
  if (problem->n_fluids == 1){
    rho = problem->rho[0];
  }
  else{
    a = s[A];
    a = fmin(1.,fmax(0.,a));
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
435
436

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

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

  for (int id = 0; id < D; ++id) {
    f[U+id] = p*n[id]+ (un> 0 ? rho*u[id]*un/c : 0);
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
446
447
448
449
450
}

static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
  const Mesh *mesh = problem->mesh;
451
452
453
454
455
456

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

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

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

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

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

567
  if(problem->n_fluids == 2){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
568
    f0[A] = (a*c-aold*cold)/dt;
Matthieu Constant's avatar
Matthieu Constant committed
569
570
  }
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
571
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
572
    double dudt = (u[i]-uold[i])/dt;
573
574
    double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==1 ? -(reduced_gravity ? 0 : 1)*c*problem->g : 0.);
    f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -(reduced_gravity ? 0 : 1)*c*rho*problem->g : 0.);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
575
    for (int j = 0; j < D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
576
      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
577
    }
Matthieu Constant's avatar
Matthieu Constant committed
578
579
    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]);
580
    if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
581
582
      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
583
    }
584
  }
585
  /*if(problem->n_fluids == 2){
Matthieu Constant's avatar
Matthieu Constant committed
586
    f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
587
  }*/
588
589
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
590
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
591
592
593
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
594
595
596
597
598
599
600
601
  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
602
  int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
603
  f_cb *f;
604
605
606
607
608
#if DIMENSION == 2
  int elbnd[3][2] = {{0,1},{1,2},{2,0}};
#else
  int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
609
  for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
Matthieu Constant's avatar
Matthieu Constant committed
610
611
    for (int i =0; i < mesh->n_nodes; ++i)
      bnd_node_local_id[i] = -1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
612
    if(mesh->phys_dimension[iphys] != D-1) continue;
Matthieu Constant's avatar
Matthieu Constant committed
613
614
615
616
617
618
619
620
621
    double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
    double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
      bnd_node_local_id[j] = i;
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
    }
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
622
623
624
625
626
    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;
Matthieu Constant's avatar
Matthieu Constant committed
627
        bnd->apply(mesh->phys_n_nodes[iphys], x, v);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
628
629
630
631
        break;
      }
    }
    if (f == NULL) {
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
632
      f = f_wall;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
633
634
      //printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
      //exit(1);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
635
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
636
637
638
639

    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];
640
641
642
      const int iebnd = bnd[1];
      int nodes[D];
      double *x[D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
643
      const unsigned int *el = &mesh->elements[iel*N_N];
644
645
646
647
      for (int i = 0; i < D; ++i){
        nodes[i] = el[elbnd[iebnd][i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
648
649
650
651
652
653
654
      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);

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
767

Matthieu Constant's avatar
Matthieu Constant committed
768
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
769
770
{
  const Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
771
  const double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
772
773
774
775

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

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

Matthieu Constant's avatar
Matthieu Constant committed
883
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
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
{
  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)
910
      printf("Boundary tag \"%s\" not found.\n", bnd->tag);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
911
912
913
914
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
    }
  }
915
  node_force_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity);
Matthieu Constant's avatar
Matthieu Constant committed
916
  compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
Matthieu Constant's avatar
Matthieu Constant committed
917
  fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix,reduced_gravity,stab_param);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed