fluid_problem.c 70.6 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
/*
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
2
 * MigFlow - Copyright (C) <2010-2020>
3
4
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
Michel Henry's avatar
test    
Michel Henry 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
 * Description and complete License: see LICENSE file.
Michel Henry's avatar
test    
Michel Henry committed
8
9
10
 *
 * This program (MigFlow) is free software:
 * 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
 * Public License as published by the Free Software Foundation, either version
 * 3 of the License, or (at your option) any later version.
Michel Henry's avatar
test    
Michel Henry committed
13
 *
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
14
15
16
 * 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.
Michel Henry's avatar
test    
Michel Henry committed
18
 *
19
 * You should have received a copy of the GNU Lesser General Public License
Michel Henry's avatar
test    
Michel Henry committed
20
 * 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

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

36
37
38
39
40
41
42
static MeshTree *get_mesh_tree(FluidProblem *problem)
{
  if (problem->mesh_tree == NULL)
    problem->mesh_tree = mesh_tree_create(problem->mesh);
  return problem->mesh_tree;
}

43
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double dsol[][D], const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
44
45
46
47
48
49
50
51
52
53
54
55
56
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
57
58
59
60
61
62
63
64
65
66
67
68
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;
Michel Henry's avatar
test    
Michel Henry committed
69
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
70
    return detj/6;
Michel Henry's avatar
test    
Michel Henry committed
71
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
72
73
74
75
76
77
78
79
80
81
}

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
82
int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}
83

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
84
void fluid_problem_node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
85
86
{
  const Mesh *mesh = problem->mesh;
87
  int drag_in_stab = problem->drag_in_stab;
88
  int n_fields = fluid_problem_n_fields(problem);
89
  const FEFields *fields = problem->fields;
90
91
  for (int ip=0; ip < problem->n_particles; ++ip) {
    int iel = problem->particle_element_id[ip];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
92
    if(iel < 0){
93
      for (int d = 0; d < D; ++d)
94
        problem->particle_force[ip*D+d] = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
95
96
      continue;
    }
97
    double *xi = problem->particle_uvw + D*ip;
98
99
100
101
102
    double sf[fields->local_size];
    fe_fields_f(fields, xi, sf);
    double sfporosity[problem->field_porosity->local_size];
    fe_fields_f(problem->field_porosity, xi, sfporosity);
    double dxidx[D][D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
103
    double detj = mesh_dxidx(mesh, iel, dxidx);
104
105
106
107
108
109
110
111
    double dsf[fields->local_size][D];
    fe_fields_df(fields, xi, dxidx, dsf);
    double dsfporosity[problem->field_porosity->local_size][D];
    fe_fields_df(problem->field_porosity, xi, dxidx, dsfporosity);
    double *local_vector = all_local_vector ? &all_local_vector[fields->local_size*iel] : NULL;
    double *local_matrix = all_local_matrix ? &all_local_matrix[fields->local_size*fields->local_size*iel] : NULL;
    double c, cold, dc[D];
    double s[n_fields], ds[n_fields][D];
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
112
    fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
113
    double sold[n_fields];
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
114
115
116
    fe_fields_eval_sf(fields, mesh, iel, sf,solution_old, sold);
    fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
    fe_fields_eval_sf(problem->field_porosity, mesh, iel, sfporosity, problem->oldporosity, &cold);
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
117
    double a = 0;
118
119
120
    if (problem->n_fluids == 2) {
      double sfconcentration[problem->field_concentration->local_size];
      fe_fields_f(problem->field_concentration, xi, sfconcentration);
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
121
      fe_fields_eval_sf(problem->field_concentration, mesh, iel, problem->concentration, sfconcentration, &a);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
122
    }
123
    double f[D],dfdu,dfddp;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
124
    particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
125
126
    if (!local_vector)
      continue;
127
128
    double rho,nu;
    fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
129
    double pspg = drag_in_stab?problem->taup[iel]/rho:0;
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
    double f0[n_fields], f1[n_fields][D], f00[n_fields*n_fields],f11[n_fields*n_fields][D][D],f10[n_fields*n_fields][D],f01[n_fields*n_fields][D];
    for (int i = 0; i < n_fields; ++i) {
      f0[i] = 0;
      for (int id = 0; id < D; ++id) {
        f1[i][id] = 0;
      }
      for (int j = 0; j < n_fields; ++j) {
        int r = i*n_fields+j;
        f00[r] = 0;
        for (int id = 0; id < D; ++id) {
          f01[r][id] = 0;
          f10[r][id] = 0;
          for (int jd = 0; jd < D; ++jd) {
            f11[r][id][jd] = 0;
          }
145
        }
146
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
147
    }
148
149
150
151
152
153
154
155
156
    for (int d = 0; d < D; ++d) {
      f0[U+d] += f[d];
      f1[P][d] += pspg*f[d];
      for (int e = 0; e < D; ++e) {
        double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
        f1[U+d][d] += supg*f[d];
      }
    }
    fe_fields_add_to_local_vector(problem->fields, f0, f1, sf, dsf, 1, local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
157
158
    if (!local_matrix)
      continue;
159
160
161
162
163
164
165
166
167
    for (int d = 0; d < D; ++d) {
      f00[(U+d)*n_fields+U+d] = dfdu;
      f01[(U+d)*n_fields+P][d] = dfddp;
      f10[P*n_fields+U+d][d] = pspg*dfdu;
      f11[P*n_fields+P][d][d] = pspg*dfddp;
      for (int e = 0; e < D; ++e) {
        double supg = drag_in_stab?sold[U+e]*problem->taup[iel]:0;
        f10[(U+d)*n_fields+U+d][e] += supg*dfdu;
        f11[(U+d)*n_fields+P][d][e] += supg*dfddp;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
168
169
      }
    }
170
    fe_fields_add_to_local_matrix(problem->fields, f00, f01, f10, f11, sf, dsf, 1, local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
171
172
173
  }
}

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

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

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

331
static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double dsol[][DIMENSION], const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
332
  const double *contact = &problem->particle_contact[ip*D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
333
  double u[D], uold[D], dp[D];
334
335
336
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
337
    dp[iD] = dsol[P][iD];
338
339
340
341
342
343
344
345
346
347
  }
  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
348
349
350
351
  double rhop = mass/vol;
  double massreduced = mass;
  if(problem->reduced_gravity) {
    massreduced = (rhop-problem->rho[0])*problem->particle_volume[ip];
352
353
  }
  double gamma;
Michel Henry's avatar
test    
Michel Henry committed
354

355
356
357
358
359
360
  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{
361
    gamma = particle_drag_coeff(Due,mu,rho,vol,c);
362
  }
Matthieu Constant's avatar
Matthieu Constant committed
363
364
  gamma = gamma*problem->drag_coeff_factor;
  double gammastar = gamma/(1+dt/mass*gamma);
365
  for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
366
    double upstar = up[i]+dt/mass*(contact[i]+massreduced*problem->g[i]-vol*dp[i]);
367
    double drag = gammastar*(upstar-u[i]/c);
368
    problem->particle_force[ip*D+i] = -drag-vol*dp[i];
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
369
    f[U+i] = -drag;//-vol*dp[i];
370
371
  }
  *dfdu = gammastar/c;
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
372
  *dfddp = gammastar*dt/mass*vol;//-vol;
373
374
}

375
static void fluid_problem_f(const FluidProblem *problem, const double *sol, double dsol[][D], const double *solold, const double dsolold[][D], const double* mesh_velocity, 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[][D], double f00[], double f10[][D], double f01[][D], double f11[][D][D]) {
376
  double p = sol[P];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
377
378
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
379
  const int n_fields = fluid_problem_n_fields(problem);
Michel Henry's avatar
ALE    
Michel Henry committed
380
  double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
381
382
383
384
385
386
387
388
389
390
391
392
393
  for (int i = 0; i < n_fields; ++i) {
    for (int j = 0; j < n_fields; ++j) {
      int r = i*n_fields+j;
      f00[r] = 0;
      for (int id = 0; id < D; ++id) {
        f01[r][id] = 0;
        f10[r][id] = 0;
        for (int jd = 0; jd < D; ++jd) {
          f11[r][id][jd] = 0;
        }
      }
    }
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
394
395
  if(problem->stab_param != 0)
    taup = tauc = 0;
396
  double divu = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
397
  double nold = 0;
398
399
400
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
401
    nold += uold[iD]*uold[iD];
Michel Henry's avatar
ALE    
Michel Henry committed
402
    umesh[iD] = mesh_velocity[iD];
403
    dp[iD] = dsol[P][iD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
404
    for (int jD = 0; jD < D; ++jD) {
405
406
      du[iD][jD] = dsol[U+iD][jD];
      duold[iD][jD] = dsolold[U+iD][jD];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
407
    }
408
    divu += du[iD][iD];
409
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
410
  nold = sqrt(nold);
Nathan Coppin's avatar
Nathan Coppin committed
411
  f0[P] = (c-cold)/dt;// + (sol[P]-solold[P])/dt*0.1;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
412
  //f00[P*n_fields+P] = 1/dt*0.1;
Matthieu Constant's avatar
Matthieu Constant committed
413
  double *g = problem->g;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
  double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
415
  double drag = problem->volume_drag + problem->quadratic_drag*nold;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
416
  for (int i = 0; i < D; ++i) {
Michel Henry's avatar
test    
Michel Henry committed
417
    f0[U+i] =
Jonathan Lambrechts's avatar
cdP    
Jonathan Lambrechts committed
418
      + c*dp[i]
Matthieu Constant's avatar
Matthieu Constant committed
419
      - g[i]*rhoreduced*c
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
420
      - bf[i]*c
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
421
422
      + drag*u[i];
      //+ 5.3e5*u[i];
423
424
    f00[(U+i)*n_fields+U+i] += drag;//5.3e5;
    f01[(U+i)*n_fields+P][i] += c;
425
426
427
428
429
430
431
432
433
434
    if(problem->temporal) {
      f0[U+i] += rho*(u[i]-uold[i])/dt;
      f00[(U+i)*n_fields+U+i] += rho/dt;
    }
    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]));
Michel Henry's avatar
ALE    
Michel Henry committed
435
        f0[U+i] -= rho/c*(umesh[j]*du[i][j]);
436
        f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
437
438
        f01[(U+i)*n_fields+(U+i)][j] += rho/c*uold[j];
        f01[(U+i)*n_fields+(U+i)][j] -= rho/c*umesh[j];
439
      }
440
441
    }
    for (int j = 0; j < D; ++j) {
442
443
444
445
446
      f1[U+i][j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
      f10[(U+i)*n_fields+U+j][j] += -mu*dc[i]/c;
      f10[(U+i)*n_fields+U+i][j] += -mu*dc[j]/c;
      f11[(U+i)*n_fields+U+j][j][i] += mu;
      f11[(U+i)*n_fields+U+i][j][j] += mu;
447
448
      // SUPG
      double supg = uold[j]*taup;
449
450
451
      f1[U+i][j] += supg*f0[U+i];
      f10[(U+i)*n_fields+U+i][j] += supg*f00[(U+i)*n_fields+U+i];
      f11[(U+i)*n_fields+P][j][i] += supg*f01[(U+i)*n_fields+P][i];
452
      for (int k = 0; k < D; ++k)
453
        f11[(U+i)*n_fields+(U+i)][j][k] += supg*f01[(U+i)*n_fields+(U+i)][k];
454
    }
455
    double lsic = tauc*rho;
456
    f1[U+i][i] += ((c-cold)/dt+divu)*lsic;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
457
    for (int j = 0; j < D; ++j) {
458
      f11[(U+i)*n_fields+(U+j)][i][j] += lsic;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
459
    }
460
461
    f1[P][i] = -u[i];
    f10[P*n_fields+U+i][i] += -1;
462
    // PSPG
Jonathan Lambrechts's avatar
pp    
Jonathan Lambrechts committed
463
    double pspg = taup/rho;
464
465
466
    f1[P][i] += pspg*(f0[U+i]+(problem->drag_in_stab?0:(1-c)*dp[i])) + problem->stab_param*dp[i];
    f11[P*n_fields+P][i][i] += pspg*(f01[(U+i)*n_fields+P][i]+(problem->drag_in_stab?0:1-c))+problem->stab_param;
    f10[P*n_fields+U][i] += pspg*f00[(U+i)*n_fields+(U+i)];
467
    for (int j = 0; j < D; ++j) {
468
      f11[P*n_fields+U+i][i][j] = pspg*f01[(U+i)*n_fields+(U+i)][j];
469
    }
Matthieu Constant's avatar
Matthieu Constant committed
470
  }
471
472
}

473
void weak_boundary_values(const FEFields *fields, const Mesh *mesh, const MeshBoundary *bnd, const WeakBoundary *wbnd, double *data);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
474
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
475
476
{
  const Mesh *mesh = problem->mesh;
477
478
479
  const size_t local_size = problem->fields->local_size;
  const FEFields *fields = problem->fields;

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
480
481
482
483
484
485
486
487
488
489
490
491
  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];
492
493
    int n_value = weak_boundary_n_values(wbnd);
    double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
494
    weak_boundary_values(problem->fields, mesh, bnd, wbnd, data);
495

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
496
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
497
      
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
498
499
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
500
501
      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
502
      const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
503
504
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
505
506
507
508
      double n[D];
      double detbnd = get_normal(mesh, interface, n);
      const FEElement *element = mesh->element;
      element->df(QP[0], dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
509
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
510
511
512
513
514
515
        double xiel[D];
        param_boundary_to_element(interface[1],LQP[iqp], xiel);
        double s[fields->n], sold[fields->n], ds[fields->n][D], dsold[fields->n][D];
        double sf[fields->local_size], dsf[fields->local_size][D];
        fe_fields_f(fields, xiel, sf);
        fe_fields_df(fields, QP[iqp], dxidx, dsf);
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
516
517
        fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, problem->solution, s, ds);
        fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
518
519
520
        double dc[D], c, sfporosity[problem->field_porosity->local_size], dsfporosity[problem->field_porosity->local_size][D] ;
        fe_fields_f(problem->field_porosity, QP[iqp], sfporosity);
        fe_fields_df(problem->field_porosity, QP[iqp], dxidx, dsfporosity);
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
521
        fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
522
523
        double f0[fields->n], f00[fields->n*fields->n],f01[fields->n*fields->n][D];
        for (int i = 0; i < fields->n; ++i) {
524
          f0[i] = 0;
525
526
527
528
529
530
531
          for (int j = 0; j < fields->n; ++j) {
            int r = i*fields->n+j;
            f00[r] = 0;
            for (int id = 0; id < D; ++id) {
              f01[r][id] = 0;
            }
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
532
        }
533
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
534
        double a = 0;
535
536
537
        if (problem->n_fluids==2) {
          double sfconcentration[problem->field_concentration->local_size];
          fe_fields_f(problem->field_concentration, QP[iqp], sfconcentration);
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
538
          fe_fields_eval_sf(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
539
540
541
542
543
544
545
546
        }
        const FEElement *mesh_element = problem->mesh->element;
        double meshf[mesh_element->nlocal];
        mesh_element->f(QP[iqp], meshf);
        double umesh[D] = {0};
        for (int i = 0; i < mesh_element->nlocal; ++i) {
          for (int j=0; j<D; ++j) {
            umesh[j] += meshf[i]*problem->mesh_velocity[el[i]*D+j];
Michel Henry's avatar
ALE    
Michel Henry committed
547
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
548
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
549
550
        double rho, mu;
        fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
551
        const double jw = LQW[iqp]*detbnd;
Michel Henry's avatar
ALE    
Michel Henry committed
552
        f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,umesh,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
553
554
        for (int ifield = 0; ifield < D; ++ifield) {
          problem->boundary_force[bndid*D + ifield] += f0[ifield]*jw;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
555
        }
556
557
        fe_fields_add_to_local_vector(fields, f0, NULL, sf, NULL, jw, local_vector);
        fe_fields_add_to_local_matrix(fields, f00, f01, NULL, NULL, sf, dsf, jw, local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
558
      }
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
559
    }
560
    free(data);
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
561
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
562
}
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
563

Matthieu Constant's avatar
Matthieu Constant committed
564
565
566
567
568
569
570
571
572
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)
Michel Henry's avatar
test    
Michel Henry committed
573
      for (int j = 0; j< N_N; ++j)
Matthieu Constant's avatar
Matthieu Constant committed
574
575
576
577
578
        I_a += detj * problem->porosity[el[i]]*problem->concentration[iel*N_N+j]*mass_matrix[i][j];
  }
  return I_a;
}

579
580
double fluid_problem_volume_flux(FluidProblem *problem, const char *tagname)
{
581
582
#if 0
  const FEElement *element = problem->fields->element[0];
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
  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];
599
600
      const int *cl = mesh->element->boundaries[icl];
      int nodes[element->boundary_element->nlocal];
601
602
603
604
605
606
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
607
608
      double n[D];
      double detbnd = get_normal(mesh, interface,n);
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
      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;
624
625
#endif
  return 0;
626
627
}

Matthieu Constant's avatar
Matthieu Constant committed
628
629
double fluid_problem_a_integ_bnd(FluidProblem *problem, double dt)
{
630
631
#if 0
  const FEElement *element = problem->fields->element[0];
Matthieu Constant's avatar
Matthieu Constant committed
632
633
634
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  const int n_fields = fluid_problem_n_fields(problem);
635
  const size_t local_size = element->nlocal*n_fields;
Matthieu Constant's avatar
Matthieu Constant committed
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
  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));
652
    weak_boundary_values(problem->fields, mesh, bnd, wbnd, data);
Matthieu Constant's avatar
Matthieu Constant committed
653
654
655
656
657

    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];
658
659
660
      const int *cl = element->boundaries[icl];
      int nodes[element->boundary_element->nlocal];
      double *x[element->boundary_element->nlocal];
Matthieu Constant's avatar
Matthieu Constant committed
661
662
663
664
665
666
667
668
669
      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);
670
      element->df(QP[0], dxidx, dphi);
Michel Henry's avatar
test    
Michel Henry committed
671

Matthieu Constant's avatar
Matthieu Constant committed
672
673
674
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
        double phi[DIMENSION];
675
        element->boundary_element->f(LQP[iqp],phi);
Matthieu Constant's avatar
Matthieu Constant committed
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
        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;
704
705
#endif
  return 0;
Matthieu Constant's avatar
Matthieu Constant committed
706
707
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
708
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
709
710
{
  const Mesh *mesh = problem->mesh;
Matthieu Constant's avatar
Matthieu Constant committed
711
  const double *porosity = problem->porosity;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
712
713
  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
714
715
716
717
718
719
720
721
  const FEFields *fields = problem->fields;
  double mesh_velocity[D];
  double f0[n_fields];
  double f1[n_fields][D];
  double f00[n_fields*n_fields];
  double f10[n_fields*n_fields][D];
  double f01[n_fields*n_fields][D];
  double f11[n_fields*n_fields][D][D];
Matthieu Constant's avatar
Matthieu Constant committed
722
  problem->kinetic_energy = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
723
  for (int iel=0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
724
    const int *el = &mesh->elements[iel*N_N];
725
    double dxidx[D][D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
726
    const double detj = mesh_dxidx(mesh,iel,dxidx);
727
728
729
730
731
732
733
    double *local_vector = &all_local_vector[fields->local_size*iel];
    double *local_matrix = &all_local_matrix[fields->local_size*fields->local_size*iel];
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double s[fields->n], sold[fields->n], ds[fields->n][D], dsold[fields->n][D];
      double sf[fields->local_size], dsf[fields->local_size][D];
      fe_fields_f(fields, QP[iqp], sf);
      fe_fields_df(fields, QP[iqp], dxidx, dsf);
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
734
735
      fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, solution, s, ds);
      fe_fields_eval_grad_sf(fields, mesh, iel, sf, dsf, solution_old, sold, dsold);
736
737
738
739
      double sfporosity[fields->local_size], dsfporosity[fields->local_size][D];
      fe_fields_f(problem->field_porosity, QP[iqp], sfporosity);
      fe_fields_df(problem->field_porosity, QP[iqp], dxidx, dsfporosity);
      double c, dc[D], cold;
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
740
741
      fe_fields_eval_grad_sf(problem->field_porosity, mesh, iel, sfporosity, dsfporosity, problem->porosity, &c, &dc);
      fe_fields_eval_sf(problem->field_porosity, mesh, iel, sfporosity, problem->porosity, &cold);
742
      double a = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
743
      if (problem->n_fluids==2) {
744
745
        double sfconcentration[problem->field_concentration->local_size];
        fe_fields_f(problem->field_concentration, QP[iqp], sfconcentration);
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
746
        fe_fields_eval_sf(problem->field_concentration, mesh, iel, sfconcentration, problem->concentration, &a);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
747
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
748
      double bf[D];
749
750
751
752
      const FEElement *mesh_element = problem->mesh->element;
      double meshf[mesh_element->nlocal];
      mesh_element->f(QP[iqp], meshf);
      for (int i = 0; i < mesh_element->nlocal; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
753
        for (int j=0; j<D; ++j) {
754
755
          mesh_velocity[j] += meshf[i]*problem->mesh_velocity[el[i]*D+j];
          bf[j] += meshf[i]*problem->bulk_force[el[i]*D+j];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
756
        }
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
757
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
758
759
760
      double mu;
      double rho;
      fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
761
      const double jw = QW[iqp]*detj;
762
763
      for (int d = 0; d < D; ++d)
        problem->kinetic_energy += rho*(s[U+d]*s[U+d])/2*jw;
Michel Henry's avatar
ALE    
Michel Henry committed
764
      fluid_problem_f(problem,s,ds,sold,dsold,mesh_velocity,c,dc,bf,cold,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
765
766
      fe_fields_add_to_local_vector(fields, f0, f1, sf, dsf, jw, local_vector);
      fe_fields_add_to_local_matrix(fields, f00, f01, f10, f11, sf, dsf, jw, local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
767
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
768
769
770
  }
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
771
static void fluid_problem_surface_tension_bnd(FluidProblem *problem, double *a_cg, double *all_local_vector){
772
#if 0
773
  const Mesh *mesh = problem->mesh;
774
  double sigma = problem->sigma;
775
  const int n_fields = fluid_problem_n_fields(problem);
776
777
  const FEElement *element = problem->fields->element[0];
  const size_t local_size = element->nlocal*n_fields;
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
  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];
795
      const int *cl = element->boundaries[icl];
796
797
798
799
      int nodes[D];
      double *x[D];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
800
        nodes[i] = el[element->boundaries[icl][i]];
801
802
803
804
805
806
        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);
807
      element->df(QP[0], dxidx, dphi);
808
809
810
      double da[D] = {0};
      double nda = 0;
      for (int k=0; k<D; ++k){
811
        for (int i = 0; i < element->nlocal; ++i) {
812
813
814
815
816
817
818
819
          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];
820
        element->boundary_element->f(LQP[iqp],phi);
821
822
823
824
825
826
827
828
829
830
        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){
831
              local_vector[cl[iphi]+element->nlocal*iD] += phi[iphi]*jw*c*c*da[iD]*da[jD]*n[jD]/nda*sigma;
832
            }
833
            local_vector[cl[iphi]+element->nlocal*iD] -= phi[iphi]*jw*c*c*nda*n[iD]*sigma;
834
835
836
837
838
          }
        }
      }
    }
  }
839
#endif
840
841
}

842
static void fluid_problem_surface_tension(FluidProblem *problem, const double *solution_old, double *grad_a_cg, double *all_local_vector){
843
  const FEElement *element = problem->fields->element[0];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
844
845
  if (problem->n_fluids == 1) return;
  const Mesh *mesh = problem->mesh;
846
  double sigma = problem->sigma;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
847
  int n_fields = fluid_problem_n_fields(problem);
848
  double *all_kappa = malloc(sizeof(double)*mesh->n_elements);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
849
850
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
851
852
853
    double dxidx[D][D], dphi[N_N][D];
    const double det = mesh_dxidx(mesh,iel,dxidx);
    const double vol = element_volume_from_detj(det);
854
    element->df(QP[0], dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
855
    double da[D] = {0};
Matthieu Constant's avatar
Matthieu Constant committed
856
857
858
    double divnda[D] = {0};
    double divext_da[D] = {0};
    double nda;
Matthieu Constant's avatar
Matthieu Constant committed
859
    double c = 0;
860
861
862
    double kappa = 0;
    double uold[D] = {0};
    double taup = problem->taup[iel];
Matthieu Constant's avatar
Matthieu Constant committed
863
864
    for (int i = 0; i < N_N; ++i){
      c += problem->porosity[el[i]]/N_N;
Matthieu Constant's avatar
Matthieu Constant committed
865
866
867
868
869
870
871
872
      nda = 0;
      for (int iD = 0; iD < D; iD++){
        nda += grad_a_cg[el[i]*D+iD]*grad_a_cg[el[i]*D+iD];
        for (int jD = 0; jD < D; jD++){
          divext_da[iD] += grad_a_cg[el[i]*D+iD]*grad_a_cg[el[i]*D+jD]*dphi[i][jD];
        }
      }
      nda = fmax(sqrt(nda),1e-8);
873
      //if (nda>1e-8) printf("nda = %.4e\n",nda);
Matthieu Constant's avatar
Matthieu Constant committed
874
875
876
      for (int iD = 0; iD < D; iD++){
        divnda[iD] += nda*dphi[i][iD];
        divext_da[iD] /= nda;
877
        kappa += -grad_a_cg[el[i]*D+iD]*dphi[i][iD]/nda;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
878
      }
879
        //if (nda>1e-8) printf("kappa[%d] = %.4e\n",i,kappa[i]);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
880
    }
881
    all_kappa[iel] = kappa;
Matthieu Constant's avatar
Matthieu Constant committed
882
    for (int iqp = 0; iqp< N_QP; ++iqp) {
883
884
      double phi[element->nlocal];
      element->f(QP[iqp],phi);
Matthieu Constant's avatar
Matthieu Constant committed
885
      double *local_vector = &all_local_vector[iel*n_fields*N_N];
886
887
888
889
      double a = 0;
      for (int id = 0; id < D; id++){
        da[id] = 0;
      }
890
      for (int i = 0; i < element->nlocal; ++i) {
891
892
893
894
895
896
897
898
899
        for (int id = 0; id < D; ++id) {
          da[id] += phi[i]*grad_a_cg[el[i]*D+id];
          uold[id] += phi[i]*solution_old[el[i]*n_fields+id];
        }
        a += problem->concentration[iel*N_N+i]*phi[i];
      }
      double mu;
      double rho;
      fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);const double jw = QW[iqp]*det;
900
      for (int iphi = 0; iphi < element->nlocal; ++iphi) {
Matthieu Constant's avatar
Matthieu Constant committed
901
        for (int id = 0; id < D; ++id) {
Matthieu Constant's avatar
Matthieu Constant committed
902
          //local_vector[(U+id)*N_N+iphi] += c*sigma*jw*(divnda[id]-divext_da[id])*phi[iphi];
903
          local_vector[(U+id)*element->nlocal+iphi] += -c*sigma*jw*kappa*da[id]*phi[iphi];