fluid_problem.c 30.1 KB
Newer Older
1
#include "tools.h"
2
3
#include <stdlib.h>
#include <stdio.h>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
#include <float.h>
5
6
#include <math.h>
#include "fluid_problem.h"
7
#include "mesh_find.h"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
#include "mbxml.h"
9

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
#define n_fields (DIMENSION+1)
11

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
12
#define newton_max_it 10
13
#define newton_atol 2.5e-6
14
#define newton_rtol 1e-5
15

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
16
#if DIMENSION==2
17
18
19

#define N_SF 3
#define N_N 3
20
21
22
#define N_QP 3
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};
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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
50
#else
51
52
53

#define N_SF 4
#define N_N 4
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
54
55
56
57
58
59
60
61
62
#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};
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
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
99
100
#endif

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
102
#if DIMENSION==2
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
103
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
104
  double d = 2*sqrt(vol/M_PI);
105
106
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
Matthieu Constant's avatar
Matthieu Constant committed
107
  double Re_O_u = sqrt(d)*c/mu*rho;
108
109
110
  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
111
  return f*Cd_u*rho/2*d;
112
}
113
#else
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
114
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
115
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
116
  double d = 2*cbrt(3./4.*vol/M_PI);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
117
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
118
119
120
121
122
123
  //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
124
  return f*Cd_u*rho/2*(M_PI*r*r);
