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

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

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
36
37
#define LOCAL_MATRIX(i,j,a,b) local_matrix[((i)*n_fields*N_SF+N_SF*(a)+(j))*n_fields+(b)]

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
38
39
40
41
42
43
44
45
46
47
48
49
50
void fluid_problem_interpolate_rho_and_nu(const FluidProblem *problem, double a, double *rho, double *mu) {
  if (problem->n_fluids == 1){
    *rho = problem->rho[0];
    *mu = problem->mu[0];
  }
  else{
    //*rho = 1/(1/problem->rho[0]*a + 1/problem->rho[1]*(1-a));
    *rho = problem->rho[0]*a + problem->rho[1]*(1-a);
    //*mu = problem->mu[0]*a + problem->mu[1]*(1-a);
    *mu = (problem->mu[0]/problem->rho[0]*a + problem->mu[1]/problem->rho[1]*(1-a))*(*rho);
  }
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
inline static double mesh_dxidx(const Mesh *mesh, int iel, double dxidx[D][D]) {
  double dxdxi[D][D];
  const int *el = &mesh->elements[iel*N_N];
  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];
  return invDD(dxdxi, dxidx);
}

inline static double element_volume_from_detj(double detj) {
#if D == 2
    return detj/2;
#else        
    return detj/6;
#endif    
}

inline static double node_volume_from_detj(double detj) {
#if DIMENSION == 2
    return detj/6;
#else
    return detj/24;
#endif
}

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

78
#if D==2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
79
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
80
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
81
  double d = 2*sqrt(vol/M_PI);
82
83
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
84
  double Re_O_u = d*c/mu*rho;
85
86
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
Matthieu Constant's avatar
Matthieu Constant committed
87
  double f = pow(c,-(1.8-0.65*exp(-(1.5-log(Re_O_u*normu))*(1.5-log(Re_O_u*normu))/2.)));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
88
  return f*Cd_u*rho/2*d;
89
}
90
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
91
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
92
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
93
  double d = 2*cbrt(3./4.*vol/M_PI);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
94
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
95
96
97
98
99
100
  //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
101
  return f*Cd_u*rho/2*(M_PI*r*r);
102
103
}
#endif
104

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
105
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
106
  for (int i = 0; i < problem->n_particles*D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
107
108
109
110
    particle_force[i] = problem->particle_force[i];
  }
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
111
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
112
113
114
115
116
117
  double u[D], uold[D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    dp[iD] = dsol[P*D+iD];
  }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
118
  double mu, rho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
119
  fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
120
121
122
123
124
125
126
  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
127
  double g = problem->g;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
128
  if (problem->reduced_gravity){
129
130
131
    double rhop = mass/vol;
    g = g*(rhop-problem->rho[0])/rhop;
  }
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
132
  double gamma;
Matthieu Constant's avatar
Matthieu Constant committed
133
  
