fluid_problem.c

/*
 * MigFlow - Copyright (C) <2010-2018>
 * <Universite catholique de Louvain (UCL), Belgium
 *  Universite de Montpellier, France>
 * 	
 * List of the contributors to the development of MigFlow: see AUTHORS file.
 * Description and complete License: see LICENSE file.
 * 	
 * This program (MigFlow) is free software: 
 * you can redistribute it and/or modify it under the terms of the GNU Lesser General 
 * 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
 * GNU Lesser General Public License for more details.
 * 
 * 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, 
 * see <http://www.gnu.org/licenses/>.
 */

#include "tools.h"
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
#include "fluid_problem.h"
#include "mesh_find.h"
#include "mbxml.h"

#define D DIMENSION

#define  deps 1e-8
#if D==2

#define N_SF 3
#define N_N 3
#define N_QP 3
#define N_LQP 2

static double LQP[] = {0.21132486540518713, 0.7886751345948129};
static double LQW[] = {0.5,0.5};
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};
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;
}
#else

#define N_SF 4
#define N_N 4
#define N_QP 5
#define N_LQP 3
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 LQP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
static double LQW[] = {1./6,1./6,1./6};
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;
};

#endif

#define P  D
#define U 0
#define A (D+1)

int fluid_problem_n_fields(const FluidProblem *p) {return D + p->n_fluids;}

#if D==2
static double particle_drag_coeff(double u[2], double mu, double rho, double vol, double c)
{
  double d = 2*sqrt(vol/M_PI);
  double normu = hypot(u[0],u[1]);
  //Reynolds/|u_p-u_f|
  double Re_O_u = d*c/mu*rho;
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
  double f = pow(c,-1.8);
  return f*Cd_u*rho/2*d;
}
#else
static double particle_drag_coeff(double u[3], double mu, double rho, double vol, double c)
{
  double d = 2*cbrt(3./4.*vol/M_PI);
  double normu = sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]);
  //Reynolds/|u_p-u_f|
  double Re_O_u = d*c/mu*rho;
  double Cd_u = 0.63*sqrt(normu)+4.8/sqrt(Re_O_u);
  Cd_u = Cd_u*Cd_u;
  double f = pow(c,-1.8);
  double r = d/2.;
  return f*Cd_u*rho/2*(M_PI*r*r);
}
#endif

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

static void particle_force_f(FluidProblem *problem, double *f0, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int ip, int in_derivative, int reduced_gravity) {

  double p = sol[P];
  double u[D], uold[D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    dp[iD] = dsol[P*D+iD];
  }
  double Ra, a, mu, rho;

  if (problem->n_fluids == 1){
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    a = sol[A];
    a = fmin(1.,fmax(0.,a));
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
    mu = problem->mu[0]*a + problem->mu[1]*(1-a);
  }

  double Du[D], Due[D];
  for (int j = 0; j < D; ++j) {
    Du[j] = problem->particle_velocity[ip*D+j]-u[j]/c;
    Due[j] = problem->particle_velocity[ip*D+j]-uold[j]/c;
  }
  double mass = problem->particle_mass[ip];
  double vol = problem->particle_volume[ip];
  double g = problem->g;
  if (reduced_gravity){
    double rhop = mass/vol;
    g = g*(rhop-problem->rho[0])/rhop;
  }
  double gamma;
  
  if(problem->n_fluids == 2){
    double gamma0 = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
    double gamma1 = particle_drag_coeff(Due,problem->mu[1],problem->rho[1],vol,c);
    gamma = gamma0*a + gamma1*(1-a);
  }
  else{
    gamma = particle_drag_coeff(Due,problem->mu[0],problem->rho[0],vol,c);
  }
  double fcoeff = mass/(mass+dt*gamma);
  for (int j = 0; j < D; ++j){
    if(!in_derivative)
      problem->particle_force[ip*D+j] = fcoeff*(-gamma*Du[j]-vol*dp[j]);
    f0[U+j] = -fcoeff*gamma*(Du[j]-dt/mass*vol*dp[j])-vol*dp[j];
  }
  if(!in_derivative)
    problem->particle_force[ip*D+1] += fcoeff*g*mass;
  f0[U+1] += -fcoeff*gamma*dt*g;
}

