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

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

Matthieu Constant's avatar
Matthieu Constant committed
33

34

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

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

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

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

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
128
129
#endif

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
130
int fluid_problem_n_fields(const FluidProblem *p) {return p->n_fluids*(D+1)+1;}
131

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

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

static void compute_node_force_implicit(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix) {
166
#if 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
167
168
169
170
  double *porosity = problem->porosity;
  Mesh *mesh = problem->mesh;
  for (int i = 0; i < problem->n_particles; ++i) {
    int iel = problem->particle_element_id[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
171
    double mass = problem->particle_mass[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
172
    if(iel < 0){
173
174
175
      for (int d = 0; d < D; ++d)
        problem->particle_force[i*D+d] = 0;
      problem->particle_force[i*D+1] = problem->g*mass;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
176
177
      continue;
    }
178
    double *xi = problem->particle_uvw + D*i;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
179
180
181
182
    double phi[N_SF];
    shape_functions(xi,phi);
    uint32_t *el=problem->mesh->elements+iel*N_SF;
    double *se = problem->solution_explicit;
183
184
185
186
    double c=0, vf[D]={0}, vfe[D]={0}, gradp[D]={0};
    double dxdxi[D][D],dxidx[D][D];
    for (int k=0; k<D; ++k)
      for (int j=0; j<D; ++j)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
187
188
        dxdxi[k][j] = mesh->x[el[j+1]*3+k]-mesh->x[el[0]*3+k];
    invDD(dxdxi,dxidx);
189
    double dphi[N_SF][D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
190
    grad_shape_functions(dxidx, dphi);
191
    double *si = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
192
193
    for (int iphi=0; iphi < N_SF; ++iphi) {
      c += porosity[el[iphi]]*phi[iphi];
194
      for (int j=0; j < D; ++j) {
195
        vf[j]  += si[el[iphi]*n_fields+j]*phi[iphi];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
196
        vfe[j] += se[el[iphi]*n_fields+j]*phi[iphi];
197
        gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+D];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
198
199
      }
    }
200
201
202
203
    double du[D], due[D];
    for (int j = 0; j < D; ++j) {
      du[j] = problem->particle_velocity[i*D+j]-vf[j]/c;
      due[j] = problem->particle_velocity[i*D+j]-vfe[j]/c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
204
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
205
    double vol = problem->particle_volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
206
207
    double rhop = mass/vol;
    double drho = rhop -problem->rho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
208
    double g = problem->g*drho/rhop;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
209
    double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
210
211
    double fcoeff = mass/(mass+dt*gamma);
    
212
213
214
    for (int d = 0; d < D; ++d)
      problem->particle_force[i*D+d] = fcoeff*(-gamma*du[d]-vol*gradp[d]);
    problem->particle_force[i*D+1] += fcoeff*g*mass;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
215
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
216
217
    double *local_vector = all_local_vector+iel*n_fields*N_SF;
    double *local_matrix = all_local_matrix+iel*n_fields*n_fields*N_SF*N_SF;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
218
    for (int iphi = 0; iphi < N_SF; ++iphi) {
219
      for (int d = 0; d < D; ++d)
220
221
222
        local_vector[iphi+N_SF*d] -= fcoeff*gamma*(du[d]-dt/mass*vol*gradp[d])*phi[iphi];
      local_vector[iphi+N_SF*1] -= fcoeff*gamma*dt*g*phi[iphi];
      
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
223
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
224
225
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int jphi = 0; jphi < N_SF; ++jphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
226
227
228
        int N2 = n_fields*N_SF;
        int IDX = (iphi*N2+jphi)*n_fields;
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
229
        double f = fcoeff*phi[iphi];
230
        for (int d = 0; d < D; ++d){
231
          LOCAL_MATRIX(d,d) -= -f/c*phi[jphi]*gamma;
232
          LOCAL_MATRIX(d,D) -= -f*gamma*dt/mass*vol*dphi[jphi][d];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
233
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
234
235
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
236
  }
237
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
238
239
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
240
typedef void f_cb(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double dt);
241

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
242
243
244
static void f_inflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt){
  int n_fluids = problem->n_fluids;
  int P = n_fluids*(D+1);
245
  const double epsilon = problem->epsilon;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
246
  double p = s[P];
247
248
249
250
  double dp[D];
  for (int id = 0; id < D; ++id) {
    dp[id] = ds[P*D+id];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
251
252
253
254
255
256
  for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
    int U = ifluid*(D+1);
    int Q = ifluid*(D+1)+D;
    double q = s[Q];
    f[Q] = 0;
    for (int id = 0; id < D; ++id) {
257
258
      f[Q] += -s[U+id]*n[id];
      f[Q] += dp[id]*n[id]*epsilon*q;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
259
260
261
262
      f[U+id] = 0;//-p*q*n[id];
    }
  }
  f[P] = 0;
263
264
265
  for (int id = 0; id < D; ++id) {
    //f[P] += dp[id]*n[id]*epsilon;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
266
267
268
269
270
}

static void f_outflow(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt){
  int n_fluids = problem->n_fluids;
  int P = n_fluids*(D+1);
271
  const double epsilon = problem->epsilon;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
272
  double p = s[P];
273
274
275
276
  double dp[D];
  for (int id = 0; id < D; ++id) {
    dp[id] = ds[P*D+id];
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
277
278
279
280
281
282
  for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
    int U = ifluid*(D+1);
    int Q = ifluid*(D+1)+D;
    double q = s[Q];
    f[Q] = 0;
    for (int id = 0; id < D; ++id) {
283
284
      f[Q] += -s[U+id]*n[id];
      f[Q] += dp[id]*n[id]*epsilon*q;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
285
286
287
288
      f[U+id] = 0;//-p*q*n[id];
    }
  }
  f[P] = 0;
289
290
291
  for (int id = 0; id < D; ++id) {
    //f[P] += dp[id]*n[id]*epsilon;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
292
293
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
static void f_wall_pressure(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds,const double c,const double dt)
{
  int n_fluids = problem->n_fluids;
  int P = n_fluids*(D+1);
  const double epsilon = problem->epsilon;
  double dp[D];
  for (int id = 0; id < D; ++id) {
    dp[id] = ds[P*D+id];
  }
  double p = s[P];
  for (int ifluid = 0; ifluid < n_fluids; ++ifluid) {
    int U = ifluid*(D+1);
    int Q = ifluid*(D+1)+D;
    double q = s[Q];
    f[Q] = 0;
    for (int id = 0; id < D; ++id) {
310
      f[Q] += dp[id]*n[id]*epsilon*q;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
311
      f[U+id] = q*p*n[id];
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
312
313
    }
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
314
315
  f[P] = 0;
  for (int id = 0; id < D; ++id) {
316
    //f[P] += dp[id]*n[id]*epsilon;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
317
  }
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
318
}
319

320
321
322
323
static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, double c, double dt) {
  int P = problem->n_fluids*(D+1);
  double p = sol[P];
  double dp[D] = {dsol[P*D+0],dsol[P*D+1]};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
324
325
  double epsilonp = problem->epsilon;
  double epsilonq = problem->epsilon/10;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
326
  f0[P] = 1-c;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
327
  for (int id = 0; id < D; ++id) {
328
    f1[P*D+id] = 0;//epsilonp*dp[id];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
329
  }
330
331
332
333
334
335
336
337
338
339
340
  for (int ifluid=0; ifluid < problem->n_fluids; ++ifluid) {
    int Q = ifluid*(D+1)+D;
    int U = ifluid*(D+1);
    double rho = problem->rho[ifluid];
    double mu = problem->mu[ifluid];
    double q = sol[Q];
    double qold = solold[Q];
    double dq[D] = {dsol[Q*D+0],dsol[Q*D+1]};
    double u[D] = {sol[U],sol[U+1]};
    double uold[D] = {solold[U],solold[U+1]};
    double du[D][D] = {{dsol[U*D+0],dsol[U*D+1]},{dsol[U*D+D+0],dsol[U*D+D+1]}};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
341

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
342
343
    double oq2 = 1./fmax(q,1e-5);
    double oq = 1./fmax(q,1e-2);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
344
345
346
347
348
349
350
351
352
353
    double divu = 0;
    double tau[D][D];
    double utau[D]={0.};
    for (int i = 0; i < D; ++i) {
      divu += du[i][i];
      for (int j = 0; j < D; ++j) {
        tau[i][j] = du[i][j]-u[i]*dq[j]*oq;
        utau[i] += u[j]*tau[i][j];
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
354
    f0[Q] = (q-qold)/dt;
355
    f0[P] -= q;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
356
    for (int i = 0; i < D; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
357
358
      //f0[Q] += epsilonp*dq[i]*dp[i];
      f0[U+i] = rho*((u[i]-uold[i])/dt + u[i]*oq2*divu + utau[i]*oq2) -p*dq[i]+ (i==(D-1) ? -q*rho*problem->g : 0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
359
360
361
      for (int j = 0; j < D; ++j) {
        f1[(U+j)*D+i] = mu*(tau[i][j]+tau[j][i]) + (i==j ? -q*p : 0);
      }
362
      f1[Q*D+i] = epsilonq*dq[i]-u[i]+q*epsilonp*dp[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
363
    }
364
365
366
  }
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
367
static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
368
369
370
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
371
372
373
374
375
376
377
378
379
380
  const int n_fields = fluid_problem_n_fields(problem);
  const size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *g0 = malloc(sizeof(double)*n_fields);
  double *h0 = malloc(sizeof(double)*n_fields*D);
  f_cb *f;
  for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
381

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
382
383
384
385
386
387
388
389
    if (strcmp(mesh->phys_name[iphys],"Bottom") == 0)
      f = f_wall_pressure;
    else if (strcmp(mesh->phys_name[iphys],"Right") == 0)
      f = f_wall_pressure;
    else if (strcmp(mesh->phys_name[iphys],"Left") == 0)
      f = f_wall_pressure;
    else if (strcmp(mesh->phys_name[iphys],"Top") == 0)
      f = f_wall_pressure;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
390
391
392
393
394
395
396
397
    else if (strcmp(mesh->phys_name[iphys],"LeftUp") == 0)
      f = f_inflow;
    else if (strcmp(mesh->phys_name[iphys],"LeftDown") == 0)
      f = f_inflow;
    else if (strcmp(mesh->phys_name[iphys],"RightUp") == 0)
      f = f_outflow;
    else if (strcmp(mesh->phys_name[iphys],"RightDown") == 0)
      f = f_outflow;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
    else
      continue;

    for (int ibnd = 0; ibnd < mesh->phys_n_boundaries[iphys]; ++ibnd) {
      const int *bnd = &mesh->phys_boundaries[iphys][ibnd*2];
      const int iel = bnd[0];
      const int i0 = bnd[1];
      const int i1 = (bnd[1]+1) % 3;
      const unsigned int *el = &mesh->elements[iel*N_N];

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

      double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
      for (int i = 0; i < D; ++i)
        for (int j = 0; j < D; ++j)
          dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
      invDD(dxdxi, dxidx);
      grad_shape_functions(dxidx, dphi);
      const double l = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
      const double n[2] = {-dx[1]/l,dx[0]/l};

      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double xi[2];
        double a = (1+LQP[iqp])/2;
        switch(i0) {
          case 0 : xi[0] = a; xi[1] = 0;   break;
          case 1 : xi[0] = 1-a; xi[1] = a;   break;
          case 2 : xi[0] = 0;   xi[1] = 1-a; break;
          default:;
        }
        double dp[D]={0};
        double *local_vector = &all_local_vector[local_size*iel];
        double *local_matrix = &all_local_matrix[local_size*local_size*iel];
        double phi[N_SF];
        shape_functions(xi,phi);
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
          sold[i] = 0;
          for (int j = 0; j < D; ++j) {
            ds[i*D+j] = 0;
          }
        }
        double c = 0;
        for (int i = 0; i < N_SF; ++i) {
          c += problem->compacity[el[i]]*phi[i];
          for (int j = 0; j < n_fields; ++j) {
            double dof = solution[el[i]*n_fields+j];
            s[j] += phi[i]*dof;
            for (int k = 0; k < D; ++k)
              ds[j*D+k] += dphi[i][k]*dof;
          }
        }
        const double jw = LQW[iqp]*l/2;
        f(problem,n,f0,s,ds,c,dt);
        for (int ifield = 0; ifield < n_fields; ++ifield) {
          for (int iphi = 0; iphi < N_SF; ++iphi){
            local_vector[iphi+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
          }
        }
        for (int jfield = 0; jfield < n_fields; ++jfield) {
          double store = s[jfield];
          s[jfield] += deps;
          f(problem,n,g0,s,ds,c,dt);
          s[jfield] = store;
          for (int jd =0; jd < D; ++jd){
            store = ds[jfield*D+jd];
            ds[jfield*D+jd] += deps;
            f(problem,n,h0+jd*n_fields,s,ds,c,dt);
            ds[jfield*D+jd] = store;
          }
          int N2 = n_fields*N_SF;
#define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)]
          for (int ifield = 0; ifield < n_fields; ++ifield){
            for (int iphi = 0; iphi < N_SF; ++iphi){
              for (int jphi = 0; jphi < N_SF; ++jphi){
                double df00 = (g0[ifield]-f0[ifield])/deps;
                double m = jw*phi[iphi]*phi[jphi]*df00;
                for (int jd = 0; jd < D; ++jd) {
                  double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
                  m += jw*phi[iphi]*dphi[jphi][jd]*df01;
                }
                int IDX = (iphi*N2)*n_fields + jphi*n_fields;
                LOCAL_MATRIX2(ifield,jfield) += m;
              }
            }
          }
        }
      }
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
488
489
    }
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
490
}
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
491

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
492
493
494
495
496
497
498
499
500
501
502

static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
  const Mesh *mesh = problem->mesh;
  const double *compacity = problem->compacity;
  const double epsilon = problem->epsilon;

  const double *solution = problem->solution;
  int n_fluids = problem->n_fluids;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
503
504
505
506
507
508
509
510
511
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *f1 = malloc(sizeof(double)*n_fields*D);
  double *g0 = malloc(sizeof(double)*n_fields);
  double *g1 = malloc(sizeof(double)*n_fields*D);
  double *h0 = malloc(sizeof(double)*n_fields*D);
  double *h1 = malloc(sizeof(double)*n_fields*D*D);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
512
  for (int iel=0; iel < mesh->n_elements; ++iel) {
513
    const unsigned int *el = &mesh->elements[iel*N_N];
514
515
516
    double dxdxi[D][D], dxidx[D][D], dphi[N_N][D];
    for (int i = 0; i < D; ++i)
      for (int j = 0; j < D; ++j)
517
518
519
        dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
    const double detj = invDD(dxdxi, dxidx);
    grad_shape_functions(dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
520
521
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
522
523
524
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
525
526
527
528
529
530
531
      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;
        }
      }
532
533
534
      double c = 0;
      for (int i = 0; i < N_SF; ++i) {
        c += problem->compacity[el[i]]*phi[i];
535
536
537
538
539
540
541
542
        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;
        }
543
544
      }
      const double jw = QW[iqp]*detj;
545
546
547
548
549
550
      fluid_problem_f(problem,f0,f1,s,ds,sold,c,dt);
      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
551
552
          }
        }
553
554
555
556
557
558
559
560
561
562
563
      }
      for (int jfield = 0; jfield < n_fields; ++jfield) {
        double store = s[jfield];
        s[jfield] += deps;
        fluid_problem_f(problem,g0,g1,s,ds,sold,c,dt);
        s[jfield] = store;
        for (int jd =0; jd < D; ++jd){
          store = ds[jfield*D+jd];
          ds[jfield*D+jd] += deps;
          fluid_problem_f(problem,h0+jd*n_fields,h1+jd*D*n_fields,s,ds,sold,c,dt);
          ds[jfield*D+jd] = store;
564
        }
565
566
567
568
569
570
571
572
573
574
575
576
577
        int N2 = n_fields*N_SF;
#define LOCAL_MATRIX2(a,b) local_matrix[IDX+N2*(a)+(b)]
        for (int ifield = 0; ifield < n_fields; ++ifield){
          for (int iphi = 0; iphi < N_SF; ++iphi){
            for (int jphi = 0; jphi < N_SF; ++jphi){
              double df00 = (g0[ifield]-f0[ifield])/deps;
              double m = jw*phi[iphi]*phi[jphi]*df00;
              for (int id = 0; id < D; ++id) {
                double df10 = (g1[ifield*D+id]-f1[ifield*D+id])/deps;
                m += jw*dphi[iphi][id]*phi[jphi]*df10;
                for (int jd = 0; jd < D; ++jd) {
                  double df11 = (h1[jd*n_fields*D+ifield*D+id]-f1[ifield*D+id])/deps;
                  m += jw*dphi[iphi][id]*dphi[jphi][jd]*df11;
578
579
                }
              }
580
581
582
              for (int jd = 0; jd < D; ++jd) {
                double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
                m += jw*phi[iphi]*dphi[jphi][jd]*df01;
Matthieu Constant's avatar
Matthieu Constant committed
583
              }
584
585
              int IDX = (iphi*N2)*n_fields + jphi*n_fields;
              LOCAL_MATRIX2(ifield,jfield) += m;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
586
            }
Matthieu Constant's avatar
Matthieu Constant committed
587
          }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
588
589
590
        }
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
  }

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

static void fluid_problem_assemble_system(FluidProblem *problem, double *rhs, const double *solution_old, double dt)
{
  HXTLinearSystem *lsys = problem->linear_system;
  const Mesh *mesh = problem->mesh;

  int n_fluids = problem->n_fluids;
  int n_fields = fluid_problem_n_fields(problem);
  
  size_t local_size = N_SF*n_fields;
  size_t N_E = mesh->n_elements;
  double *all_local_vector = (double*) malloc(N_E*local_size*sizeof(double));
  double *all_local_matrix = (double*) malloc(N_E*local_size*local_size*sizeof(double));
  for (size_t i = 0; i < N_E*local_size; ++i)
    all_local_vector[i] = 0;
  for (size_t i = 0; i < N_E*local_size*local_size; ++i)
    all_local_matrix[i] = 0;

  char *forced_row = malloc(sizeof(char)*n_fields*mesh->n_nodes);
  for (int i = 0; i < n_fields*mesh->n_nodes; ++i)
    forced_row[i] = -1;
  for (int ibnd = 0; ibnd < problem->n_strong_boundaries; ++ibnd) {
    const StrongBoundary *bnd = problem->strong_boundaries + ibnd;
    int iphys;
    for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
      if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_phys)
      printf("Boundary tag \"%s\" not found.", bnd->tag);
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      forced_row[mesh->phys_nodes[iphys][i]*n_fields + bnd->row] = bnd->field;
    }
  }

  //compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
  compute_weak_boundary_conditions(problem, dt, all_local_vector, all_local_matrix);
  fluid_problem_volume(problem, solution_old, dt, all_local_vector, all_local_matrix);
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
644
645
    for (int inode = 0; inode < N_SF; ++inode) {
      for(int ifield = 0; ifield < n_fields; ++ifield) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
646
        const unsigned int *el = &mesh->elements[iel*N_N];
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
647
648
649
650
651
652
653
654
655
656
657
658
        char forced_field = forced_row[el[inode]*n_fields+ifield];
        if (forced_field == -1) continue;
        int i = inode*n_fields + ifield;
        for (int jnode = 0; jnode< N_SF; ++jnode){
          for (int jfield= 0; jfield< n_fields; ++jfield){
            local_matrix[(inode*n_fields+ifield)*local_size+(jnode*n_fields+jfield)] =
              (inode==jnode && jfield == forced_field) ? 1. : 0.;
          }
        }
        local_vector[inode+ifield*N_SF] = 0;
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
659
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
660
661
662
    /*for (int i = 0; i < local_size; ++i)
      printf("%8.2g ", local_vector[i]);
    printf("\n");*/
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
663
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
664
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
665
  //exit(0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
666
  free(forced_row);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
667
668
  free(all_local_matrix);
  free(all_local_vector);
669
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
670

671
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int n_fields, int nBnd, StrongBoundary *bnds)
672
{
673
674
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
675
676
677
678
679
680
681
682
    int iphys;
    for (iphys = 0; iphys < mesh->n_phys; ++iphys) {
      if (strcmp(bnd->tag, mesh->phys_name[iphys]) == 0)
        break;
    }
    if (iphys == mesh->n_phys){
      printf("Boundary tag \"%s\" not found.", bnd->tag);
    }
683
684
    double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
    double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
685
686
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
687
688
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
689
690
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
691
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
692
693
694
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
695
696
697
  }
}

698
int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
699
{
700
701
  static double totaltime = 0;
  struct timespec tic, toc;
702
  const Mesh *mesh = problem->mesh;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
703
  int n_fields = fluid_problem_n_fields(problem);
704
705
706
  apply_boundary_conditions(mesh, problem->solution, n_fields, problem->n_strong_boundaries, problem->strong_boundaries);
  double *solution_old = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
  double *rhs = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
707
708
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
709
  problem->solution_explicit = solution_old;
710
  double *dx = (double*)malloc(sizeof(double)*mesh->n_nodes*n_fields);
711
  double norm0 = 0;
712
713
  int i;
  for (i=0; i<newton_max_it; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
714
715
    double norm = 0;
    hxtLinearSystemZeroMatrix(problem->linear_system);
716
717

    for (int i=0; i<mesh->n_nodes*n_fields; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
718
      rhs[i] = 0;
719
720
721
      //printf("solution=%g\n",problem->solution[i]);
    }
    //break;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
722
723
    fluid_problem_assemble_system(problem, rhs, solution_old, dt);
    hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
724
725
726
727
728
729
730
    printf("iter %i norm = %g\n", i, norm);
    if (norm < newton_atol)
      break;
    if (i == 0)
      norm0 = norm;
    else if(norm/norm0 < newton_rtol)
      break;
731
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
732
  hxtLinearSystemSolve(problem->linear_system, rhs, dx);
733
  clock_get(&toc);
734
  //break;
735
  totaltime += clock_diff(&tic, &toc);
736
737
738
739
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
740
  printf("total solve time : %g\n", totaltime);
741
  free(dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
742
  free(rhs);
743
  free(solution_old);
744
  return i;
745
746
}

Matthieu Constant's avatar
Matthieu Constant committed
747
748
749
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
  free(problem->node_volume);
  Mesh *mesh = problem->mesh;
750
  problem->node_volume = (double*)malloc(sizeof(double)*mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
751
752
753
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->node_volume[i] = 0;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
754
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
755
    double dxdxi[D][D];
756
    const uint32_t *el = &mesh->elements[iel*N_N];
757
758
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
759
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
760
#if D == 2
761
    const double vol = detDD(dxdxi)/6;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
762
#else
763
    const double vol = detDD(dxdxi)/24;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
764
#endif
765
766
    for (int i = 0; i < N_N; ++i)
      problem->node_volume[el[i]] += vol;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
767
768
769
  }
}

770
static void fluid_problem_compute_compacity(FluidProblem *problem)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
771
772
773
774
{
  Mesh *mesh = problem->mesh;
  double *volume = problem->particle_volume;
  for (int i = 0; i < mesh->n_nodes; ++i){
775
    problem->compacity[i] = 0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
776
777
778
  }
  for (int i = 0; i < problem->n_particles; ++i) {
    if(problem->particle_element_id[i] == -1) continue;
779
    double sf[N_SF];
780
    shape_functions(&problem->particle_uvw[i*D],sf);
781
    const uint32_t *el = &mesh->elements[problem->particle_element_id[i]*N_N];
782
    for (int j=0; j<N_N;++j)
783
      problem->compacity[el[j]] += sf[j]*volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
784
785
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
786
    if(problem->node_volume[i]==0.){
787
      problem->compacity[i] = 0.;
Matthieu Constant's avatar
Matthieu Constant committed
788
789
    }
    else{
790
      problem->compacity[i] = problem->compacity[i]/problem->node_volume[i];
Matthieu Constant's avatar
Matthieu Constant committed
791
    }
Matthieu Constant's avatar
Matthieu Constant committed
792
793
794
  }
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
795
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double *mu, double *rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries, int n_fluids)
796
{
797
798
799
800
801
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
802
    hxtPETScInsertOptions("-pc_type lu -ksp_max_it 30", "fluid");
803
804
#endif
  }
805
  FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
806
807
808
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
809
  problem->mesh_tree = mesh_tree_create(mesh);
Matthieu Constant's avatar
Matthieu Constant committed
810

Matthieu Constant's avatar
Matthieu Constant committed
811
  problem->n_fluids = n_fluids;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
812
  int n_fields = fluid_problem_n_fields(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
813
814
  problem->mu = (double*)malloc(sizeof(double)*n_fluids);
  problem->rho = (double*)malloc(sizeof(double)*n_fluids);
Matthieu Constant's avatar
Matthieu Constant committed
815
816
817
818
  for (int i = 0; i < n_fluids; ++i){
    problem->rho[i] = rho[i];
    problem->mu[i] = mu[i];
  }
819
  problem->g = g;
820
  problem->mesh = mesh;
Matthieu Constant's avatar
Matthieu Constant committed
821
  problem->epsilon = epsilon;
822
  problem->strong_boundaries = (StrongBoundary*)malloc(sizeof(StrongBoundary)*n_strong_boundaries);
823
824
825
  problem->n_strong_boundaries = n_strong_boundaries;
  for (int i = 0; i < n_strong_boundaries; ++i) {
    problem->strong_boundaries[i].tag = strdup(strong_boundaries[i].tag);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
826
    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
827
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
Jonathan Lambrechts's avatar
weak bc    
Jonathan Lambrechts committed
828
    problem->strong_boundaries[i].row = strong_boundaries[i].row;
829
  }
830
  problem->compacity = (double*)malloc(mesh->n_nodes*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
831
832
  for(int i = 0; i < mesh->n_nodes; ++i){
    problem->compacity[i] = 0.;
Matthieu Constant's avatar
Matthieu Constant committed
833
  }
834
  problem->solution = (double*)malloc(mesh->n_nodes*n_fields*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
835
  for (int i = 0; i < problem->mesh->n_nodes*n_fields; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
836
    problem->solution[i] = 0.;
837
  }
838
  problem->node_volume = NULL;
Matthieu Constant's avatar
Matthieu Constant committed
839
  fluid_problem_compute_node_volume(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
840
  problem->linear_system = NULL;
841
#ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
842
    hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
843
#else
844
    hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
845
#endif
846
  problem->n_particles = 0;
847
848
849
850
851
852
853
  problem->particle_uvw = NULL;
  problem->particle_element_id = NULL;
  problem->particle_position = NULL;
  problem->particle_mass = NULL;
  problem->particle_force = NULL;
  problem->particle_volume = NULL;
  problem->particle_velocity = NULL;
854
855
856
857
858
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
859
860
  free(problem->mu);
  free(problem->rho);
861
  free(problem->solution);
862
  free(problem->compacity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
863
  hxtLinearSystemDelete(&problem->linear_system);
864
865
866
867
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
868
  free(problem->particle_force);
869
870
871
872
873
  free(problem->particle_volume);
  free(problem->particle_velocity);
  for (int i = 0; i < problem->n_strong_boundaries; ++i)
    free(problem->strong_boundaries[i].tag);
  free(problem->strong_boundaries);
874
  mesh_tree_free(problem->mesh_tree);
875
876
877
  mesh_free(problem->mesh);
}

Matthieu Constant's avatar
Matthieu Constant committed
878
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
879
{
Matthieu Constant's avatar
Matthieu Constant committed
880
  double ngradM = 0.;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
881
882
  Mesh *mesh = problem->mesh;
  double *solution = problem->solution;
883
884
  double *new_mesh_size = (double*)malloc(sizeof(double)*mesh->n_nodes);
  double *num_lc = (double*)malloc(sizeof(double)*mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
885
886
887
888
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] = 0.;
    num_lc[i] = 0.;
  }
Matthieu Constant's avatar
Matthieu Constant committed
889
890
891
892
893
894
895
  double gradmax = 0.;
  double gradPmax = 0.;
  double gradmin = 1e8;
  double gradPmin = 1e8;
  double gradM = 0.;
  double gradPM = 0.;
  double volT = 0.;
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
896
  int n_fields = fluid_problem_n_fields(problem);
897
  double *compacity = problem->compacity;
Matthieu Constant's avatar
Matthieu Constant committed
898
899
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const unsigned int *el = &mesh->elements[iel*N_N];
900
    double C[N_N],U[D][N_N], P[N_N];
Matthieu Constant's avatar
Matthieu Constant committed
901
    for (int i = 0; i < N_N; ++i){
902
      C[i] = 1-compacity[el[i]];
903
904
      P[i] = solution[el[i]*n_fields+D];
      for (int j=0; j<D; ++j) {
Matthieu Constant's avatar
Matthieu Constant committed
905
906
907
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
908
909
910
    double dxdxi[D][D],dxidx[D][D];
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
Matthieu Constant's avatar
Matthieu Constant committed
911
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
912
#if D == 2
Matthieu Constant's avatar
Matthieu Constant committed
913
914
915
916
917
    const double vol = detDD(dxdxi)/2;
#else        
    const double vol = detDD(dxdxi)/6;
#endif    
    invDD(dxdxi, dxidx);
918
    double dphi[N_SF][D];
Matthieu Constant's avatar
Matthieu Constant committed
919
920
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
921
922
    for (int i=0; i<D; ++i){
      for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
923
924
925
926
927
928
929
930
931
932
933
934
        double gradU = 0;
        for (int k=0; k<N_SF; ++k){
          gradU += dphi[k][j]*U[i][k]/C[k];
        }
        ngrad2 += gradU*gradU;
      }
    }
    gradM += sqrt(ngrad2)*vol;
    volT += vol;
    gradmax = fmax(gradmax,ngrad2);
    gradmin = fmin(gradmin,ngrad2);
    ngrad2 = 0;
935
    for (int j=0; j<D; ++j){
Matthieu Constant's avatar
Matthieu Constant committed
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
      double gradP = 0;
      for (int k=0; k<N_SF; ++k){
        gradP += dphi[k][j]*P[k];
      }
      ngrad2 += gradP*gradP;
    }
    gradPM += sqrt(ngrad2)*vol;
    gradPmax = fmax(gradPmax,ngrad2);
    gradPmin = fmin(gradPmin,ngrad2);
  }
  gradPM /= volT;
  gradM /= volT;
  gradPmin = 20*sqrt(gradPmin)/gradPM;
  gradPmax = .1*sqrt(gradPmax)/gradPM;
  gradmin = 20*sqrt(gradmin)/gradM;
  gradmax = .1*sqrt(gradmax)/gradM;
  printf("gradPmax=%e\n",gradPmax);
  printf("gradPmin=%e\n",gradPmin);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
954
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
955
    const unsigned int *el = &mesh->elements[iel*N_N];
956
    double C[N_N],U[D][N_N], P[N_N];
957
    for (int i = 0; i < N_N; ++i){
958
      C[i] = 1-compacity[el[i]];
959
960
      P[i] = solution[el[i]*n_fields+D];
      for (int j=0; j<D; ++j) {
961
962
963
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
964
965
966
    double dxdxi[D][D],dxidx[D][D];
    for (int i=0; i<D; ++i)
      for (int j=0; j<D; ++j)
967
968
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
    invDD(dxdxi, dxidx);
969
    double dphi[N_SF][D];
970
971
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
972
973
    for (int i=0; i<D; ++i){
      for (int j=0; j<D; ++j){
974
975
976
977
978
979
980
        double gradU = 0;
        for (int k=0; k<N_SF; ++k){
          gradU += dphi[k][j]*U[i][k]/C[k];
        }
        ngrad2 += gradU*gradU;
      }
    }
Matthieu Constant's avatar
Matthieu Constant committed
981
/*    double ngrad = pow(ngrad2, 0.25); */
982
/*    for (int j=0; j<D; ++j){*/
Matthieu Constant's avatar
Matthieu Constant committed
983
984
985
986
987
988
/*      double gradP = 0;*/
/*      for (int k=0; k<N_SF; ++k){*/
/*        gradP += dphi[k][j]*P[k];*/
/*      }*/
/*      ngrad2 += gradP*gradP;*/
/*    }*/
Matthieu Constant's avatar
Matthieu Constant committed
989
    double ngrad = sqrt(ngrad2)/gradM;
Matthieu Constant's avatar
Matthieu Constant committed
990
    double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
Matthieu Constant's avatar
Matthieu Constant committed
991