fluid_problem.c 17 KB
Newer Older
1
2
3
4
5
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "fluid_problem.h"
6
#include "mesh_find.h"
7
8
9
10
11

#define n_fields 3

#define newton_max_it 20
#define newton_atol 1e-12
12
#define newton_rtol 1e-5
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

#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};

int fluid_problem_export(const FluidProblem *problem, const char *output_dir, double t, double iter)
{
  if(mesh_write_msh_scalar(problem->mesh, output_dir, "porosity", t, iter, problem->porosity, 1, 0))
    return -1;
  if(mesh_write_msh_vector(problem->mesh, output_dir, "uv", t, iter, problem->solution, n_fields, 0, 1))
    return -1;
  if(mesh_write_msh_scalar(problem->mesh, output_dir, "p", t, iter, problem->solution, n_fields, 2))
    return -1;
  return 0;
}

29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
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
static void particle_drag(double u[2], double mu, double rho, double d, double c, double drag[2], double dt, double mass)
{
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
  double Re_O_u = d*c/mu*rho;
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
  double f = pow(c,-1.8);
  double GoU = f*Cd_u*rho/2*d;
  drag[0] = GoU*u[0]/(1+dt/mass*GoU);
  drag[1] = GoU*u[1]/(1+dt/mass*GoU);
}

static void fluid_problem_add_particle_force(FluidProblem *problem, const double *solution, double dt) {
  double *porosity = problem->porosity;
  Mesh *mesh = problem->mesh;
  for (int i = 0; i < problem->n_particles; ++i) {
    double *xi = problem->particle_uvw + 2*i;
    double phi[3] = {1-xi[0]-xi[1],xi[0],xi[1]};
    int iel = problem->particle_element_id[i];
    int *tri = problem->mesh->triangles+iel*3;
    double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
    double U[] = {solution[tri[0]*3+0],solution[tri[1]*3+0],solution[tri[2]*3+0]};
    double V[] = {solution[tri[0]*3+1],solution[tri[1]*3+1],solution[tri[2]*3+1]};
    double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
    double vf[2] = {
      phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2],
      phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2]
    };
    double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2];
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
    const double dxdxi[2][2] = {
      {x[1][0]-x[0][0],x[2][0]-x[0][0]},
      {x[1][1]-x[0][1],x[2][1]-x[0][1]}};
    const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
    const double dxidx[2][2] = {
      {dxdxi[1][1]/detj,-dxdxi[0][1]/detj},
      {-dxdxi[1][0]/detj,dxdxi[0][0]/detj}
    };
    const double dphi[][2] = {
      {-dxidx[0][0]-dxidx[1][0],-dxidx[0][1]-dxidx[1][1]},
      {dxidx[0][0],dxidx[0][1]},
      {dxidx[1][0],dxidx[1][1]}
    };
    double gradp[2] = {
      dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2],
      dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]
    };
    double du[2] = {
      problem->particle_velocity[i*2+0]-vf[0]/c,
      problem->particle_velocity[i*2+1]-vf[1]/c
    };
    double vol = problem->particle_volume[i];
    double mass = problem->particle_mass[i];
    double d = 2*sqrt(vol/M_PI);
    double drag[2];
    particle_drag(du,problem->mu,problem->rho,d,c,drag, dt, mass);
    double ddragdu[2],ddragdv[2];
    double eps = 1e-8;
    double duu[2] = {du[0]+eps, du[1]};
    particle_drag(duu,problem->mu,problem->rho,d,c,ddragdu, dt, mass);
    double duv[2] = {du[0], du[1]+eps};
    particle_drag(duv,problem->mu,problem->rho,d,c,ddragdv, dt, mass);
    ddragdu[0] = -(ddragdu[0]-drag[0])/(eps*c);
    ddragdu[1] = -(ddragdu[1]-drag[1])/(eps*c);
    ddragdv[0] = -(ddragdv[0]-drag[0])/(eps*c);
    ddragdv[1] = -(ddragdv[1]-drag[1])/(eps*c);
    problem->particle_force[i*2+0] = -drag[0]-vol*gradp[0];
    problem->particle_force[i*2+1] = -drag[1]-9.81*(mass-problem->rho*vol)-vol*gradp[1];
    double local_vector[3*n_fields];
    double local_matrix[9*n_fields*n_fields];
    for (int iphi = 0; iphi < 3; ++iphi) {
      local_vector[iphi+0] = -drag[0]*phi[iphi];
      local_vector[iphi+3] = -drag[1]*phi[iphi];
      local_vector[iphi+6] = 0;
      for (int jphi = 0; jphi < 3; ++jphi) {
        local_matrix[(iphi+0)*9+jphi+0] = -phi[iphi]*ddragdu[0]*phi[jphi];
        local_matrix[(iphi+0)*9+jphi+3] = -phi[iphi]*ddragdv[0]*phi[jphi];
        local_matrix[(iphi+0)*9+jphi+6] = 0;//phi[iphi]*vol*dphi[jphi][0];
        local_matrix[(iphi+3)*9+jphi+0] = -phi[iphi]*ddragdu[1]*phi[jphi];
        local_matrix[(iphi+3)*9+jphi+3] = -phi[iphi]*ddragdv[1]*phi[jphi];
        local_matrix[(iphi+3)*9+jphi+6] = 0;//phi[iphi]*vol*dphi[jphi][1];
        local_matrix[(iphi+6)*9+jphi+0] = 0;
        local_matrix[(iphi+6)*9+jphi+3] = 0;
        local_matrix[(iphi+6)*9+jphi+6] = 0;
      }
    }
    linear_system_add_to_matrix(problem->linear_system,iel,iel,local_matrix);
    linear_system_add_to_rhs(problem->linear_system,iel,local_vector);
  }
}