125
126
}
#endif
127

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
128
void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
129
130
131
132
133
134
135
136
137
138
139
140
141
  for (int i = 0; i < problem->n_particles*DIMENSION; ++i) {
    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) {
  double *porosity = problem->porosity;
  Mesh *mesh = problem->mesh;
  for (int i = 0; i < problem->n_particles; ++i) {
    int iel = problem->particle_element_id[i];
    if(iel < 0){
      continue;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
142
143
144
145
146
    double *xi = problem->particle_uvw + DIMENSION*i;
    double phi[N_SF];
    shape_functions(xi,phi);
    uint32_t *el=problem->mesh->elements+iel*N_SF;
    double *se = problem->solution_explicit;
147
    double c=0, vf[DIMENSION]={0}, vfe[DIMENSION]={0}, gradp[DIMENSION]={0};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
149
150
151
152
153
154
    double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
    for (int k=0; k<DIMENSION; ++k)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[k][j] = mesh->x[el[j+1]*3+k]-mesh->x[el[0]*3+k];
    invDD(dxdxi,dxidx);
    double dphi[N_SF][DIMENSION];
    grad_shape_functions(dxidx, dphi);
155
    double *si = problem->solution;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
156
157
158
    for (int iphi=0; iphi < N_SF; ++iphi) {
      c += porosity[el[iphi]]*phi[iphi];
      for (int j=0; j < DIMENSION; ++j) {
159
        vf[j]  += si[el[iphi]*n_fields+j]*phi[iphi];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
160
        vfe[j] += se[el[iphi]*n_fields+j]*phi[iphi];
161
        gradp[j] += dphi[iphi][j]*si[el[iphi]*n_fields+DIMENSION];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
162
163
164
165
166
167
168
      }
    }
    double du[DIMENSION], due[DIMENSION];
    for (int j = 0; j < DIMENSION; ++j) {
      du[j] = problem->particle_velocity[i*DIMENSION+j]-vf[j]/c;
      due[j] = problem->particle_velocity[i*DIMENSION+j]-vfe[j]/c;
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
169
170
    double vol = problem->particle_volume[i];
    double mass = problem->particle_mass[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
171
172
    double rhop = mass/vol;
    double drho = rhop -problem->rho;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
173
174
175
176
177
178
179
180
    double gamma = particle_drag_coeff(due,problem->mu,problem->rho,vol,c);
    double g = problem->g*drho/rhop;
    double fcoeff = mass/(mass+dt*gamma);
    
    for (int d = 0; d < DIMENSION; ++d)
      problem->particle_force[i*DIMENSION+d] = fcoeff*(-gamma*du[d]-vol*gradp[d]);
    problem->particle_force[i*DIMENSION+1] += fcoeff*g*mass;
    
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
181
182
    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
183
184
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int d = 0; d < DIMENSION; ++d)
185
186
187
        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
188
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
189
190
    for (int iphi = 0; iphi < N_SF; ++iphi) {
      for (int jphi = 0; jphi < N_SF; ++jphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
191
192
193
        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
194
195
        double f = fcoeff*phi[iphi];
        for (int d = 0; d < DIMENSION; ++d){
196
197
          LOCAL_MATRIX(d,d) -= -f/c*phi[jphi]*gamma;
          LOCAL_MATRIX(d,DIMENSION) -= -f*gamma*dt/mass*vol*dphi[jphi][d];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
198
        }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
199
200
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
201
202
203
204
205
206
207
208
209
210
211
212
  }
}

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;
  const double *solution = problem->solution;
  const double *porosity = problem->porosity;
  const double *old_porosity = problem->old_porosity;
  const double mu = problem->mu;
  const double rho = problem->rho;
Matthieu Constant's avatar
Matthieu Constant committed
213
  const double *epsilon = problem->epsilon;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
214
215
216
217
218
219
220
221
222
  size_t local_size = N_SF*n_fields;
  size_t N_E = mesh->n_elements;
  double *all_local_vector = malloc(N_E*local_size*sizeof(double));
  double *all_local_matrix = 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;
  compute_node_force_implicit(problem, dt, all_local_vector, all_local_matrix);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
223
  for (int iel=0; iel < mesh->n_elements; ++iel) {
224
225
226
227
228
229
230
    const unsigned int *el = &mesh->elements[iel*N_N];
    double dxdxi[DIMENSION][DIMENSION], dxidx[DIMENSION][DIMENSION], dphi[N_N][DIMENSION];
    for (int i = 0; i < DIMENSION; ++i)
      for (int j = 0; j < DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
    const double detj = invDD(dxdxi, dxidx);
    grad_shape_functions(dxidx, dphi);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
231
232
    double *local_vector = &all_local_vector[N_SF*n_fields*iel];
    double *local_matrix = &all_local_matrix[N_SF*N_SF*n_fields*n_fields*iel];
233
234
235
236
237
238
239
240
241
    double dc[DIMENSION]={0}, du[DIMENSION][DIMENSION]={{0}}, dp[DIMENSION]={0};
    for (int i = 0; i< N_SF; ++i) {
      for (int j = 0; j < DIMENSION; ++j) {
        dc[j] += dphi[i][j]*porosity[el[i]];
        dp[j] += dphi[i][j]*solution[el[i]*n_fields+DIMENSION];
        for (int k = 0; k < DIMENSION; ++k)
          du[k][j] += dphi[i][j]*solution[el[i]*n_fields+k];
      }
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
242
    for (int iqp = 0; iqp< N_QP; ++iqp) {
243
244
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
245
      double u[DIMENSION]={0}, dudt[DIMENSION]={0}, p=0, c=0, dcdt=0;
246
247
248
249
250
251
252
253
254
      for (int i = 0; i < N_SF; ++i) {
        for (int j = 0; j < DIMENSION; ++j){
          u[j] += solution[el[i]*n_fields+j]*phi[i];
          dudt[j] += (solution[el[i]*n_fields+j]-solution_old[el[i]*n_fields+j])/dt*phi[i];
        }
        p += solution[el[i]*n_fields+DIMENSION]*phi[i];
        c += porosity[el[i]]*phi[i];
        dcdt += (porosity[el[i]]-old_porosity[el[i]])/dt*phi[i];
      }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
255
      const double jw = QW[iqp]*detj;
256
      for (int iphi = 0; iphi < N_SF; ++iphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
257
        double phii = phi[iphi];
258
259
260
261
262
263
264
265
266
        const double *dphii = dphi[iphi];
        double tau[DIMENSION][DIMENSION];
        for (int i = 0; i < DIMENSION; ++i)
          for (int j = 0; j < DIMENSION; ++j)
            tau[i][j] = du[i][j]-u[i]*dc[j]/c;
        for (int j = 0; j < DIMENSION; ++j) {
          double utau = 0;
          double dphisigma=0;
          for (int k = 0; k < DIMENSION; ++k) {
267
268
            //utau += u[k]*dphii[k]*u[j]/c;
            utau += u[k]*tau[j][k];
269
270
271
272
            dphisigma += 2*mu*dphii[k]*0.5*(tau[k][j]+tau[j][k]);
          }
          local_vector[iphi+N_SF*j] += jw*(
              dphisigma
273
274
275
              //+rho*phii*dudt[j]-rho*utau
              +rho*phii*(dudt[j]-u[j]/c*dcdt+utau/c)
              -p*dphii[j]
276
277
278
279
280
281
282
              );
        }
        double dphidp = 0, divu = 0;
        for (int k = 0; k < DIMENSION; ++k){
          dphidp += dphii[k]*dp[k];
          divu += du[k][k];
        }
Matthieu Constant's avatar
Matthieu Constant committed
283
        local_vector[iphi+N_SF*DIMENSION] += jw*(epsilon[el[iphi]]*dphidp+phii*(divu+dcdt));
284
        for (int jphi = 0; jphi < N_SF; ++jphi) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
285
          double phij = phi[jphi];
286
          const double *dphij = dphi[jphi];
287
288
289
290
291
          int P = DIMENSION;
          //int N2 = N_SF*N_SF;
          //int IDX = iphi*N2+jphi;
          int N2 = n_fields*N_SF;
          int IDX = (iphi*N2+jphi)*n_fields;
292
293
294
295
296
297
298
299
#define LOCAL_MATRIX(a,b) local_matrix[IDX+N2*a+b]
          double dphiidphij = 0;
          double dtau[DIMENSION];
          double udtau = 0;
          double dphiidtau = 0;
          for (int j = 0; j < DIMENSION; ++j) {
            dtau[j] = dphij[j]-phij*dc[j]/c;
            dphiidphij += dphii[j]*dphij[j];
300
            udtau += u[j]*dtau[j];
301
302
303
            dphiidtau += dphii[j]*dtau[j];
          }
          for (int j = 0; j < DIMENSION; ++j) {
304
            int U = j;
305
306
307
            LOCAL_MATRIX(U,P) += -jw*dphii[j]*phij;
            LOCAL_MATRIX(P,U) +=  jw*phii*dphij[j];
            for (int k = 0; k < DIMENSION; ++k) {
308
              int V = k;
309
310
              //LOCAL_MATRIX(U,V) += -jw*rho*u[j]*dphii[k]*phij/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
              LOCAL_MATRIX(U,V) += jw*rho*phii*phij*tau[j][k]/c+jw*2*mu*dphii[k]*dtau[j]*0.5;
311
            }
312
313
            //LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*(phii*phij/dt-phij*dphii[j]*u[j]/c);
            LOCAL_MATRIX(U,U) += jw*mu*2*0.5*dphiidtau + jw*rho*phii*(phij/dt-phij/c*dcdt+udtau/c);
314
          }
Matthieu Constant's avatar
Matthieu Constant committed
315
          LOCAL_MATRIX(P,P) +=  jw*epsilon[el[iphi]]*dphiidphij;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
316
317
318
319
320
        }
      }
    }
    hxtLinearSystemAddToMatrix(lsys,iel,iel,local_matrix);
    hxtLinearSystemAddToRhs(lsys,rhs,iel,local_vector);
321
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
322
323
  free(all_local_matrix);
  free(all_local_vector);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
324

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
325
326
327
328
329
330
331
332
333
334
  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){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
335
336
      hxtLinearSystemSetMatrixRowIdentity(lsys, mesh->phys_nodes[iphys][i], bnd->field);
      hxtLinearSystemSetRhsEntry(lsys, rhs,mesh->phys_nodes[iphys][i], bnd->field, 0.);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
337
338
    }
  }
339
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
340
341


342

343
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
344
{
345
346
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
347
348
349
350
351
352
353
354
    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);
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
355
    double *x = malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*DIMENSION);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
356
357
358
    double *v = malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
359
360
      for (int k = 0; k < DIMENSION; ++k)
        x[i*DIMENSION+k] = mesh->x[j*3+k];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
361
362
    }
    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
363
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
364
365
366
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
    free(x);
    free(v);
367
368
369
  }
}

