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

24
#include "tools.h"
25
26
#include <stdlib.h>
#include <stdio.h>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
27
#include <float.h>
28
29
#include <math.h>
#include "fluid_problem.h"
30
#include "mesh_find.h"
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)]

38
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);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
39
40
41
42
43
44
45
46
47
48
49
50
51
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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
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
77
int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}
78

79

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
80
void fluid_problem_node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
81
82
83
84
{
  const Mesh *mesh = problem->mesh;
  const double *porosity = problem->porosity;
  const double *solution = problem->solution;
85
  int drag_in_stab = problem->drag_in_stab;
86
87
88
89
90
91
92
  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
93
    if(iel < 0){
94
      for (int d = 0; d < D; ++d)
95
        problem->particle_force[ip*D+d] = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
96
97
      continue;
    }
98
    double *xi = problem->particle_uvw + D*ip;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
99
100
    double phi[N_SF];
    shape_functions(xi,phi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
102
103
    double dxidx[D][D], dphi[N_SF][D];
    double detj = mesh_dxidx(mesh, iel, dxidx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
104
    grad_shape_functions(dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
105
106
    double *local_vector = all_local_vector ? &all_local_vector[local_size*iel] : NULL;
    double *local_matrix = all_local_matrix ? &all_local_matrix[local_size*local_size*iel] : NULL;
107
108
109
110
111
    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
112
113
      }
    }
114
    double c = 0;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
115
    double a = 0;
116
117
118
119
120
    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
121
122
123
      if (problem->n_fluids == 2) {
        a += problem->concentration[iel*N_N+i]*phi[i];
      }
124
125
126
127
128
129
130
131
132
133
134
      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
135
    }
136
    double f[D],dfdu,dfddp;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
137
    particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
138
139
    if (!local_vector)
      continue;
140
141
    double rho,nu;
    fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
142
    double pspg = drag_in_stab?problem->taup[iel]/rho:0;
143
    for (int iD = 0; iD < D; ++iD) {
144
      for (int iphi = 0; iphi < N_SF; ++iphi){
145
        local_vector[iphi+N_SF*(U+iD)] += phi[iphi]*f[iD];
146
147
        local_vector[iphi+N_SF*P] += pspg*dphi[iphi][iD]*f[iD];
        for (int jD = 0; jD < D; ++jD) {
148
          double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
149
150
          local_vector[iphi+N_SF*(U+iD)] += supg*dphi[iphi][jD]*f[iD];
        }
151
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
152
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
153
154
    if (!local_matrix)
      continue;
155
156
157
158
159
    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;
160
161
162
          LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu;
          LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][iD]*dfddp;
          for (int jD = 0; jD < D; ++jD) {
163
            double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
164
165
166
            LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += supg*dphi[iphi][jD]*phi[jphi]*dfdu;
            LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][jD]*dphi[jphi][iD]*dfddp;
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
167
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
168
169
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
170
  }
171
172
173
  free(s);
  free(ds);
  free(sold);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
174
175
}