125
126
127
128
129
130
static void fluid_problem_assemble_system(FluidProblem *problem, const double *solution_old, double dt)
{
  LinearSystem *lsys = problem->linear_system;
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  const double *porosity = problem->porosity;
131
  const double *old_porosity = problem->old_porosity;
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
  double *local_vector = malloc(sizeof(double)*n_fields*3);
  double *local_matrix = malloc(sizeof(double)*n_fields*n_fields*9);
  const double mu = problem->mu;
  const double rho = problem->rho;
  const double epsilon = problem->epsilon;
  for (int iel=0; iel < mesh->n_triangles; ++iel) {
    const int *tri = &mesh->triangles[iel*3];
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
    const double dxdxi[2][2] = {
      {x[1][0]-x[0][0],x[2][0]-x[0][0]},
      {x[1][1]-x[0][1],x[2][1]-x[0][1]}};
    const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
    const double dxidx[2][2] = {
      {dxdxi[1][1]/detj,-dxdxi[0][1]/detj},
      {-dxdxi[1][0]/detj,dxdxi[0][0]/detj}
    };
    const double dphi[][2] = {
      {-dxidx[0][0]-dxidx[1][0],-dxidx[0][1]-dxidx[1][1]},
      {dxidx[0][0],dxidx[0][1]},
      {dxidx[1][0],dxidx[1][1]}
    };
156
157
158
159
    for (int i=0; i < 3*n_fields; ++i)
      local_vector[i] = 0;
    for (int i=0; i < 9*n_fields*n_fields; ++i)
      local_matrix[i] = 0;
160
    double C[] = {porosity[tri[0]],porosity[tri[1]],porosity[tri[2]]};
161
162
    double OLDC[] = {old_porosity[tri[0]],old_porosity[tri[1]],old_porosity[tri[2]]};
    double DCDT[] = {(C[0]-OLDC[0])/dt,(C[1]-OLDC[1])/dt,(C[2]-OLDC[2])/dt};
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
    double U[] = {solution[tri[0]*3+0],solution[tri[1]*3+0],solution[tri[2]*3+0]};
    double V[] = {solution[tri[0]*3+1],solution[tri[1]*3+1],solution[tri[2]*3+1]};
    double U_old[] = {solution_old[tri[0]*3+0],solution_old[tri[1]*3+0],solution_old[tri[2]*3+0]};
    double V_old[] = {solution_old[tri[0]*3+1],solution_old[tri[1]*3+1],solution_old[tri[2]*3+1]};
    double P[] = {solution[tri[0]*3+2],solution[tri[1]*3+2],solution[tri[2]*3+2]};
    double dc[2] = {
      dphi[0][0]*C[0]+dphi[1][0]*C[1]+dphi[2][0]*C[2],
      dphi[0][1]*C[0]+dphi[1][1]*C[1]+dphi[2][1]*C[2]
    };
    double du[2] = {
      dphi[0][0]*U[0]+dphi[1][0]*U[1]+dphi[2][0]*U[2],
      dphi[0][1]*U[0]+dphi[1][1]*U[1]+dphi[2][1]*U[2]
    };
    double dv[2] = {
      dphi[0][0]*V[0]+dphi[1][0]*V[1]+dphi[2][0]*V[2],
      dphi[0][1]*V[0]+dphi[1][1]*V[1]+dphi[2][1]*V[2]
    };
    double dp[2] = {
      dphi[0][0]*P[0]+dphi[1][0]*P[1]+dphi[2][0]*P[2],
      dphi[0][1]*P[0]+dphi[1][1]*P[1]+dphi[2][1]*P[2]
    };
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      const double xi = QP[iqp][0], eta = QP[iqp][1];
      const double phi[] = {1-xi-eta, xi, eta};
      double u = phi[0]*U[0]+phi[1]*U[1]+phi[2]*U[2];
      double v = phi[0]*V[0]+phi[1]*V[1]+phi[2]*V[2];
      double u_old = phi[0]*U_old[0]+phi[1]*U_old[1]+phi[2]*U_old[2];
      double v_old = phi[0]*V_old[0]+phi[1]*V_old[1]+phi[2]*V_old[2];
      double p = phi[0]*P[0]+phi[1]*P[1]+phi[2]*P[2];
      double c = phi[0]*C[0]+phi[1]*C[1]+phi[2]*C[2];
193
      double dcdt = phi[0]*DCDT[0]+phi[1]*DCDT[1]+phi[2]*DCDT[2];
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
      const double jw = QW[iqp]*detj;
      for (int iphi = 0; iphi < 3; ++iphi) {
        double phii = phi[iphi];
        double dphii[2] = {dphi[iphi][0],dphi[iphi][1]};
        double tau[2][2] = {{du[0]-u*dc[0]/c,du[1]-u*dc[1]/c},
                            {dv[0]-v*dc[0]/c,dv[1]-v*dc[1]/c}};
        local_vector[iphi+0] += jw*(
           mu*(dphii[0]*tau[0][0]+dphii[1]*(tau[0][1]+tau[1][0])*0.5)
          +rho*phii*((u-u_old)/dt+(u*tau[0][0]+v*tau[0][1])/c)
          -p*dphii[0]
        );
        local_vector[iphi+3] += jw*(
           mu*(dphii[0]*(tau[0][1]+tau[1][0])*0.5+dphii[1]*tau[1][1])
          +rho*phii*((v-v_old)/dt+(u*tau[1][0]+v*tau[1][1])/c)
          -p*dphii[1]
        );
        local_vector[iphi+6] += jw*(
           epsilon*(dphii[0]*dp[0]+dphii[1]*dp[1])
212
           +phii*((du[0]+dv[1])+dcdt)
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
        );
        for (int jphi = 0; jphi < 3; ++jphi) {
          double phij = phi[jphi];
          double dphij[2] = {dphi[jphi][0],dphi[jphi][1]};
          double dtau[2] = {dphij[0]-phij*dc[0]/c,dphij[1]-phij*dc[1]/c};
          local_matrix[(iphi+0)*9+jphi+0] += jw*mu*(dphii[0]*dtau[0]+dphii[1]*dtau[1]*0.5)    // U U
                                            +jw*rho*phii*(phij/dt+(phij*tau[0][0]+u*dtau[0]+v*dtau[1])/c);
          local_matrix[(iphi+0)*9+jphi+3] += jw*mu*dphii[1]*dtau[0]*0.5                       // U V
                                            +jw*rho*phii*phij*tau[0][1]/c;
          local_matrix[(iphi+0)*9+jphi+6] +=-jw*dphii[0]*phij;                                // U P
          local_matrix[(iphi+3)*9+jphi+0] += jw*mu*dphii[0]*dtau[1]*0.5                       // V U
                                            +jw*rho*phii*phij*tau[1][0]/c;
          local_matrix[(iphi+3)*9+jphi+3] += jw*mu*(dphii[0]*dtau[0]*0.5+dphii[1]*dtau[1])    // V V
                                            +jw*rho*phii*(phij/dt+(phij*tau[1][1]+u*dtau[0]+v*dtau[1])/c);
          local_matrix[(iphi+3)*9+jphi+6] +=-jw*dphii[1]*phij;                                // V P
          local_matrix[(iphi+6)*9+jphi+0] += jw*phii*dphij[0];                                // P U
          local_matrix[(iphi+6)*9+jphi+3] += jw*phii*dphij[1];                                // P V
          local_matrix[(iphi+6)*9+jphi+6] += jw*epsilon*(dphii[0]*dphij[0]+dphii[1]*dphij[1]);// P P
        }
      }
    }
    linear_system_add_to_matrix(lsys,iel,iel,local_matrix);
    linear_system_add_to_rhs(lsys,iel,local_vector);
  }
  free(local_vector);
  free(local_matrix);