static void node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity)
{
  const Mesh *mesh = problem->mesh;
  const double *porosity = problem->porosity;

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *g0 = malloc(sizeof(double)*n_fields);
  double *h0 = malloc(sizeof(double)*n_fields*D);

  for (int i=0;i<n_fields;++i){
    f0[i] = 0;
    g0[i] = 0;
    for (int j = 0; j < D; ++j)
      h0[i*D+j] = 0;
  }
  for (int ip=0; ip < problem->n_particles; ++ip) {
    int iel = problem->particle_element_id[ip];
    if(iel < 0){
      for (int d = 0; d < D; ++d)
        problem->particle_force[ip*D+d] = 0;
      problem->particle_force[ip*D+1] = problem->g*problem->particle_mass[ip];
      continue;
    }
    double *xi = problem->particle_uvw + D*ip;
    double phi[N_SF];
    shape_functions(xi,phi);
    const unsigned int *el = &mesh->elements[iel*N_N];
    double dxdxi[D][D], dxidx[D][D], dphi[N_SF][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];
    const double detj = invDD(dxdxi, dxidx);
    grad_shape_functions(dxidx, dphi);
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
    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;
    double cold = 0;
    double dc[D] = {0};
    for (int i = 0; i < N_SF; ++i) {
      c += problem->porosity[el[i]]*phi[i];
      cold += problem->oldporosity[el[i]]*phi[i];
      for (int j = 0; j < n_fields; ++j) {
        double dof = solution[el[i]*n_fields+j];
        double dofold = solution_old[el[i]*n_fields+j];
        s[j] += phi[i]*dof;
        sold[j] += phi[i]*dofold;
        for (int k = 0; k < D; ++k)
          ds[j*D+k] += dphi[i][k]*dof;
      }
      for (int k=0; k<D; ++k){
        dc[k] += problem->porosity[el[i]]*dphi[i][k];
      }
    }
    particle_force_f(problem,f0,s,ds,sold,c,dc,cold,dt,iel,ip,0,reduced_gravity);
    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];
      }
    }
    for (int jfield = 0; jfield < n_fields; ++jfield) {
      double store = s[jfield];
      s[jfield] += deps;
      particle_force_f(problem,g0,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
      s[jfield] = store;
      for (int jd =0; jd < D; ++jd){
        store = ds[jfield*D+jd];
        ds[jfield*D+jd] += deps;
        particle_force_f(problem,h0+jd*n_fields,s,ds,sold,c,dc,cold,dt,iel,ip,1,reduced_gravity);
        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 = phi[iphi]*phi[jphi]*df00;
            for (int jd = 0; jd < D; ++jd) {
              double df01 = (h0[jd*n_fields+ifield]-f0[ifield])/deps;
              m += phi[iphi]*dphi[jphi][jd]*df01;
            }
            int IDX = (iphi*N2)*n_fields + jphi*n_fields;
            LOCAL_MATRIX2(ifield,jfield) += m;
          }
        }
      }
    }
  }
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
}

static void f_pressure(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
  double p = s[P];
  double u[D], du[D][D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
  }

  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
  double un = 0;
  double rho, mu;
  for (int i = 0; i < D; ++i) {
    un += u[i]*n[i];
    divu += du[i][i];
    for (int j = 0; j < D; ++j) {
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
      utau[i] += u[j]*tau[i][j];
    }
  }
  
  f[P] = un;
  if(problem->n_fluids == 1){
    mu = problem->mu[0];
    rho = problem->rho[0];
  }
  else{
    double a = s[A];
    a = fmin(1.,fmax(0.,a));
    rho = problem->rho[0]*a+problem->rho[1]*(1-a);
    mu = problem->mu[0]*a+problem->mu[1]*(1-a);
    f[A] = un*(un > 0 ? a : v[1]);
  }
  for (int id = 0; id < D; ++id) {
    f[U+id] = v[0]*n[id] + (un> 0 ? rho*u[id]*un/c : 0);
  }
}

static void f_velocity(FluidProblem *problem, const double *n, double *f, const double *s, const double *ds, const double c, const double *dc, const double dt, int eid, const double *v){
  double p = s[P];
  double u[D], du[D][D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
  }

  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
  double rho, mu, rho_ext;
  double unint = 0;
  
  double un = 0;
  for (int i = 0; i < D; ++i) {
    un += c*v[i]*n[i];
    unint += u[i]*n[i];
    divu += du[i][i];
    for (int j = 0; j < D; ++j) {
      tau[i][j] = du[i][j]-u[i]*dc[j]/c;
      utau[i] += u[j]*tau[i][j];
    }
  }
  f[P] = un;
  if(problem->n_fluids == 1){
    mu = problem->mu[0];
    rho = problem->rho[0];
    rho_ext = rho;
  }
  else{
    double a = s[A];
    a = fmin(1.,fmax(0.,a));
    rho = problem->rho[0]*a+problem->rho[1]*(1-a);
    rho_ext = problem->rho[0]*v[D]+problem->rho[1]*(1-v[D]);
    mu = problem->mu[0]*a+problem->mu[1]*(1-a);
    f[A] = un*(un > 0 ? a : v[D]);
    //printf("%g\n",v[D]);
  }

  for (int id = 0; id < D; ++id) {
    f[U+id] = p*n[id] + (un> 0 ? rho*u[id]*unint : rho_ext*un*un*n[id]);
  }
}