176
static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double *f01)
177
178
179
180
181
{
  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
182
  double p = s[P];
183
  const int n_fields = fluid_problem_n_fields(problem);
184
  double u[D], dp[D];
185
  double s_c = 0;
186
  double unold = 0;
187
  double unext = 0;
188
189
190
191
  double dcn = 0;
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
192
    unold += sold[U+iD]*n[iD];
193
    unext += (vid<0?s[U+iD]:data[vid+iD])*n[iD];
194
195
    dcn += dc[iD]*n[iD];
    s_c += vid<0?0:(u[iD]-data[vid+iD])*n[iD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
196
  }
197
  double h = problem->element_size[eid];
198
199
200
  if (wbnd->type == BND_WALL && pid >= 0) {
      f0[P] = -(s[P]-data[pid])/h*problem->taup[eid];
      f00[P*n_fields+P] += -1/h*problem->taup[eid];
201
  }
202
  if(wbnd->type == BND_OPEN) {
203
    for (int i = 0; i < D; ++i) {
204
      f00[P*n_fields+U+i] += vid<0?n[i]:0;
205
    }
206
207
208
    f0[P] = unext;
  }
  for (int id = 0; id < D; ++id) {
209
    f0[U+id] = c*(pid<0?0:data[pid]-p)*n[id];
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
210
    f00[(U+id)*n_fields+P] += c*(pid<0?0:-n[id]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
211
212
  }
  double c_du_o_c[D][D];
Matthieu Constant's avatar
Matthieu Constant committed
213
  double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
214
215
216
217
218
219
220
221
222
223
224
225
226
  for (int iD = 0; iD < D; ++iD) {
    for (int jD = 0; jD < D; ++jD)
      c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
  }
  for (int id = 0; id < D; ++id) {
    f0[U+id] += vid<0?0:sigma*(u[id]-data[vid+id]+s_c*n[id]);
    f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
    for (int jd = 0; jd <D; ++jd) {
      f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
      f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
      f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
      f01[(U+id)*n_fields*D+(U+id)*D+jd] -= mu*n[jd];
      f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
227
228
229
230
231
232
233
234
235
236
237
    }
  }
  if (problem->advection) {
    double unbnd = unold;
    if (wbnd->type == BND_WALL) unbnd = unold/2;
    else if(vid>0) unbnd = (unold+unext)/2;
    if (unbnd<0) {
      for (int id = 0; id < D; ++id) {
        f0[U+id] += (unbnd*(vid<0?0:data[vid+id])-unold*u[id])*rho/c;
        f00[(U+id)*n_fields+U+id] -= unold*rho/c;
      }
238
    }
Matthieu Constant's avatar
Matthieu Constant committed
239
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
240
241
242
}

static double pow2(double a) {return a*a;}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
243
void fluid_problem_compute_stab_parameters(FluidProblem *problem, double dt) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
244
  const Mesh *mesh = problem->mesh;
245

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
246
  problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
247
248
249
  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
250
251
252
  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
253
254
  double maxtaa = 0;
  double mintaa = DBL_MAX;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
255
256
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    double normu = {0.};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
257
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
258
    double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
259
260
    double amax = 0.5;
    double amin = 0.5;
Matthieu Constant's avatar
Matthieu Constant committed
261
    double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
262
263
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
Matthieu Constant's avatar
Matthieu Constant committed
264
      c += problem->porosity[el[i]]/N_N;
265
      if(problem->n_fluids == 2){
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
266
267
268
        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
269
      }
270
      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
271
272
      normu = fmax(normu,sqrt(normup));
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
273
274
    double dxidx[DIMENSION][DIMENSION];
    double detj = mesh_dxidx(mesh,iel,dxidx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
275
#if DIMENSION == 2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
276
    const double h = 2*sqrt(detj/(2*M_PI));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
277
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
278
    const double h = 2*cbrt(detj/(8*M_PI));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
279
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
280
    problem->element_size[iel] = h;
Matthieu Constant's avatar
Matthieu Constant committed
281
    double rho, mu;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
282
283
284
    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
285
    double nu = mu/rho;
286
287
288
289
290
    problem->taup[iel] = 1/sqrt(
        (problem->temporal?pow2(2/dt):0.)
        +(problem->advection ?pow2(2*normu/h):0.)
        +pow2(4*nu/pow2(h)));
    problem->tauc[iel] = problem->advection?h*normu*fmin(h*normu/(6*nu),0.5):0.;
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
291
292
  }
}
293

294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
#if D==2
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
{
  double d = 2*sqrt(vol/M_PI);
  double normu = hypot(u[0],u[1]);
  //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-0.65*exp(-(1.5-log(Re_O_u*normu))*(1.5-log(Re_O_u*normu))/2.)));
  return f*Cd_u*rho/2*d;
}
#else
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
{
  double d = 2*cbrt(3./4.*vol/M_PI);
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
  //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.;
  return f*Cd_u*rho/2*(M_PI*r*r);
}
#endif

void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
  for (int i = 0; i < problem->n_particles*D; ++i) {
    particle_force[i] = problem->particle_force[i];
  }
}

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) {
  const double *contact = &problem->particle_contact[ip*D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
329
  double u[D], uold[D], dp[D];
330
331
332
333
334
335
336
337
338
339
340
341
342
343
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    dp[iD] = dsol[P*D+iD];
  }
  double mu, rho;
  fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
  double Due[D];
  const double *up = &problem->particle_velocity[ip*D];
  for (int j = 0; j < D; ++j) {
    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
344
345
346
347
  double rhop = mass/vol;
  double massreduced = mass;
  if(problem->reduced_gravity) {
    massreduced = (rhop-problem->rho[0])*problem->particle_volume[ip];
348
349
350
351
352
353
354
355
356
  }
  double gamma;
  
  if(problem->n_fluids == 2){
    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{
357
    gamma = particle_drag_coeff(Due,mu,rho,vol,c);
358
  }
Matthieu Constant's avatar
Matthieu Constant committed
359
360
  gamma = gamma*problem->drag_coeff_factor;
  double gammastar = gamma/(1+dt/mass*gamma);
361
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
362
    double upstar = up[i]+dt/mass*(contact[i]+massreduced*problem->g[i]-vol*dp[i]);
363
    double drag = gammastar*(upstar-u[i]/c);
364
    problem->particle_force[ip*D+i] = -drag-vol*dp[i];
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
365
    f[U+i] = -drag;//-vol*dp[i];
366
367
  }
  *dfdu = gammastar/c;
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
368
  *dfddp = gammastar*dt/mass*vol;//-vol;
369
370
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
371
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double c, const double *dc, const double *bf, const double cold, const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
372
  double p = sol[P];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
373
374
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
375
  const int n_fields = fluid_problem_n_fields(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
376
377
378
  double u[D], uold[D], du[D][D], duold[D][D],dp[D];
  if(problem->stab_param != 0)
    taup = tauc = 0;
379
  double divu = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
380
  double nold = 0;
381
382
383
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
384
    nold += uold[iD]*uold[iD];
385
    dp[iD] = dsol[P*D+iD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
386
387
388
389
    for (int jD = 0; jD < D; ++jD) {
      du[iD][jD] = dsol[(U+iD)*D+jD];
      duold[iD][jD] = dsolold[(U+iD)*D+jD];
    }
390
    divu += du[iD][iD];
391
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
392
  nold = sqrt(nold);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
393
394
  f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
  //f00[P*n_fields+P] = 1/dt*0.1;
Matthieu Constant's avatar
Matthieu Constant committed
395
  double *g = problem->g;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
396
  double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
397
  double drag = problem->volume_drag + problem->quadratic_drag*nold;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
398
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
399
    f0[U+i] = 
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
400
      + c*dp[i]
Matthieu Constant's avatar
Matthieu Constant committed
401
      - g[i]*rhoreduced*c
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
402
      - bf[i]*c
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
403
404
405
      + drag*u[i];
      //+ 5.3e5*u[i];
    f00[(U+i)*n_fields+U+i] = drag;//5.3e5;
406
407
408
409
    if(problem->temporal) {
      f0[U+i] += rho*(u[i]-uold[i])/dt;
      f00[(U+i)*n_fields+U+i] += rho/dt;
    }
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
410
    f01[(U+i)*n_fields*D+P*D+i] = c;
411
412
413
414
415
416
417
418
419
    if (problem->advection) {
      // advection :
      // dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
      //              = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
      for (int j = 0; j < D; ++j) {
        f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
        f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
        f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j];
      }
420
421
    }
    for (int j = 0; j < D; ++j) {
422
      f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
423
424
425
426
427
428
      f10[((U+i)*D+j)*n_fields+U+j] += -mu*dc[i]/c;
      f10[((U+i)*D+j)*n_fields+U+i] += -mu*dc[j]/c;
      f11[((U+i)*D+j)*n_fields*D+(U+j)*D+i] += mu; 
      f11[((U+i)*D+j)*n_fields*D+(U+i)*D+j] += mu;
      // SUPG
      double supg = uold[j]*taup;
429
      f1[(U+i)*D+j] += supg*f0[U+i];
430
431
432
433
434
      f10[((U+i)*D+j)*n_fields+U+i] += supg*f00[(U+i)*n_fields+U+i];
      f11[((U+i)*D+j)*(n_fields*D)+P*D+i] += supg*f01[(U+i)*n_fields*D+P*D+i]; 
      for (int k = 0; k < D; ++k)
        f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
    }
435
436
    double lsic = tauc*rho;
    f1[(U+i)*D+i] += ((c-cold)/dt+divu)*lsic;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
437
    for (int j = 0; j < D; ++j) {
438
      f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
439
    }
440
441
442
    f1[P*D+i] = -u[i];
    f10[(P*D+i)*n_fields+U+i] += -1;
    // PSPG
Jonathan Lambrechts's avatar
pp    
Jonathan Lambrechts committed
443
    double pspg = taup/rho;
444
445
    f1[P*D+i] += pspg*(f0[U+i]+(problem->drag_in_stab?0:(1-c)*dp[i])) + problem->stab_param*dp[i];
    f11[(P*D+i)*(n_fields*D)+P*D+i] += pspg*(f01[(U+i)*n_fields*D+P*D+i]+(problem->drag_in_stab?0:1-c))+problem->stab_param;
Jonathan Lambrechts's avatar
pp    
Jonathan Lambrechts committed
446
    f10[(P*D+i)*n_fields+U+i] += pspg*f00[(U+i)*n_fields+(U+i)];
447
    for (int j = 0; j < D; ++j) {
Jonathan Lambrechts's avatar
pp    
Jonathan Lambrechts committed
448
      f11[(P*D+i)*n_fields*D+(U+i)*D+j] = pspg*f01[(U+i)*n_fields*D+(U+i)*D+j];
449
    }
Matthieu Constant's avatar
Matthieu Constant committed
450
  }
451
452
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
453
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
454
455
456
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
457
458
459
460
461
  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
462
  double *dsold = malloc(sizeof(double)*n_fields*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
463
  double *f0 = malloc(sizeof(double)*n_fields);
464
  double *f00 = malloc(sizeof(double)*n_fields*n_fields);
Jonathan Lambrechts's avatar
tmp    
Jonathan Lambrechts committed
465
  double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
Matthieu Constant's avatar
Matthieu Constant committed
466
467
468
  double i_bnd = 0;
  double iv_bnd_up = 0;
  double iv_bnd_down = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
469
470
471
472
473
474
475
476
477
478
479
480
  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];
481
482
483
484
    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
485
486
487
488
    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
489
490
491
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
492
      const int *el = &mesh->elements[iel*N_N];
493
      for (int i = 0; i < D; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
494
        nodes[i] = el[cl[i]];
495
496
        x[i] = &mesh->x[nodes[i]*3];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
497
498
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
499
500
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
501
      grad_shape_functions(dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
502
503
504
505
506
507
508
      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
509
510
          double dofold = solution_old[el[i]*n_fields+j];
          for (int k = 0; k < D; ++k){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
511
            ds[j*D+k] += dphi[i][k]*dof;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
512
513
            dsold[j*D+k] += dphi[i][k]*dofold;
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
514
515
516
517
518
        }
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
519
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
520
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
521
522
523
        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
524
525
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
526
527
528
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
          sold[i] = 0;
529
          f0[i] = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
530
        }
531
532
533
        for (int i = 0; i < n_fields*n_fields; ++i) {
          f00[i] = 0;
        }
Jonathan Lambrechts's avatar
tmp    
Jonathan Lambrechts committed
534
535
536
        for (int i = 0; i < n_fields*n_fields*D; ++i) {
          f01[i] = 0;
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
537
        double c = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
538
        double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
539
540
        for (int i = 0; i < DIMENSION; ++i) {
          c += problem->porosity[nodes[i]]*phi[i];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
541
          if (problem->n_fluids == 2) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
542
            a += problem->concentration[iel*N_N+cl[i]]*phi[i];
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
543
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
544
          for (int j = 0; j < n_fields; ++j) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
545
            double dof = solution[nodes[i]*n_fields+j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
546
            double dofold = solution_old[nodes[i]*n_fields+j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
547
            s[j] += phi[i]*dof;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
548
            sold[j] += phi[i]*dofold;
Matthieu Constant's avatar
Matthieu Constant committed
549
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
550
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
551
552
        double rho, mu;
        fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
553
        const double jw = LQW[iqp]*detbnd;
Matthieu Constant's avatar
Matthieu Constant committed
554
555
556
557
558
559
560
561
562
563
564
565
        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
tmp    
Jonathan Lambrechts committed
566
        f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
567
        for (int ifield = 0; ifield < n_fields; ++ifield) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
568
          for (int iphi = 0; iphi < D; ++iphi){
569
570
571
            if (ifield<D){
              problem->boundary_force[bndid*D + ifield] += phi[iphi]*f0[ifield]*jw;
            }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
572
            local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
573
574
575
576
          }
        }
        for (int jfield = 0; jfield < n_fields; ++jfield) {
          for (int ifield = 0; ifield < n_fields; ++ifield){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
577
578
            for (int iphi = 0; iphi < N_LSF; ++iphi){
              for (int jphi = 0; jphi < N_LSF; ++jphi){
Jonathan Lambrechts's avatar
tmp    
Jonathan Lambrechts committed
579
580
581
582
583
584
585
586
587
                double d = phi[jphi]*f00[ifield*n_fields+jfield];
                LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*d;
              }
              for (int jphi = 0; jphi < N_SF; ++jphi){
                double d = 0;
                for (int jD = 0; jD < D; ++jD) {
                  d += dphi[jphi][jD]*f01[ifield*n_fields*D+jfield*D+jD];
                }
                LOCAL_MATRIX(cl[iphi],jphi,ifield,jfield) += jw*phi[iphi]*d;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
588
589
590
591
592
              }
            }
          }
        }
      }
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
593
    }
594
    free(data);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
595
  }
596
597
598
  free(s);
  free(ds);
  free(sold);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
599
  free(dsold);
600
  free(f0);
601
  free(f00);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
602
}
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
603

Matthieu Constant's avatar
Matthieu Constant committed
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
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;
}

619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
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
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
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
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
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
742
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
743
744
{
  const Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
745
  const double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
746
747
748
749

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
750
751
752
  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
753
  double *dsold = malloc(sizeof(double)*n_fields*D);
754
755
  double *f0 = malloc(sizeof(double)*n_fields);
  double *f1 = malloc(sizeof(double)*n_fields*D);
756
757
758
759
  double *f00 = malloc(sizeof(double)*n_fields*n_fields);
  double *f10 = malloc(sizeof(double)*n_fields*n_fields*D);
  double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
  double *f11 = malloc(sizeof(double)*n_fields*n_fields*D*D);
Matthieu Constant's avatar
Matthieu Constant committed
760
  problem->kinetic_energy = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
761
  for (int iel=0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
762
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
763
764
    double dxidx[D][D], dphi[N_N][D];
    const double detj = mesh_dxidx(mesh,iel,dxidx);
765
    grad_shape_functions(dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
766
767
    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
768
769
770
    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
771
        dsold[i*D+j] = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
772
773
774
775
776
777
      }
    }
    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
778
779
        double dofold = solution_old[el[i]*n_fields+j];
        for (int k = 0; k < D; ++k) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
780
          ds[j*D+k] += dphi[i][k]*dof;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
781
782
          dsold[j*D+k] += dphi[i][k]*dofold;
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
783
784
785
786
787
788
789
790
791
792
      }
      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];
        }
      }
    }