239
  fluid_problem_add_particle_force(problem, solution, dt);
240
241
}

242
static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nBnd, StrongBoundary *bnds)
243
{
244
245
  for (int ibnd = 0; ibnd < nBnd; ++ibnd) {
    StrongBoundary *bnd = bnds + ibnd;
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
    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)
      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = bnd->v;
  }
}

void fluid_problem_implicit_euler(FluidProblem *problem, double dt)
{
  const Mesh *mesh = problem->mesh;
262
  apply_boundary_conditions(mesh, problem->solution, problem->n_strong_boundaries, problem->strong_boundaries);
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
  double *solution_old = malloc(mesh->n_nodes*n_fields*sizeof(double));
  for (int i=0; i<mesh->n_nodes*n_fields; ++i)
    solution_old[i] = problem->solution[i];
  double *dx = malloc(sizeof(double)*mesh->n_nodes*n_fields);
  double norm0 = 0;
  for (int i=0; i<newton_max_it; ++i) {
    linear_system_zero_matrix_and_rhs(problem->linear_system);
    fluid_problem_assemble_system(problem, solution_old, dt);
    double norm = linear_system_get_rhs_norm(problem->linear_system);
    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;
    linear_system_solve(problem->linear_system, dx);
    for (int j=0; j<mesh->n_nodes*n_fields; ++j) {
      problem->solution[j] -= dx[j];
    }
  }
  free(dx);
  free(solution_old);
}