static void f_wall(FluidProblem *problem,const double *n, double *f,const double *s,const double *ds, const double c, const double *dc, const double dt, int eid, const double *v)
{
  double p = s[P];
  double u[D], du[D][D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
  }
  double un = 0;
  double a;
  for (int i = 0; i < D; ++i) {
    un += u[i]*n[i];
  }
  double rho;
  if (problem->n_fluids == 1){
    rho = problem->rho[0];
  }
  else{
    a = s[A];
    a = fmin(1.,fmax(0.,a));
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
  }

  f[P] = 0;

  if(problem->n_fluids == 2){
    double a = s[A];
    f[A] = 0;
  }

  for (int id = 0; id < D; ++id) {
    f[U+id] = (v==NULL ? p : v[0])*n[id]+ (un> 0 ? rho*u[id]*un/c : 0);
  }
}

static double pow2(double a) {return a*a;}
static void compute_stab_parameters(FluidProblem *problem, double dt) {
  const Mesh *mesh = problem->mesh;

  problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
  problem->tauc = realloc(problem->tauc,sizeof(double)*mesh->n_elements);
  problem->taua = realloc(problem->taua,sizeof(double)*mesh->n_elements);
  problem->extraa = realloc(problem->extraa,sizeof(double)*mesh->n_elements);

  double maxtaup = 0;
  double mintaup = DBL_MAX;
  const int n_fields = fluid_problem_n_fields(problem);
  double maxtaa = 0;
  double mintaa = DBL_MAX;
  for (int iel = 0; iel < mesh->n_elements; ++iel) {
    double dxdxi[DIMENSION][DIMENSION];
    double normu = {0.};
    const uint32_t *el = &mesh->elements[iel*N_N];
    double a = 0;
    double amax = 0;
    double amin = 0;
    double c = 0;
    for (int i=0; i< N_N; ++i) {
      double normup = 0;
      c += problem->porosity[el[i]]/N_N;
      if(problem->n_fluids == 2){
        a += problem->solution[el[i]*n_fields+A];
        amax = fmax(amax,problem->solution[el[i]*n_fields+A]);
        amin = fmin(amin,problem->solution[el[i]*n_fields+A]);
      }

      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
      normu = fmax(normu,sqrt(normup));
    }
    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 h = 2*sqrt(detDD(dxdxi)/(2*M_PI));
#else
    const double h = 2*cbrt(detDD(dxdxi)/(8*M_PI));
#endif
    double rho, mu;
    if(problem->n_fluids == 1){
      rho = problem->rho[0];
      mu = problem->mu[0];
    }else{
      a /= N_N;
      a = fmin(1.,fmax(0.,a));
      rho = problem->rho[0]*a + problem->rho[1]*(1-a);
      mu = problem->mu[0]*a + problem->mu[1]*(1-a);
    }
    double nu = mu/rho;
    problem->taup[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h)+pow2(4*nu/pow2(h)));
    problem->tauc[iel] = h*normu*fmin(h*normu/(6*nu),0.5);
    problem->taua[iel] = 1/sqrt(pow2(2/dt)+pow2(2*normu/h));
    problem->extraa[iel] = fmax(0.,fmax(-amin+0.01,amax-1+0.01))*problem->coeffStab;
  }
}