370
int fluid_problem_implicit_euler(FluidProblem *problem, double dt)
371
{
372
373
  static double totaltime = 0;
  struct timespec tic, toc;
374
  const Mesh *mesh = problem->mesh;
375
  apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
376
  double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
377
  double *rhs = malloc(mesh->n_nodes*n_fields*sizeof(double));
378
379
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
380
  problem->solution_explicit = solution_old;
381
382
  double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
  double norm0 = 0;
383
384
  int i;
  for (i=0; i<newton_max_it; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
385
386
387
388
389
390
    double norm = 0;
    hxtLinearSystemZeroMatrix(problem->linear_system);
    for (int i=0; i<mesh->n_nodes*n_fields; ++i)
      rhs[i] = 0;
    fluid_problem_assemble_system(problem, rhs, solution_old, dt);
    hxtLinearSystemGetRhsNorm(problem->linear_system,rhs,&norm);
391
392
393
394
395
396
397
    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;
398
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
399
    hxtLinearSystemSolve(problem->linear_system, rhs, dx);
400
401
  clock_get(&toc);
  totaltime += clock_diff(&tic, &toc);
402
403
404
405
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
406
  printf("total solve time : %g\n", totaltime);
407
  free(dx);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
408
  free(rhs);
409
  free(solution_old);
410
  return i;
411
412
}

Matthieu Constant's avatar
Matthieu Constant committed
413
414
415
416
417
418
419
static void fluid_problem_compute_node_volume(FluidProblem *problem) {
  free(problem->node_volume);
  Mesh *mesh = problem->mesh;
  problem->node_volume = malloc(sizeof(double)*mesh->n_nodes);
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->node_volume[i] = 0;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
420
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
421
422
423
424
425
    double dxdxi[DIMENSION][DIMENSION];
    const int *el = &mesh->elements[iel*N_N];
    for (int i=0; i<DIMENSION; ++i)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
426
#if DIMENSION == 2
427
    const double vol = detDD(dxdxi)/6;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
428
#else
429
    const double vol = detDD(dxdxi)/24;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
430
#endif
431
432
    for (int i = 0; i < N_N; ++i)
      problem->node_volume[el[i]] += vol;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
433
434
435
436
437
438
439
440
441
442
443
444
  }
}