288
FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rho, double epsilon, int n_strong_boundaries, StrongBoundary *strong_boundaries)
289
290
291
292
293
294
295
296
297
{
  FluidProblem *problem = malloc(sizeof(FluidProblem));
  problem->epsilon = epsilon;
  problem->rho = rho;
  problem->mu = mu;
  Mesh *mesh = mesh_load_msh(mesh_file_name);
  if (!mesh)
    return NULL;
  problem->mesh = mesh;
298
299
300
301
302
303
304
  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);
    problem->strong_boundaries[i].v = strong_boundaries[i].v;
    problem->strong_boundaries[i].field = strong_boundaries[i].field;
  }
305
  problem->porosity = malloc(mesh->n_nodes*sizeof(double));
306
  problem->old_porosity = malloc(mesh->n_nodes*sizeof(double));
307
308
  problem->solution = malloc(mesh->n_nodes*n_fields*sizeof(double));
  // begin to remove
309
310
  //static StrongBoundary strong_boundaries[] = {{"Bottom",1,VEXT},{"Bottom",0,0.},{"Top",0,0.},{"PtFix",2,0.},{"Box",0,0.},{NULL, 0, 0.}};
  //static StrongBoundary strong_boundaries[] = {{"Top",2,0},{"Top",1,0.},{"Top",0,0.},{"Box",1,0.},{"Box",0,0.},{NULL, 0, 0.}};