134
  if(problem->n_fluids == 2){
Matthieu Constant's avatar
adapt    
Matthieu Constant committed
135
136
137
138
139
140
141
    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);
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
142
143
144
145
146
147
148
149
150
151
152
153
#define USE_PATANKAR 1
#if USE_PATANKAR
  double fcoeff = mass/(mass+dt*gamma);
  for (int i = 0; i < D; ++i){
      problem->particle_force[ip*D+i] = fcoeff*(-gamma*Du[i]-vol*dp[i]);
    f[U+i] = -fcoeff*gamma*Du[i]+(gamma*dt/(dt*gamma+mass)-1)*vol*dp[i];
  }
  problem->particle_force[ip*D+1] += fcoeff*g*mass;
  f[U+1] += -fcoeff*gamma*dt*g;
  *dfdu = fcoeff*gamma/c;
  *dfddp = (gamma*dt/(dt*gamma+mass)-1)*vol;
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
154
155
  gamma = fmin(gamma,mass/dt);
  for (int j = 0; j < D; ++j){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
156
    problem->particle_force[ip*D+j] = (-gamma*Du[j]-vol*dp[j]);
157
    f[U+j] = -gamma*Du[j]-vol*dp[j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
158
  }
159
160
  *dfdu = gamma/c;
  *dfddp = -vol;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
161
#endif
162
163
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
164
static void node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
165
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
166

167
168
169
170
171
172
173
174
175
176
177
  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);
  for (int ip=0; ip < problem->n_particles; ++ip) {
    int iel = problem->particle_element_id[ip];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
178
    if(iel < 0){
179
      for (int d = 0; d < D; ++d)
180
181
        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
182
183
      continue;
    }
184
    double *xi = problem->particle_uvw + D*ip;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
185
186
    double phi[N_SF];
    shape_functions(xi,phi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
187
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
188
189
    double dxidx[D][D], dphi[N_SF][D];
    double detj = mesh_dxidx(mesh, iel, dxidx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
190
    grad_shape_functions(dxidx, dphi);
191
192
193
194
195
196
197
    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
198
199
      }
    }
200
    double c = 0;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
201
    double a = 0;
202
203
204
205
206
    double cold = 0;
    double dc[D] = {0};
    for (int i = 0; i < N_SF; ++i) {
      c += problem->porosity[el[i]]*phi[i];
      cold += problem->oldporosity[el[i]]*phi[i];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
207
208
209
      if (problem->n_fluids == 2) {
        a += problem->concentration[iel*N_N+i]*phi[i];
      }
210
211
212
213
214
215
216
217
218
219
220
      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
221
    }
222
    double f[D],dfdu,dfddp;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
223
    particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
224
    for (int iD = 0; iD < D; ++iD) {
225
      for (int iphi = 0; iphi < N_SF; ++iphi){
226
        local_vector[iphi+N_SF*(U+iD)] += phi[iphi]*f[iD];
227
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
228
    }
229
230
231
232
233
    for (int iD = 0; iD < D; ++iD){
      for (int iphi = 0; iphi < N_SF; ++iphi){
        for (int jphi = 0; jphi < N_SF; ++jphi){
          LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += phi[iphi]*phi[jphi]*dfdu;
          LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][iD]*dfddp;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
234
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
235
236
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
237
  }
238
239
240
  free(s);
  free(ds);
  free(sold);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
241
242
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
243
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, const double mu, const double dt, int eid, const double *data)
244
245
246
247
248
{
  const int vid = wbnd->vid;
  const int pid = wbnd->pid;
  const int cid = wbnd->cid;
  const int aid = wbnd->aid;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
249
  double p = s[P];
250
  if (wbnd->type == BND_WALL) {
251
    double u[D], du[D][D], dp[D], dun[D], dunt[D];
252
253
254
    for (int iD = 0; iD < D; ++iD) {
      u[iD] = s[U+iD];
      dp[iD] = ds[P*D+iD];
255
      dun[iD] = 0;
256
      dunt[iD] = 0;
257
      for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
Matthieu Constant's avatar
Matthieu Constant committed
258
    }
259
    double un = 0;
260
    double dcn = 0;
261
    double s_c = 0;
262
263
    for (int i = 0; i < D; ++i) {
      un += u[i]*n[i];
264
      dcn += dc[i]*n[i];
265
266
267
268
269
      s_c += vid<0?0:(u[i]-data[vid+i])*n[i];
      for (int j = 0; j < D; ++j) {
        dun[i] += du[i][j]*n[j];
        dunt[i] += du[j][i]*n[j];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
270
    }
271
272
    double h = problem->element_size[eid];
    f[P] = (pid<0 ? 0 : -(s[P]-data[pid])/h*problem->taup[eid]);
Matthieu Constant's avatar
Matthieu Constant committed
273

274
    for (int id = 0; id < D; ++id) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
275
276
      double sigma = (1+D)/(D*h)*mu*N_N;
      f[U+id] = (pid<0?0:data[pid]-p)*n[id] + (vid<0?0:sigma*((u[id]-data[vid+id])+s_c*n[id]));
277
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
278
  }
279
  else {
280
    double u[D], du[D][D], dp[D], dun[D], dunt[D];
281
282
283
    for (int iD = 0; iD < D; ++iD) {
      u[iD] = s[U+iD];
      dp[iD] = ds[P*D+iD];
284
285
      dun[iD] = 0;
      dunt[iD] = 0;
286
287
288
289
      for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
    }
    double un = 0;
    double unext = 0;
290
    double dcn = 0;
291
    const double *uext = vid < 0 ? &s[U] : &data[vid];
292
293
    double cext = cid<0?c:data[cid];
    double s_c = 0;
294
295
296
    for (int i = 0; i < D; ++i) {
      un += u[i]*n[i];
      unext += uext[i]*n[i];
297
298
      dcn += dc[i]*n[i];
      s_c += vid<0?0:(u[i]-data[vid+i])*n[i];
299
      for (int j = 0; j < D; ++j) {
300
301
        dun[i] += du[i][j]*n[j];
        dunt[i] += du[j][i]*n[j];
302
303
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
304

305
    double rho_ext =(problem->n_fluids == 1 ? rho : problem->rho[0]*data[aid]+problem->rho[1]*(1-data[aid]));
306
    double h = problem->element_size[eid];
Matthieu Constant's avatar
Matthieu Constant committed
307
    f[P] = unext;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
308
309
310
    if(vid<0){
      unext = 0;
    }
311
    for (int id = 0; id < D; ++id) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
312
      f[U+id] = (pid<0?0:data[pid]-p)*n[id] + (vid<0?0:(1+D)/(D*h)*mu*N_N*((u[id]-uext[id])+s_c*n[id])); 
313
    }
Matthieu Constant's avatar
Matthieu Constant committed
314
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
315
316
317
318
319
}

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
321
  problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
322
323
324
  problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
  problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
325
326
327
  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
328
329
  double maxtaa = 0;
  double mintaa = DBL_MAX;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
330
331
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    double normu = {0.};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
332
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
333
    double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
334
335
    double amax = 0.5;
    double amin = 0.5;
Matthieu Constant's avatar
Matthieu Constant committed
336
    double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
337
338
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
Matthieu Constant's avatar
Matthieu Constant committed
339
      c += problem->porosity[el[i]]/N_N;
340
      if(problem->n_fluids == 2){
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
341
342
343
        a += problem->concentration[iel*N_N+i];
        amax = fmax(amax,problem->concentration[iel*N_N+i]);
        amin = fmin(amin,problem->concentration[iel*N_N+i]);
Matthieu Constant's avatar
Matthieu Constant committed
344
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
345

346
      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
347
348
      normu = fmax(normu,sqrt(normup));
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
349
350
    double dxidx[DIMENSION][DIMENSION];
    double detj = mesh_dxidx(mesh,iel,dxidx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
351
#if DIMENSION == 2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
352
    const double h = 2*sqrt(detj/(2*M_PI));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
353
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
354
    const double h = 2*cbrt(detj/(8*M_PI));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
355
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
356
    problem->element_size[iel] = h;
Matthieu Constant's avatar
Matthieu Constant committed
357
    double rho, mu;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
358
359
360
    a /= N_N;
    a = fmin(1.,fmax(0.,a));
    fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
361
    double nu = mu/rho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
362
    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)))*3;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
363
    problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
364
365
  }
}
366

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
367
static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double cold, const double rho, const double mu, double dt, int iel) {
368
  double p = sol[P];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
369
370
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
371
372
373
  double u[D], uold[D], du[D][D], duold[D][D],dp[D];
  if(problem->stab_param != 0)
    taup = tauc = 0;
374
375
376
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
377
    dp[iD] = dsol[P*D+iD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
378
379
380
381
    for (int jD = 0; jD < D; ++jD) {
      du[iD][jD] = dsol[(U+iD)*D+jD];
      duold[iD][jD] = dsolold[(U+iD)*D+jD];
    }
382
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
383

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
384
385
386
  double divu = 0, divuold = 0;
  double tau[D][D],tauold[D][D];
  double utauold[D]={0};
Matthieu Constant's avatar
Matthieu Constant committed
387
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
388
389
  for (int i = 0; i < D; ++i) {
    divu += du[i][i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
390
    divuold += duold[i][i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
391
    for (int j = 0; j < D; ++j) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
392
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
393
      tauold[i][j] = duold[i][j]-uold[i]*dc[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
394
      utauold[i] += uold[j]*tauold[i][j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
395
    }
396
  }
Matthieu Constant's avatar
Matthieu Constant committed
397

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
398
  f0[P] = (c-cold)/dt;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
399
400
401
  double G[D] = {0};
  G[1] = -problem->g;
  double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
Matthieu Constant's avatar
Matthieu Constant committed
402

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
403
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
404
405
406
407
408
409
    f0[U+i] = 
      rho*(u[i]-uold[i])/dt 
      + dp[i] 
      + G[i]*rhoreduced*c
      + rho/c*(uold[i]*divuold + utauold[i])
      + problem->volume_drag*u[i]*mu;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
410
    double R = f0[U+i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
411
    for (int j = 0; j < D; ++j) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
412
      f1[(U+i)*D+j] = mu*(tau[i][j]+tau[j][i]) + R*uold[j]*taup;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
413
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
415
    f1[(U+i)*D+i] += ((c-cold)/dt+divu)*tauc*rho;
    f1[P*D+i] = -u[i]+ taup*R/rho + problem->stab_param*dp[i];
Matthieu Constant's avatar
Matthieu Constant committed
416
  }
417
418
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
419
static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
420
421
422
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
423
424
425
426
427
  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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
428
  double *dsold = malloc(sizeof(double)*n_fields*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
429
430
  double *f0 = malloc(sizeof(double)*n_fields);
  double *g0 = malloc(sizeof(double)*n_fields);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
431
  double *h0 = malloc(sizeof(double)*n_fields*N_LSF);
Matthieu Constant's avatar
Matthieu Constant committed
432
433
434
  double i_bnd = 0;
  double iv_bnd_up = 0;
  double iv_bnd_down = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
435
436
437
438
439
440
441
442
443
444
445
446
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
447
448
449
450
    int n_value = weak_boundary_n_values(wbnd);
    double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
    weak_boundary_values(mesh, bnd, wbnd, data);

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
451
452
453
454
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
455
456
457
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
458
      const int *el = &mesh->elements[iel*N_N];
459
      for (int i = 0; i < D; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
460
        nodes[i] = el[cl[i]];
461
462
        x[i] = &mesh->x[nodes[i]*3];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
463
464
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
465
466
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
467
      grad_shape_functions(dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
468
469
470
471
472
473
474
      double dc[D] = {0};
      for (int j = 0; j < n_fields*D; ++j) {
        ds[j] = 0;
      }
      for (int i = 0; i < N_SF; ++i) {
        for (int j = 0; j < n_fields; ++j) {
          double dof = solution[el[i]*n_fields+j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
475
476
          double dofold = solution_old[el[i]*n_fields+j];
          for (int k = 0; k < D; ++k){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
477
            ds[j*D+k] += dphi[i][k]*dof;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
478
479
            dsold[j*D+k] += dphi[i][k]*dofold;
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
480
481
482
483
484
        }
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
485
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
486
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
487
488
489
        double dp[D]={0};
        double *local_vector = &all_local_vector[local_size*iel];
        double *local_matrix = &all_local_matrix[local_size*local_size*iel];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
490
491
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
492
493
494
495
496
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
          sold[i] = 0;
        }
        double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
497
        double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
498
499
        for (int i = 0; i < DIMENSION; ++i) {
          c += problem->porosity[nodes[i]]*phi[i];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
500
          if (problem->n_fluids == 2) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
501
            a += problem->concentration[iel*N_N+cl[i]]*phi[i];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
502
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
503
          for (int j = 0; j < n_fields; ++j) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
504
            double dof = solution[nodes[i]*n_fields+j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
505
            double dofold = solution_old[nodes[i]*n_fields+j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
506
            s[j] += phi[i]*dof;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
507
            sold[j] += phi[i]*dofold;
Matthieu Constant's avatar
Matthieu Constant committed
508
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
509
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
510
511
        double rho, mu;
        fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
512
        const double jw = LQW[iqp]*detbnd;
Matthieu Constant's avatar
Matthieu Constant committed
513
514
515
516
517
518
519
520
521
522
523
524
        if (wbnd->type != BND_WALL){
          for (int i = 0; i < D; ++i){
            if (wbnd->vid<0) {
              i_bnd -= a*s[U+i]*n[i]*jw*dt;
              iv_bnd_up -= s[U+i]*n[i]*jw;
            }
            else {
              iv_bnd_down -= dataqp[wbnd->vid+i]*n[i]*jw;
              i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
            }
          }
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
525
        f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
526
        for (int ifield = 0; ifield < n_fields; ++ifield) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
527
528
          for (int iphi = 0; iphi < D; ++iphi){
            local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
529
530
531
532
533
          }
        }
        for (int jfield = 0; jfield < n_fields; ++jfield) {
          double store = s[jfield];
          s[jfield] += deps;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
534
          f_boundary(wbnd,problem,n,g0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
535
536
537
538
          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
539
            f_boundary(wbnd,problem,n,h0+jd*n_fields,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
540
541
542
            ds[jfield*D+jd] = store;
          }
          for (int ifield = 0; ifield < n_fields; ++ifield){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
543
544
            for (int iphi = 0; iphi < N_LSF; ++iphi){
              for (int jphi = 0; jphi < N_LSF; ++jphi){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
545
                double df00 = (g0[ifield]-f0[ifield])/deps;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
546
547
548
549
550
                LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*phi[jphi]*df00;
              }
              for (int jphi = 0; jphi < N_SF; ++jphi){
                double m = 0;
                for (int jd = 0; jd < N_LSF; ++jd) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
551
552
553
                  double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
                  m += jw*phi[iphi]*dphi[jphi][jd]*df01;
                }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
554
                LOCAL_MATRIX(cl[iphi],jphi,ifield,jfield) += m;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
555
556
557
558
559
              }
            }
          }
        }
      }
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
560
    }
561
    free(data);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
562
  }
563
564
565
  free(s);
  free(ds);
  free(sold);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
566
  free(dsold);
567
568
569
  free(f0);
  free(g0);
  free(h0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
570
}
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
571

Matthieu Constant's avatar
Matthieu Constant committed
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
double fluid_problem_a_integ_volume(FluidProblem *problem)
{
  const Mesh *mesh = problem->mesh;
  double I_a = 0;
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D];
    const double detj = mesh_dxidx(mesh,iel,dxidx);
    for (int i = 0; i < N_N; ++i)
      for (int j = 0; j< N_N; ++j) 
        I_a += detj * problem->porosity[el[i]]*problem->concentration[iel*N_N+j]*mass_matrix[i][j];
  }
  return I_a;
}

587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
double fluid_problem_volume_flux(FluidProblem *problem, const char *tagname)
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  double Q = 0;
  double wtotal = 0;
  int n_fields = D+1;
  int found = 0;
  for (int i = 0; i < N_LQP; ++i) wtotal += LQW[i];
  for (int ibnd = 0; ibnd < mesh->n_boundaries; ++ibnd) {
    if (strcmp(mesh->boundary_names[ibnd],tagname) != 0)
      continue;
    found = 1;
    MeshBoundary *bnd = &problem->boundaries[ibnd];
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      double vnmean = 0;
      for (int i = 0; i < D; ++i) {
        for (int iD = 0; iD < D; ++iD) {
          vnmean += n[iD]*problem->solution[nodes[i]*n_fields+U+iD];
        }
      }
      Q += vnmean*detbnd*wtotal/D;
    }
  }
  if(found == 0) {
    printf("boundary '%s' not found\n", tagname);
    exit(1);
  }
  printf("Q = %g wtotal = %g\n", Q,wtotal);
  return Q;
}

Matthieu Constant's avatar
Matthieu Constant committed
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
double fluid_problem_a_integ_bnd(FluidProblem *problem, double dt)
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  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 i_bnd = 0;
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    int n_value = weak_boundary_n_values(wbnd);
    double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
    weak_boundary_values(mesh, bnd, wbnd, data);

    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      grad_shape_functions(dxidx, dphi);
      
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
        }
        double a = 0;
        for (int i = 0; i < DIMENSION; ++i) {
          if (problem->n_fluids == 2) {
            a += problem->concentration[iel*N_N+cl[i]]*phi[i];
          }
          for (int j = 0; j < n_fields; ++j) {
            double dof = solution[nodes[i]*n_fields+j];
            s[j] += phi[i]*dof;
          }
        }
        const double jw = LQW[iqp]*detbnd;
        if (wbnd->type != BND_WALL){
          for (int i = 0; i < D; ++i){
            if (wbnd->vid<0) {
              i_bnd -= a*s[U+i]*n[i]*jw*dt;
            }
            else {
              i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
            }
          }
        }
      }
    }
  }
  return i_bnd;
}

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

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
844
static void fluid_problem_surface_tension_bnd(FluidProblem *problem, double *a_cg, double *all_local_vector){
845
  const Mesh *mesh = problem->mesh;
846
  double sigma = problem->sigma;
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
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
910
911
  const int n_fields = fluid_problem_n_fields(problem);
  const size_t local_size = N_SF*n_fields;
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    int n_value = weak_boundary_n_values(wbnd);
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[D];
      double *x[D];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[elbnd[icl][i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      grad_shape_functions(dxidx, dphi);
      double da[D] = {0};
      double nda = 0;
      for (int k=0; k<D; ++k){
        for (int i = 0; i < N_SF; ++i) {
          da[k] += a_cg[el[i]]*dphi[i][k];
        }
        nda += da[k]*da[k];
      }
      nda = fmax(sqrt(nda),1e-8);
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *local_vector = &all_local_vector[local_size*iel];
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
        const double jw = LQW[iqp]*detbnd;
        double a = 0;
        double c = 0;
        for (int i = 0; i < D; ++i){
          a += a_cg[nodes[i]]*phi[i];
          c += problem->porosity[nodes[i]]*phi[i];
        }
        for (int iphi = 0; iphi < D; ++iphi) {
          for (int iD = 0; iD < D; ++iD) {
            for (int jD = 0; jD < D; ++jD){
              local_vector[cl[iphi]+N_SF*iD] += phi[iphi]*jw*c*da[iD]*da[jD]*n[jD]/nda*sigma;
            }
            local_vector[cl[iphi]+N_SF*iD] -= phi[iphi]*jw*c*nda*n[iD]*sigma;
          }
        }
      }
    }
  }
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
912
static void fluid_problem_surface_tension(FluidProblem *problem, double *a_cg, double *all_local_vector){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
913
914
  if (problem->n_fluids == 1) return;
  const Mesh *mesh = problem->mesh;
915
  double sigma = problem->sigma;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
916
917
918
  int n_fields = fluid_problem_n_fields(problem);
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
919
920
921
    double dxidx[D][D], dphi[N_N][D];
    const double det = mesh_dxidx(mesh,iel,dxidx);
    const double vol = element_volume_from_detj(det);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
922
923
924
    grad_shape_functions(dxidx, dphi);
    double da[D] = {0};
    double nda = 0;
Matthieu Constant's avatar
Matthieu Constant committed
925
926
927
928
    double c = 0;
    for (int i = 0; i < N_N; ++i){
      c += problem->porosity[el[i]]/N_N;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
929
930
    for (int k=0; k<D; ++k){
      for (int i = 0; i < N_SF; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
931
932
        //da[k] += a_cg[el[i]]*dphi[i][k];
        da[k] += problem->concentration[iel*N_N+i]*dphi[i][k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
933
934
935
936
      }
      nda += da[k]*da[k];
    }
    nda = fmax(sqrt(nda),1e-8);
937
    double *local_vector = &all_local_vector[iel*n_fields*N_N];
938
939
940
    for (int iphi = 0; iphi < N_N; ++iphi) {
      for (int id = 0; id < D; ++id) {
        for (int jd = 0; jd < D; ++jd) {
Matthieu Constant's avatar
Matthieu Constant committed
941
          local_vector[(U+id)*N_N+iphi] -= c*vol*da[id]*da[jd]*dphi[iphi][jd]/nda*sigma;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
942
        }
Matthieu Constant's avatar
Matthieu Constant committed
943
        local_vector[(U+id)*N_N+iphi] += c*vol*nda*dphi[iphi][id]*sigma;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
944
945
946
      }
    }
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
}


static void fluid_problem_dg_to_cg(FluidProblem *problem, const double *dg, double *cg) {
  const Mesh *mesh = problem->mesh;
  for (int i = 0; i< mesh->n_nodes; ++i) {
    cg[i] = 0;
  }
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D];
    const double det = mesh_dxidx(mesh,iel,dxidx);
    const double vol = node_volume_from_detj(det);
    for (int i = 0; i< N_N; ++i) {
      cg[el[i]] += vol * dg[iel*N_N+i];
    }
  }
  for (int i = 0; i< mesh->n_nodes; ++i) {
    cg[i] /= problem->node_volume[i];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
967
968
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
969
static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
970
9