static void fluid_problem_compute_porosity(FluidProblem *problem)
{
  Mesh *mesh = problem->mesh;
  double *volume = problem->particle_volume;
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->porosity[i] = 0;
  }
  for (int i = 0; i < problem->n_particles; ++i) {
    if(problem->particle_element_id[i] == -1) continue;
445
446
    double sf[N_SF];
    shape_functions(&problem->particle_uvw[i*DIMENSION],sf);
447
    const uint32_t *el = &mesh->elements[problem->particle_element_id[i]*N_N];
448
449
    for (int j=0; j<N_N;++j)
      problem->porosity[el[j]] += sf[j]*volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
450
451
452
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->porosity[i] = (1-problem->porosity[i]/problem->node_volume[i]);
Matthieu Constant's avatar
Matthieu Constant committed
453
454
455
  }
}

Matthieu Constant's avatar
Matthieu Constant committed
456
FluidProblem *fluid_problem_new(const char *mesh_file_name, double g, double mu, double rho, double alpha, int n_strong_boundaries, StrongBoundary *strong_boundaries, int notComputeEpsilon)
457
{
458
459
460
461
462
  static int initialized = 0;
  if (!initialized) {
    hxtInitializeLinearSystems(0, NULL);
    initialized = 1;
#ifdef HXT_HAVE_PETSC
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
463
    hxtPETScInsertOptions("-pc_type ilu -ksp_max_it 30", "fluid");
464
465
#endif
  }
466
  FluidProblem *problem = malloc(sizeof(FluidProblem));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
467
468
469
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
470
  problem->mesh_tree = mesh_tree_create(mesh);
471
  problem->rho = rho;
Matthieu Constant's avatar
Matthieu Constant committed
472

473
  problem->mu = mu;
474
  problem->g = g;
475
  problem->mesh = mesh;
476
477
478
479
  problem->strong_boundaries = malloc(sizeof(StrongBoundary)*n_strong_boundaries);
  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
480
    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
481
482
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
483
  problem->porosity = malloc(mesh->n_nodes*sizeof(double));
Matthieu Constant's avatar
Matthieu Constant committed
484
  problem->epsilon = malloc(mesh->n_nodes*sizeof(double));
485
  problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
486
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
Matthieu Constant's avatar
Matthieu Constant committed
487
488
489
490
491
  for (int i = 0; i < problem->mesh->n_nodes; ++i)
  {
    problem->porosity[i] = 1.;
    problem->old_porosity[i] = 1.;
  }
492
  // begin to remove
493
  for (int i = 0; i < problem->mesh->n_nodes*N_N; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
494
    problem->solution[i] = 0.;
495
  }
Matthieu Constant's avatar
Matthieu Constant committed
496
497
  problem->node_volume = malloc(0);
  fluid_problem_compute_node_volume(problem);
Matthieu Constant's avatar
Matthieu Constant committed
498
499
500
#if DIMENSION == 2
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
501
  }