311
312
313
314
315
316
  for (int i = 0; i < problem->mesh->n_nodes; ++i){
    problem->solution[i*3+0] = 0.;
    problem->solution[i*3+1] = 0;
    problem->solution[i*3+2] = 0.;
  }
  // end to remove
317
318
319
320
321
322
323
324
325
  problem->linear_system = linear_system_new(mesh, n_fields, problem->n_strong_boundaries, problem->strong_boundaries); 
  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);
  problem->particle_force = malloc(0);
  problem->particle_volume = malloc(0);
  problem->particle_velocity = malloc(0);
326
327
328
329
330
331
332
  return problem;
}

void fluid_problem_free(FluidProblem *problem)
{
  free(problem->solution);
  free(problem->porosity);
333
  free(problem->old_porosity);
334
  linear_system_free(problem->linear_system);
335
336
337
338
339
340
341
342
343
344
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
  free(problem->particle_force);
  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);
345
346
347
  mesh_free(problem->mesh);
}

348
void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity)
349
{
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
  problem->n_particles = n;
  free(problem->particle_uvw);
  free(problem->particle_element_id);
  free(problem->particle_position);
  free(problem->particle_mass);
  free(problem->particle_force);
  free(problem->particle_volume);
  free(problem->particle_velocity);
  problem->particle_uvw = malloc(sizeof(double)*2*n);
  problem->particle_element_id = malloc(sizeof(int)*n);
  problem->particle_position = malloc(sizeof(double)*2*n);
  problem->particle_mass = malloc(sizeof(double)*n);
  problem->particle_force = malloc(sizeof(double)*2*n);
  problem->particle_volume = malloc(sizeof(double)*n);
  problem->particle_velocity = malloc(sizeof(double)*n*2);
  mesh_particles_to_mesh(problem->mesh, n, position, problem->particle_element_id, problem->particle_uvw);
  for (int i = 0; i < n; ++i) {
    problem->particle_position[i*2+0] = position[i*2+0];
    problem->particle_position[i*2+1] = position[i*2+1];
    problem->particle_mass[i] = mass[i];
    problem->particle_volume[i] = volume[i];
    problem->particle_velocity[i*2+0] = velocity[i*2+0];
    problem->particle_velocity[i*2+1] = velocity[i*2+1];
  }
  Mesh *mesh = problem->mesh;
  double *vertex_volume = malloc(sizeof(double)*mesh->n_nodes);
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->old_porosity[i] = problem->porosity[i];
    problem->porosity[i] = 0;
    vertex_volume[i] = 0;
  }
  for (int iel = 0; iel < mesh->n_triangles; ++iel) {
    const int *tri = &mesh->triangles[iel*3];
    const double x[3][2] = {
      {mesh->x[tri[0]*3+0],mesh->x[tri[0]*3+1]},
      {mesh->x[tri[1]*3+0],mesh->x[tri[1]*3+1]},
      {mesh->x[tri[2]*3+0],mesh->x[tri[2]*3+1]}};
    const double dxdxi[2][2] = {
      {x[1][0]-x[0][0],x[2][0]-x[0][0]},
      {x[1][1]-x[0][1],x[2][1]-x[0][1]}};
    const double detj = dxdxi[0][0]*dxdxi[1][1]-dxdxi[0][1]*dxdxi[1][0];
    vertex_volume[tri[0]] += detj/6;
    vertex_volume[tri[1]] += detj/6;
    vertex_volume[tri[2]] += detj/6;
  }
  for (int i = 0; i < n; ++i) {
    const int *tri = &mesh->triangles[problem->particle_element_id[i]*3];
    double u = problem->particle_uvw[i*2+0], v = problem->particle_uvw[i*2+1];
    problem->porosity[tri[0]] += (1-u-v)*volume[i];
    problem->porosity[tri[1]] += u*volume[i];
    problem->porosity[tri[2]] += v*volume[i];
  }
  for (int i = 0; i < mesh->n_nodes; ++i){
    problem->porosity[i] = (1-problem->porosity[i]/vertex_volume[i]);
  }
  free(vertex_volume);
406
}