793
794
795
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
796
797
798
799
      double bf[D];
      for (int j = 0; j < D; ++j) {
        bf[j] = 0;
      }
800
801
802
803
      for (int i = 0; i < n_fields; ++i) {
        s[i] = 0;
        sold[i] = 0;
      }
804
805
806
807
808
809
810
811
812
813
      for (int i = 0; i < n_fields*n_fields; ++i) {
        f00[i] = 0;
      }
      for (int i = 0; i < D*n_fields*n_fields; ++i) {
        f10[i] = 0;
        f01[i] = 0;
      }
      for (int i = 0; i < D*D*n_fields*n_fields; ++i) {
        f11[i] = 0;
      }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
814
815
      double c = 0, a = 0;
      double cold = 0, aold = 0;
816
      for (int i = 0; i < N_SF; ++i) {
Matthieu Constant's avatar
Matthieu Constant committed
817
818
        c += problem->porosity[el[i]]*phi[i];
        cold += problem->oldporosity[el[i]]*phi[i];
819
820
821
822
823
        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
824
        }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
825
826
827
        if (problem->n_fluids==2) {
          a += problem->concentration[iel*N_N+i]*phi[i];
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
828
829
830
        for (int j=0; j<D; ++j) {
          bf[j] += phi[i]*problem->bulk_force[el[i]*D+j];
        }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
831
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
832
833
834
      double mu;
      double rho;
      fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
835
      const double jw = QW[iqp]*detj;
Matthieu Constant's avatar
Matthieu Constant committed
836
      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;
837

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
838
      fluid_problem_f(problem,s,ds,sold,dsold,c,dc,bf,cold,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
839
840
841
842
843
      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
844
845
          }
        }