Matthieu Constant's avatar
Matthieu Constant committed
502
503
504
505
506
507
508
509
510
511
#else
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
  }
#endif
  if(notComputeEpsilon){
    for (int i = 0; i < problem->mesh->n_nodes; ++i){
      problem->epsilon[i] = alpha;
    }
  }
512
  // end to remove
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
513
  problem->linear_system = NULL;
514
#ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
515
    hxtLinearSystemCreatePETSc(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements,"fluid");
516
#else
517
    hxtLinearSystemCreateLU(&problem->linear_system,mesh->n_elements,N_N,n_fields,mesh->elements);
518
#endif
519
520
521
522
523
  problem->n_particles = 0;
  problem->particle_uvw = malloc(0);
  problem->particle_element_id = malloc(0);
  problem->particle_position = malloc(0);
  problem->particle_mass = malloc(0);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
524
  problem->particle_force = malloc(0);
525
526
  problem->particle_volume = malloc(0);
  problem->particle_velocity = malloc(0);
527
528
529
530
531
532
533
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
  free(problem->porosity);
Matthieu Constant's avatar
Matthieu Constant committed
534
  free(problem->epsilon);
535
  free(problem->old_porosity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
536
  hxtLinearSystemDelete(&problem->linear_system);
537
538
539
540
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
541
  free(problem->particle_force);
542
543
544
545
546
  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);
547
  mesh_tree_free(problem->mesh_tree);
548
549
550
  mesh_free(problem->mesh);
}