static void fluid_problem_f(FluidProblem *problem, double *f0, double *f1, double *sol, double *dsol, const double *solold, const double c, const double *dc, const double cold, double dt, int iel, int reduced_gravity, double stab_param) {
  double p = sol[P];
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
  double taua = problem->taua[iel];
  double extraa = problem->extraa[iel];
  double u[D], uold[D], du[D][D], dp[D];
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    dp[iD] = dsol[P*D+iD];
    for (int jD = 0; jD < D; ++jD) du[iD][jD] = dsol[(U+iD)*D+jD];
  }
  double Ra, a, aold, mu, rho, rhoold, drhodt=0;

  if (problem->n_fluids == 1){
    rho = problem->rho[0];
    mu = problem->mu[0];
  }
  else{
    a = sol[A];
    a = fmin(1.,fmax(0.,a));
    aold = solold[A];
    aold = fmin(1.,fmax(0.,aold));
    Ra = (a*c-aold*cold)/dt;
    rho = problem->rho[0]*a + problem->rho[1]*(1-a);
    rhoold = problem->rho[0]*aold + problem->rho[1]*(1-aold);
    mu = problem->mu[0]*a + problem->mu[1]*(1-a);
  }
  double divu = 0;
  double tau[D][D];
  double utau[D]={0.};
  double da[D];
  double uugradrho[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]*dc[j]/c;
      utau[i] += u[j]*tau[i][j];
    }
    if(problem->n_fluids == 2){
      da[i] = dsol[A*D+i];
      Ra += u[i]*da[i];
    }
  }
  if(problem->n_fluids == 2){
    Ra += a*divu;
    drhodt = (rho-rhoold)/dt;
    for (int i = 0; i < D; ++i) {
      for (int j=0;j<D;++j){
        uugradrho[i] += u[i]*u[j]/c *(problem->rho[0]*da[j] - problem->rho[1]*da[j]);
      }
    }
  }

  f0[P] = (c-cold)/dt;

  if(problem->n_fluids == 2){
    f0[A] = (a*c-aold*cold)/dt;
  }
  
  for (int i = 0; i < D; ++i) {
    double dudt = (u[i]-uold[i])/dt;
    double R = dudt + u[i]*drhodt/rho + dp[i]/rho + u[i]*divu/c + utau[i]/c + uugradrho[i]/rho + (i==1 ? -(reduced_gravity ? 0 : 1)*c*problem->g : 0.);
    f0[U+i] = rho*dudt + u[i]*drhodt + (i==1 ? -(reduced_gravity ? 0 : 1)*c*rho*problem->g : 0.);
    for (int j = 0; j < D; ++j) {
      f1[(U+i)*D+j] = -rho*u[i]*u[j]/c + mu*(tau[i][j]+tau[j][i]) + (i==j ? -p : 0) + (stab_param==0 ? 1 : 0)*R*u[j]*rho*taup;
    }
    f1[(U+i)*D+i] += (stab_param==0 ? 1 : 0)*divu*tauc*rho;
    f1[P*D+i] = -u[i] + (stab_param==0 ? taup*R : stab_param*dp[i]);
    if(problem->n_fluids == 2){
      f1[P*D+i] +=  taua*(divu+(c-cold)/dt)*u[i];
      f1[A*D+i] = -u[i]*a + taua*Ra*u[i] + taup*R*a + da[i]*(extraa);
    }
  }
  if(problem->n_fluids == 2){
    f1[A*D+1] += -fmax(a,0)*fmax((1-a),0)*1e-4;
  }
}

static void compute_weak_boundary_conditions(FluidProblem *problem, double dt, double *all_local_vector, double *all_local_matrix)
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  const int n_fields = fluid_problem_n_fields(problem);
  const size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *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);
  int *bnd_node_local_id = malloc(sizeof(int)*mesh->n_nodes);
  f_cb *f;
#if DIMENSION == 2
  int elbnd[3][2] = {{0,1},{1,2},{2,0}};
#else
  int elbnd[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,3,2}};
#endif
  for (int iphys = 0; iphys < mesh->n_phys; ++iphys) {
    for (int i =0; i < mesh->n_nodes; ++i)
      bnd_node_local_id[i] = -1;
    if(mesh->phys_dimension[iphys] != D-1) continue;
    double *x = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*D);
    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
      int j = mesh->phys_nodes[iphys][i];
      bnd_node_local_id[j] = i;
      for (int k = 0; k < D; ++k)
        x[i*D+k] = mesh->x[j*3+k];
    }
    
    f = NULL;
    int n_value;
    double *v = (double*)malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
    for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
      WeakBoundary *bnd = problem->weak_boundaries + ibnd;
      if (strcmp(mesh->phys_name[iphys],bnd->tag) == 0){
        f = bnd->cb;
        n_value = bnd->n_value;
        if (n_value > 0){
          v = realloc(v,mesh->phys_n_nodes[iphys]*sizeof(double)*n_value);
          bnd->apply(mesh->phys_n_nodes[iphys], x, v);
        }
        break;
      }
    }
    if (f == NULL) {
      printf("no weak boundary condition defined for tag \"%s\".\n", mesh->phys_name[iphys]);
      exit(1);
    }

    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 iebnd = bnd[1];
      int nodes[D];
      double *x[D];
      const unsigned int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[elbnd[iebnd][i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      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);

#if DIMENSION == 2 
      const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
      const double detbnd = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
      const double n[2] = {dx[1]/detbnd,-dx[0]/detbnd};
#else
      const double dx[2][3] = {{x[1][0]-x[0][0], x[1][1]-x[0][1], x[1][2]-x[0][2]},{x[2][0]-x[0][0], x[2][1]-x[0][1], x[2][2]-x[0][2]}};
      const double nn[3] = {
        dx[0][1]*dx[1][2] - dx[0][2]*dx[1][1],
        dx[0][2]*dx[1][0] - dx[0][0]*dx[1][2],
        dx[0][0]*dx[1][1] - dx[0][1]*dx[1][0]};
      const double detbnd = sqrt(nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2]);
      const double n[3] = {-nn[0]/detbnd,-nn[1]/detbnd,-nn[2]/detbnd};