846
847
848
849
850
      }
      for (int jfield = 0; jfield < n_fields; ++jfield) {
        for (int ifield = 0; ifield < n_fields; ++ifield){
          for (int iphi = 0; iphi < N_SF; ++iphi){
            for (int jphi = 0; jphi < N_SF; ++jphi){
851
              double m = jw*phi[iphi]*phi[jphi]*f00[ifield*n_fields+jfield];
852
              for (int id = 0; id < D; ++id) {
853
                m += jw*dphi[iphi][id]*phi[jphi]*f10[(ifield*D+id)*n_fields+jfield];
854
                for (int jd = 0; jd < D; ++jd) {
855
                  m += jw*dphi[iphi][id]*dphi[jphi][jd]*f11[(ifield*D+id)*n_fields*D+jfield*D+jd];
856
857
                }
              }
858
              for (int jd = 0; jd < D; ++jd) {
859
                m += jw*phi[iphi]*dphi[jphi][jd]*f01[ifield*(n_fields*D)+jfield*D+jd];
Matthieu Constant's avatar
Matthieu Constant committed
860
              }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
861
              LOCAL_MATRIX(iphi,jphi,ifield,jfield) += m;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
862
            }
Matthieu Constant's avatar
Matthieu Constant committed
863
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
864
865
866
        }
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
867
868
869
  }
  free(s);
  free(ds);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
870
  free(dsold);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
871
872
873
  free(sold);
  free(f0);
  free(f1);
874
875
876
877
  free(f00);
  free(f01);
  free(f10);
  free(f11);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
878
879
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
880
static void fluid_problem_surface_tension_bnd(FluidProblem *problem, double *a_cg, double *all_local_vector){
881
  const Mesh *mesh = problem->mesh;
882
  double sigma = problem->sigma;
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
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
  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){
Matthieu Constant's avatar
Matthieu Constant committed
938
              local_vector[cl[iphi]+N_SF*iD] += phi[iphi]*jw*c*c*da[iD]*da[jD]*n[jD]/nda*sigma;
939
            }
Matthieu Constant's avatar
Matthieu Constant committed
940
            local_vector[cl[iphi]+N_SF*iD] -= phi[iphi]*jw*c*c*nda*n[iD]*sigma;
941
942
943
944
945
946
947
          }
        }
      }
    }
  }
}

948
static void fluid_problem_surface_tension(FluidProblem *problem, const double *solution_old, double *grad_a_cg, double *all_local_vector){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
949
950
  if (problem->n_fluids == 1) return;
  const Mesh *mesh = problem->mesh;
951
  double sigma = problem->sigma;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
952
  int n_fields = fluid_problem_n_fields(problem);
953
  double *all_kappa = malloc(sizeof(double)*mesh->n_elements);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
954
955
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
956
957
958
    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
959
960
    grad_shape_functions(dxidx, dphi);
    double da[D] = {0};
Matthieu Constant's avatar
Matthieu Constant committed
961
962
963
    double divnda[D] = {0};
    double divext_da[D] = {0};
    double nda;
Matthieu Constant's avatar
Matthieu Constant committed
964
    double c = 0;
965
966
967
    double kappa = 0;
    double uold[D] = {0};
    double taup = problem->taup[iel];
Matthieu Constant's avatar
Matthieu Constant committed
968
969
    for (int i = 0; i < N_N