Matthieu Constant's avatar
Matthieu Constant committed
551
void fluid_problem_adapt_gen_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
552
{
Matthieu Constant's avatar
Matthieu Constant committed
553
  double ngradM = 0.;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
554
555
556
  Mesh *mesh = problem->mesh;
  double *solution = problem->solution;
  double *new_mesh_size = malloc(sizeof(double)*mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
557
558
559
560
561
  double *num_lc = malloc(sizeof(double)*mesh->n_nodes);
  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
562
563
564
565
566
567
568
  double gradmax = 0.;
  double gradPmax = 0.;
  double gradmin = 1e8;
  double gradPmin = 1e8;
  double gradM = 0.;
  double gradPM = 0.;
  double volT = 0.;
Matthieu Constant's avatar
Matthieu Constant committed
569
  double *porosity = problem->porosity;
Matthieu Constant's avatar
Matthieu Constant committed
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    const unsigned int *el = &mesh->elements[iel*N_N];
    double C[N_N],U[DIMENSION][N_N], P[N_N];
    for (int i = 0; i < N_N; ++i){
      C[i] = porosity[el[i]];
      P[i] = solution[el[i]*n_fields+DIMENSION];
      for (int j=0; j<DIMENSION; ++j) {
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
    double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
    for (int i=0; i<DIMENSION; ++i)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
#if DIMENSION == 2
    const double vol = detDD(dxdxi)/2;
#else        
    const double vol = detDD(dxdxi)/6;
#endif    
    invDD(dxdxi, dxidx);
    double dphi[N_SF][DIMENSION];
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
    for (int i=0; i<DIMENSION; ++i){
      for (int j=0; j<DIMENSION; ++j){
        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;
    for (int j=0; j<DIMENSION; ++j){
      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
626
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
627
    const unsigned int *el = &mesh->elements[iel*N_N];
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
    double C[N_N],U[DIMENSION][N_N], P[N_N];
    for (int i = 0; i < N_N; ++i){
      C[i] = porosity[el[i]];
      P[i] = solution[el[i]*n_fields+DIMENSION];
      for (int j=0; j<DIMENSION; ++j) {
        U[j][i] = solution[el[i]*n_fields+j];
      }
    }
    double dxdxi[DIMENSION][DIMENSION],dxidx[DIMENSION][DIMENSION];
    for (int i=0; i<DIMENSION; ++i)
      for (int j=0; j<DIMENSION; ++j)
        dxdxi[i][j] = mesh->x[el[j+1]*3+i]-mesh->x[el[0]*3+i];
    invDD(dxdxi, dxidx);
    double dphi[N_SF][DIMENSION];
    grad_shape_functions(dxidx, dphi);
    double ngrad2 = 0;
Matthieu Constant's avatar
Matthieu Constant committed
644
    for (int i=0; i<DIMENSION; ++i){
645
646
647
648
649
650
651
652
      for (int j=0; j<DIMENSION; ++j){
        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
653
654
655
656
657
658
659
660
/*    double ngrad = pow(ngrad2, 0.25); */
/*    for (int j=0; j<DIMENSION; ++j){*/
/*      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
661
    double ngrad = sqrt(ngrad2)/gradM;
Matthieu Constant's avatar
Matthieu Constant committed
662
    double lc = lcmin /fmin(1,fmax(ngrad/(gradmax),gradmin/(gradmax)));
Matthieu Constant's avatar
Matthieu Constant committed
663
    
Matthieu Constant's avatar
Matthieu Constant committed
664
    ngrad2 = 0.; 
Matthieu Constant's avatar
Matthieu Constant committed
665
666
667
668
669
670
671
    for (int j=0; j<DIMENSION; ++j){
      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
672
673
    ngrad = sqrt(ngrad2)/gradPM;
    lc = fmin(fmin(lcmin /fmin(1,fmax(ngrad/(gradPmax),gradPmin/(gradPmax))),lc),lcmax);
Matthieu Constant's avatar
Matthieu Constant committed
674
675
    
    
Matthieu Constant's avatar
Matthieu Constant committed
676
677
678
679
680
    for (int j = 0; j < N_N; ++j){
      new_mesh_size[el[j]] += lc;
      num_lc[el[j]] += 1.;
    }
  }
Matthieu Constant's avatar
Matthieu Constant committed
681
  double nOfEl = 0.0;
Matthieu Constant's avatar
Matthieu Constant committed
682
#if DIMENSION==2
Matthieu Constant's avatar
Matthieu Constant committed
683
684
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] /= num_lc[i];
Matthieu Constant's avatar
Matthieu Constant committed
685
686
687
688
689
690
691
    double nOfElByNode = 4*problem->node_volume[i]/(new_mesh_size[i]*new_mesh_size[i]*sqrt(3));
    nOfEl += nOfElByNode;
  }
  double ratioEl = n_el/nOfEl;
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] = fmax(new_mesh_size[i]/sqrt(ratioEl),lcmin);
  }
Matthieu Constant's avatar
Matthieu Constant committed
692
#else
Matthieu Constant's avatar
Matthieu Constant committed
693
  for (int i = 0; i < mesh->n_nodes; ++i){
Matthieu Constant's avatar
Matthieu Constant committed
694
695
    new_mesh_size[i] /= num_lc[i];
    double nOfElByNode = 12*problem->node_volume[i]/(new_mesh_size[i]*new_mesh_size[i]*new_mesh_size[i]*sqrt(2));
Matthieu Constant's avatar
Matthieu Constant committed
696
    nOfEl += nOfElByNode;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
697
  }
Matthieu Constant's avatar
Matthieu Constant committed
698
  double ratioEl = n_el/nOfEl;
Matthieu Constant's avatar
Matthieu Constant committed
699
  printf("RatioEl = %e\n",ratioEl);
Matthieu Constant's avatar
Matthieu Constant committed
700
701
702
703
  for (int i = 0; i < mesh->n_nodes; ++i){
    new_mesh_size[i] = fmax(new_mesh_size[i]/cbrt(ratioEl),lcmin);
  }
#endif
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
704
  FILE *f = fopen("adapt/lc.pos", "w");
705
706
  if(!f)
    printf("cannot open file adapt/lc.pos for writing\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
707
  fprintf(f, "View \"lc\" {\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
708
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
709
710
711
712
713
714
    unsigned int *el = problem->mesh->elements+iel*N_N;
    fprintf(f, DIMENSION == 2 ? "ST(" : "SS(");
    for (int i = 0; i < N_N; ++i)
      fprintf(f, "%.8g, %.8g, %.8g%s", mesh->x[el[i]*3+0], mesh->x[el[i]*3+1],mesh->x[el[i]*3+2], i != N_N-1 ? "," : "){");
    for (int i = 0; i < N_N; ++i)
      fprintf(f, "%.8g%s", new_mesh_size[el[i]],i != N_N-1?",":"};\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
715
  }
716
  fprintf(f,"};\n");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
717
718
  fclose(f);
  free(new_mesh_size);
Matthieu Constant's avatar
Matthieu Constant committed
719
  free(num_lc);
720
  system("gmsh -"xstr(DIMENSION)" adapt/mesh.geo");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
721
722
723
}


Matthieu Constant's avatar
Matthieu Constant committed
724
void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double alpha, int notComputeEpsilon)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
725
{
Matthieu Constant's avatar
Matthieu Constant committed
726
727
728
  struct timespec tic, toc;
  fluid_problem_adapt_gen_mesh(problem, lcmax, lcmin, n_el);
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
729
  Mesh *new_mesh = mesh_load_msh("adapt/mesh.msh");
Matthieu Constant's avatar
Matthieu Constant committed
730
731
  clock_get(&toc);
  printf("Time mesh_load : %.6e \n",clock_diff(&tic,&toc));
732
  double *new_solution = malloc(sizeof(double)*new_mesh->n_nodes*N_N);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
733
  double *new_xi = malloc(sizeof(double)*new_mesh->n_nodes*DIMENSION);
Matthieu Constant's avatar
Matthieu Constant committed
734
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
735
  int *new_eid = malloc(sizeof(int)*new_mesh->n_nodes);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
736
737
  for (int i = 0; i < new_mesh->n_nodes; ++i)
    new_eid[i] = -1;
Matthieu Constant's avatar
Matthieu Constant committed
738
739
740
  clock_get(&toc);
  printf("Time new_eid : %.6e \n",clock_diff(&tic,&toc));
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
741
742
743
744
  double *newx = malloc(sizeof(double)*DIMENSION*new_mesh->n_nodes);
  for (int i = 0;i < new_mesh->n_nodes;++i) 
    for (int k = 0; k < DIMENSION; ++k)
      newx[i*DIMENSION+k] = new_mesh->x[i*3+k];
Matthieu Constant's avatar
Matthieu Constant committed
745
746
747
  clock_get(&toc);
  printf("Time newx : %.6e \n",clock_diff(&tic,&toc));
  clock_get(&tic);
748
  mesh_tree_particle_to_mesh(problem->mesh_tree,new_mesh->n_nodes, newx, new_eid, new_xi);
Matthieu Constant's avatar
Matthieu Constant committed
749
750
  clock_get(&toc);
  printf("Time find new point in old mesh : %.6e \n",clock_diff(&tic,&toc));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
751
  for (int i = 0; i < new_mesh->n_nodes; ++i) {
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
752
753
754
755
756
    unsigned int *el = problem->mesh->elements+new_eid[i]*N_N;
    if(new_eid[i] < 0) {
      printf("new mesh point not found in old mesh\n");
      exit(1);
    }
757
758
759
760
761
    double phi[N_SF];
    shape_functions(new_xi+i*DIMENSION, phi);
    for (int j=0; j<n_fields; ++j) {
      new_solution[i*n_fields+j] = 0;
      for (int k = 0; k < N_SF; ++k)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
762
        new_solution[i*n_fields+j] += problem->solution[el[k]*n_fields+j]*phi[k];
763
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
764
765
766
767
768
769
  }
  free(new_eid);
  free(new_xi);
  free(newx);
  free(problem->solution);
  problem->solution = new_solution;
Matthieu Constant's avatar
Matthieu Constant committed
770
771
  free(problem->epsilon);
  problem->epsilon = malloc(new_mesh->n_nodes*sizeof(double));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
772
773
774
775
  free(problem->porosity);
  problem->porosity = malloc(sizeof(double)*new_mesh->n_nodes);
  free(problem->old_porosity);
  problem->old_porosity = malloc(sizeof(double)*new_mesh->n_nodes);
Matthieu Constant's avatar
Matthieu Constant committed
776
777
778
779
780
  for (int i = 0; i < new_mesh->n_nodes; ++i)
  {
    problem->porosity[i] = 1.;
    problem->old_porosity[i] = 1.;
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
781
782
  mesh_free(problem->mesh);
  problem->mesh = new_mesh;
783
784
  mesh_tree_free(problem->mesh_tree);
  problem->mesh_tree = mesh_tree_create(problem->mesh);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
785
786
  for (int i = 0; i < problem->n_particles; ++i)
    problem->particle_element_id[i] = -1;
787
  mesh_tree_particle_to_mesh(problem->mesh_tree, problem->n_particles, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
788
  fluid_problem_compute_node_volume(problem);
Matthieu Constant's avatar
Matthieu Constant committed
789
790
791
#if DIMENSION == 2
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->epsilon[i] = alpha*sqrt(problem->node_volume[i])*sqrt(problem->node_volume[i])/(problem->mu/problem->rho);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
792
  }
Matthieu Constant's avatar
Matthieu Constant committed
793
794
795
796
797
798
799
#else
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->epsilon[i] = alpha*cbrt(problem->node_volume[i])*cbrt(problem->node_volume[i])/(problem->mu/problem->rho);
  }
#endif
  if(notComputeEpsilon){
    for (int i = 0; i < problem->mesh->n_nodes; ++i){
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
800
801
      problem->epsilon[i] = alpha;
    }
Matthieu Constant's avatar
Matthieu Constant committed
802
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
803
  fluid_problem_compute_porosity(problem);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
804
805
  hxtLinearSystemDelete(&problem->linear_system);
  #ifdef HXT_HAVE_PETSC
Matthieu Constant's avatar
Matthieu Constant committed
806
    hxtLinearSystemCreatePETSc(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements,"fluid");
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
807
  #else
808
    hxtLinearSystemCreateLU(&problem->linear_system,new_mesh->n_elements,N_N,n_fields,new_mesh->elements);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
809
  #endif 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
810
811
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
812
813

void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
814
{
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
815
816
817
  
  struct timespec tic, toc;
  
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
818
819
820
821
822
823
  if(problem->n_particles != n) {
    problem->n_particles = n;
    free(problem->particle_uvw);
    free(problem->particle_element_id);
    free(problem->particle_position);
    free(problem->particle_mass);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
824
    free(problem->particle_force);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
825
826
    free(problem->particle_volume);
    free(problem->particle_velocity);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
827
    problem->particle_uvw = malloc(sizeof(double)*DIMENSION*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
828
    problem->particle_element_id = malloc(sizeof(int)*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
829
830
831
832
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=-1;
    }
    problem->particle_position = malloc(sizeof(double)*DIMENSION*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
833
    problem->particle_mass = malloc(sizeof(double)*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
834
    problem->particle_force = malloc(sizeof(double)*n*DIMENSION);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
835
    problem->particle_volume = malloc(sizeof(double)*n);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
836
837
838
839
840
841
    problem->particle_velocity = malloc(sizeof(double)*n*DIMENSION);
  }
  if (elid) {
    for (int i = 0; i < n; ++i) {
      problem->particle_element_id[i]=elid[i];
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
842
  }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
843
  clock_get(&tic);
844
  mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
845
846
  clock_get(&toc);
  printf("Time mesh_particles_to_mesh : %.6e \n",clock_diff(&tic,&toc));
Matthieu Constant's avatar
Matthieu Constant committed
847
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
848
849
850
  for (int i = 0; i < n; ++i) {
    problem->particle_mass[i] = mass[i];
    problem->particle_volume[i] = volume[i];
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
851
852
853
854
    for (int k = 0; k < DIMENSION; ++k) {
      problem->particle_position[i*DIMENSION+k] = position[i*DIMENSION+k];
      problem->particle_velocity[i*DIMENSION+k] = velocity[i*DIMENSION+k];
    }
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
855
  }
Matthieu Constant's avatar
Matthieu Constant committed
856
857
  clock_get(&toc);
  printf("Time initialize particle_fluid : %.6e \n",clock_diff(&tic,&toc));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
858
859
860
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->old_porosity[i] = problem->porosity[i];
  }
Matthieu Constant's avatar
Matthieu Constant committed
861
  clock_get(&tic);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
862
  fluid_problem_compute_porosity(problem);
Matthieu Constant's avatar
Matthieu Constant committed
863
864
  clock_get(&toc);
  printf("Time compute_porosity : %.6e \n",clock_diff(&tic,&toc));
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
865
866
867
868
869
870
871
872
873
874
}

int *fluid_problem_particle_element_id(FluidProblem *problem)
{
  return problem->particle_element_id;
}

int fluid_problem_n_particles(FluidProblem *problem)
{
  return problem->n_particles;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
875
}