#endif
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double xi[D];
#if DIMENSION == 2 
        double a = LQP[iqp];
        switch(iebnd) {
          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:;
        }
#else
        double a = LQP[iqp][0];
        double b = LQP[iqp][1];
        switch(iebnd) {
          case 0 : xi[0] = a;     xi[1] = b; xi[2] = 0; break;
          case 1 : xi[0] = 0;     xi[1] = a; xi[2] = b; break;
          case 2 : xi[0] = b;     xi[1] = 0; xi[2] = a; break;
          case 3 : xi[0] = 1-a-b; xi[1] = b; xi[2] = a; break;
          default:;
        }
#endif
        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;
        double *vbnd = NULL;
        if (n_value>0){
          vbnd = malloc(sizeof(double)*n_value);
        }
        for (int i = 0; i < n_value; ++i){
          vbnd[i] = 0;
        }
        double dc[D] = {0};
        for (int i = 0; i < N_SF; ++i) {
          c += problem->porosity[el[i]]*phi[i];
          int node_local_id = bnd_node_local_id[el[i]];
          if (node_local_id != -1)
            for (int j = 0; j < n_value; j++){
              vbnd[j] += v[node_local_id*n_value+j]*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;
          }
          for (int k=0; k<D; ++k){
            dc[k] += problem->porosity[el[i]]*dphi[i][k];
          }
        }
        const double jw = LQW[iqp]*detbnd;
        f(problem,n,f0,s,ds,c,dc,dt,iel,vbnd);
        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,dc,dt,iel,vbnd);
          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,dc,dt,iel,vbnd);
            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;
              }
            }
          }
        }
        free(vbnd);
      }
    }
    free(x);
    free(v);
  }
  free(s);
  free(ds);
  free(sold);
  free(f0);
  free(g0);
  free(h0);
}


static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix, int reduced_gravity, double stab_param)
{
  const Mesh *mesh = problem->mesh;
  const double *porosity = problem->porosity;

  const double *solution = problem->solution;
  int n_fields = fluid_problem_n_fields(problem);
  size_t local_size = N_SF*n_fields;
  double *s = malloc(sizeof(double)*n_fields);
  double *sold = malloc(sizeof(double)*n_fields);
  double *ds = malloc(sizeof(double)*n_fields*D);
  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);
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    const unsigned int *el = &mesh->elements[iel*N_N];
    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];
    const double detj = invDD(dxdxi, dxidx);
    grad_shape_functions(dxidx, dphi);
    double *local_vector = &all_local_vector[local_size*iel];
    double *local_matrix = &all_local_matrix[local_size*local_size*iel];
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],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;
      double cold = 0;
      double dc[D] = {0};
      for (int i = 0; i < N_SF; ++i) {
        c += problem->porosity[el[i]]*phi[i];
        cold += problem->oldporosity[el[i]]*phi[i];
        for (int j = 0; j < n_fields; ++j) {
          double dof = solution[el[i]*n_fields+j];
          double dofold = solution_old[el[i]*n_fields+j];
          s[j] += phi[i]*dof;
          sold[j] += phi[i]*dofold;
          for (int k = 0; k < D; ++k)
            ds[j*D+k] += dphi[i][k]*dof;
        }
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
      }
      const double jw = QW[iqp]*detj;
      fluid_problem_f(problem,f0,f1,s,ds,sold,c,dc,cold,dt,iel, reduced_gravity, stab_param);
      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;
          }
        }
      }
      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,dc,cold,dt,iel,reduced_gravity,stab_param);
        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,dc,cold,dt,iel,reduced_gravity,stab_param);
          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 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;
                }
              }
              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;
            }
          }