fluid_problem.c

/*
 * MigFlow - Copyright (C) <2010-2020>
 * <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 "fluid_problem_fem.h"

#define P  D
#define U 0

#define LOCAL_MATRIX(i,j,a,b) local_matrix[((i)*n_fields*N_SF+N_SF*(a)+(j))*n_fields+(b)]

static MeshTree *get_mesh_tree(FluidProblem *problem)
{
  if (problem->mesh_tree == NULL)
    problem->mesh_tree = mesh_tree_create(problem->mesh);
  return problem->mesh_tree;
}

static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip);
void fluid_problem_interpolate_rho_and_nu(const FluidProblem *problem, double a, double *rho, double *mu) {
  if (problem->n_fluids == 1){
    *rho = problem->rho[0];
    *mu = problem->mu[0];
  }
  else{
    //*rho = 1/(1/problem->rho[0]*a + 1/problem->rho[1]*(1-a));
    *rho = problem->rho[0]*a + problem->rho[1]*(1-a);
    //*mu = problem->mu[0]*a + problem->mu[1]*(1-a);
    *mu = (problem->mu[0]/problem->rho[0]*a + problem->mu[1]/problem->rho[1]*(1-a))*(*rho);
  }
}

inline static double mesh_dxidx(const Mesh *mesh, int iel, double dxidx[D][D]) {
  double dxdxi[D][D];
  const int *el = &mesh->elements[iel*N_N];
  for (int i = 0; i < D; ++i)
    for (int j = 0; j < D; ++j)
      dxdxi[i][j] = mesh->x[el[j+1]*3+i] - mesh->x[el[0]*3+i];
  return invDD(dxdxi, dxidx);
}

inline static double element_volume_from_detj(double detj) {
#if D == 2
    return detj/2;
#else
    return detj/6;
#endif
}

inline static double node_volume_from_detj(double detj) {
#if DIMENSION == 2
    return detj/6;
#else
    return detj/24;
#endif
}

int fluid_problem_n_fields(const FluidProblem *p) {return D+1;}


void fluid_problem_node_force_volume(FluidProblem *problem, const double *solution_old, double dt, double *all_local_vector, double *all_local_matrix)
{
  const Mesh *mesh = problem->mesh;
  const double *porosity = problem->porosity;
  const double *solution = problem->solution;
  int drag_in_stab = problem->drag_in_stab;
  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);
  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;
      continue;
    }
    double *xi = problem->particle_uvw + D*ip;
    double phi[N_SF];
    shape_functions(xi,phi);
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D], dphi[N_SF][D];
    double detj = mesh_dxidx(mesh, iel, dxidx);
    grad_shape_functions(dxidx, dphi);
    double *local_vector = all_local_vector ? &all_local_vector[local_size*iel] : NULL;
    double *local_matrix = all_local_matrix ? &all_local_matrix[local_size*local_size*iel] : NULL;
    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 a = 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];
      if (problem->n_fluids == 2) {
        a += problem->concentration[iel*N_N+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];
      }
    }
    double f[D],dfdu,dfddp;
    particle_force_f(problem,f,&dfdu,&dfddp,s,ds,sold,c,dc,a,dt,iel,ip);
    if (!local_vector)
      continue;
    double rho,nu;
    fluid_problem_interpolate_rho_and_nu(problem,a, &rho,&nu);
    double pspg = drag_in_stab?problem->taup[iel]/rho:0;
    for (int iD = 0; iD < D; ++iD) {
      for (int iphi = 0; iphi < N_SF; ++iphi){
        local_vector[iphi+N_SF*(U+iD)] += phi[iphi]*f[iD];
        local_vector[iphi+N_SF*P] += pspg*dphi[iphi][iD]*f[iD];
        for (int jD = 0; jD < D; ++jD) {
          double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
          local_vector[iphi+N_SF*(U+iD)] += supg*dphi[iphi][jD]*f[iD];
        }
      }
    }
    if (!local_matrix)
      continue;
    for (int iD = 0; iD < D; ++iD){
      for (int iphi = 0; iphi < N_SF; ++iphi){
        for (int jphi = 0; jphi < N_SF; ++jphi){
          LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += phi[iphi]*phi[jphi]*dfdu;
          LOCAL_MATRIX(iphi,jphi,U+iD,P) += phi[iphi]*dphi[jphi][iD]*dfddp;
          LOCAL_MATRIX(iphi,jphi,P,U+iD) += pspg*dphi[iphi][iD]*phi[jphi]*dfdu;
          LOCAL_MATRIX(iphi,jphi,P,P) += pspg*dphi[iphi][iD]*dphi[jphi][iD]*dfddp;
          for (int jD = 0; jD < D; ++jD) {
            double supg = drag_in_stab?sold[U+jD]*problem->taup[iel]:0;
            LOCAL_MATRIX(iphi,jphi,U+iD,U+iD) += supg*dphi[iphi][jD]*phi[jphi]*dfdu;
            LOCAL_MATRIX(iphi,jphi,U+iD,P) += supg*dphi[iphi][jD]*dphi[jphi][iD]*dfddp;
          }
        }
      }
    }
  }
  free(s);
  free(ds);
  free(sold);
}

static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n, double *f0,const double *s,const double *ds, const double *sold, const double *dsold, const double *mesh_velocity, const double c, const double *dc, const double rho, double mu, const double dt, int eid, const double *data, double *f00, double *f01)
{
  const int vid = wbnd->vid;
  const int pid = wbnd->pid;
  const int cid = wbnd->cid;
  const int aid = wbnd->aid;
  double p = s[P];
  const int n_fields = fluid_problem_n_fields(problem);
  double u[D], dp[D], uext[D];
  double s_c = 0;
  double unold = 0;
  double cext = cid<0?c:data[cid];
  double unext = 0;
  double unmesh = 0;
  double dcn = 0;
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = s[U+iD];
    dp[iD] = ds[P*D+iD];
    unold += sold[U+iD]*n[iD];
    unmesh += mesh_velocity[iD]*n[iD];
    uext[iD] = (vid<0?s[U+iD]:data[vid+iD]*c);
    unext += uext[iD]*n[iD];
    dcn += dc[iD]*n[iD];
    s_c += (u[iD]-uext[iD])*n[iD];
  }
  double h = problem->element_size[eid];
  if (wbnd->type == BND_WALL && pid >= 0) {
      f0[P] = -(s[P]-data[pid])/h*problem->taup[eid];
      f00[P*n_fields+P] += -1/h*problem->taup[eid];
  }
  if(wbnd->type == BND_OPEN) {
    f0[P] = unext;
    for (int i = 0; i < D; ++i) {
      f00[P*n_fields+U+i] += vid<0?n[i]:0;
    }
  }
  for (int id = 0; id < D; ++id) {
    f0[U+id] = c*(pid<0?0:data[pid]-p)*n[id];
    f00[(U+id)*n_fields+P] += c*(pid<0?0:-n[id]);
  }
  double c_du_o_c[D][D];
  double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
  for (int iD = 0; iD < D; ++iD) {
    for (int jD = 0; jD < D; ++jD)
      c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
  }
  if(wbnd->compute_viscous_term == 1){
    for (int id = 0; id < D; ++id) {
      f0[U+id] += sigma*(u[id]-uext[id]+s_c*n[id]);
      f00[(U+id)*n_fields+U+id] += (vid<0?0:sigma);
      for (int jd = 0; jd <D; ++jd) {
        f0[U+id] -= mu*(c_du_o_c[id][jd]+c_du_o_c[jd][id])*n[jd];
        f00[(U+id)*n_fields+U+id] += (vid<0?0:n[id]*n[jd]*sigma) + mu/c*dc[jd]*n[jd];
        f00[(U+id)*n_fields+U+jd] += mu/c*dc[id]*n[jd];
        f01[(U+id)*n_fields*D+(U+id)*D+jd] -= mu*n[jd];
        f01[(U+id)*n_fields*D+(U+jd)*D+id] -= mu*n[jd];
      }
    }
  }
  if (problem->advection) {
    double unbnd = unold;
    if (wbnd->type == BND_WALL) unbnd = unold/2;
    else if(vid>0) unbnd = (unold+unext)/2;
    if (unbnd - unmesh < 0) {
      for (int id = 0; id < D; ++id) {
        f0[U+id] += ((unbnd - unmesh)*(vid<0?0:uext[id])-(unold - unmesh)*u[id])*rho/c;
        f00[(U+id)*n_fields+U+id] -= (unold-unmesh)*rho/c;
      }
    }
  }
}

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

  problem->element_size = realloc(problem->element_size,sizeof(double)*mesh->n_elements);
  problem->taup = realloc(problem->taup,sizeof(double)*mesh->n_elements);
  problem->tauc = realloc(problem->tauc,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 normu = {0.};
    const int *el = &mesh->elements[iel*N_N];
    double a = 0;
    double amax = 0.5;
    double amin = 0.5;
    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->concentration[iel*N_N+i];
        amax = fmax(amax,problem->concentration[iel*N_N+i]);
        amin = fmin(amin,problem->concentration[iel*N_N+i]);
      }
      for(int j = 0; j < DIMENSION; ++j) normup += pow2(problem->solution[el[i]*n_fields+j]);
      normu = fmax(normu,sqrt(normup));
    }
    double dxidx[DIMENSION][DIMENSION];
    double detj = mesh_dxidx(mesh,iel,dxidx);
#if DIMENSION == 2
    const double h = 2*sqrt(detj/(2*M_PI));
#else
    const double h = 2*cbrt(detj/(8*M_PI));
#endif
    problem->element_size[iel] = h;
    double rho, mu;
    a /= N_N;
    a = fmin(1.,fmax(0.,a));
    fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
    double nu = mu/rho;
    problem->taup[iel] = 1/sqrt(
        (problem->temporal?pow2(2/dt):0.)
        +(problem->advection ?pow2(2*normu/h):0.)
        +pow2(4*nu/pow2(h)));
    problem->tauc[iel] = problem->advection?h*normu*fmin(h*normu/(6*nu),0.5):0.;
  }
}

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

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

static void particle_force_f(FluidProblem *problem, double *f, double *dfdu, double *dfddp, double *sol, double *dsol, const double *solold, const double c, const double *dc, double a, double dt, int iel, int ip) {
  const double *contact = &problem->particle_contact[ip*D];
  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 mu, rho;
  fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
  double Due[D];
  const double *up = &problem->particle_velocity[ip*D];
  for (int j = 0; j < D; ++j) {
    Due[j] = problem->particle_velocity[ip*D+j]-uold[j]/c;
  }
  double mass = problem->particle_mass[ip];
  double vol = problem->particle_volume[ip];
  double rhop = mass/vol;
  double massreduced = mass;
  if(problem->reduced_gravity) {
    massreduced = (rhop-problem->rho[0])*problem->particle_volume[ip];
  }
  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,mu,rho,vol,c);
  }
  gamma = gamma*problem->drag_coeff_factor;
  double gammastar = gamma/(1+dt/mass*gamma);
  for (int i = 0; i < D; ++i) {
    double upstar = up[i]+dt/mass*(contact[i]+massreduced*problem->g[i]-vol*dp[i]);
    double drag = gammastar*(upstar-u[i]/c);
    problem->particle_force[ip*D+i] = -drag-vol*dp[i];
    f[U+i] = -drag;//-vol*dp[i];
  }
  *dfdu = gammastar/c;
  *dfddp = gammastar*dt/mass*vol;//-vol;
}

static void fluid_problem_f(const FluidProblem *problem, const double *sol, double *dsol, const double *solold, const double *dsolold, const double* mesh_velocity, const double c, const double *dc, const double *bf, const double cold, const double dus[D][D], const double rho, const double mu, double dt, int iel, double *f0, double *f1, double *f00, double *f10, double *f01, double *f11) {
  double p = sol[P];
  double taup = problem->taup[iel];
  double tauc = problem->tauc[iel];
  const int n_fields = fluid_problem_n_fields(problem);
  double u[D], uold[D], du[D][D], duold[D][D], umesh[D], dp[D];
  if(problem->stab_param != 0)
    taup = tauc = 0;
  double divu = 0;
  double divus = 0;
  double nold = 0;
  for (int iD = 0; iD < D; ++iD) {
    u[iD] = sol[U+iD];
    uold[iD] = solold[U+iD];
    nold += uold[iD]*uold[iD];
    umesh[iD] = mesh_velocity[iD];
    dp[iD] = dsol[P*D+iD];
    for (int jD = 0; jD < D; ++jD) {
      du[iD][jD] = dsol[(U+iD)*D+jD];
      duold[iD][jD] = dsolold[(U+iD)*D+jD];
    }
    divu += du[iD][iD];
  }
  if (problem->usolid) {
    for (int iD = 0; iD < D; ++iD)
      divus += dus[iD][iD];
  }
  else {
    divus = (c-cold)/dt;
  }
  nold = sqrt(nold);
  f0[P] = divus;
  double *g = problem->g;
  double rhoreduced = (problem->reduced_gravity ? (rho-problem->rho[0]) : rho);
  double drag = problem->volume_drag + problem->quadratic_drag*nold;
  for (int i = 0; i < D; ++i) {
    f0[U+i] =
      + c*(dp[i]- g[i]*rhoreduced- bf[i])
      + drag*u[i];
      //+ 5.3e5*u[i];
    f00[(U+i)*n_fields+U+i] = drag;//5.3e5;
    if(problem->temporal) {
      f0[U+i] += rho*(u[i]-uold[i])/dt;
      f00[(U+i)*n_fields+U+i] += rho/dt;
    }
    f01[(U+i)*n_fields*D+P*D+i] = c;
    if (problem->advection) {
      // advection :
      // dj(uj ui /c) = uj/c dj(ui) + ui/c dj(uj) - uj ui /(c*c) dj(c)
      //              = ui/c * (dj(uj)- uj/c*dj(c)) + uj/c*dj(ui)
      for (int j = 0; j < D; ++j) {
        f0[U+i] += rho/c*(uold[j]*du[i][j] + u[i]*(duold[j][j]-uold[j]/c*dc[j]));
        f0[U+i] -= rho/c*(umesh[j]*du[i][j]);
        f00[(U+i)*n_fields+U+i] += rho/c*(duold[j][j]-uold[j]/c*dc[j]);
        f01[(U+i)*n_fields*D+(U+i)*D+j] += rho/c*uold[j];
        f01[(U+i)*n_fields*D+(U+i)*D+j] -= rho/c*umesh[j];
      }
    }
    for (int j = 0; j < D; ++j) {
      f1[(U+i)*D+j] = mu*(du[i][j]-u[i]*dc[j]/c+ du[j][i]-u[j]*dc[i]/c);
      f10[((U+i)*D+j)*n_fields+U+j] += -mu*dc[i]/c;
      f10[((U+i)*D+j)*n_fields+U+i] += -mu*dc[j]/c;
      f11[((U+i)*D+j)*n_fields*D+(U+j)*D+i] += mu;
      f11[((U+i)*D+j)*n_fields*D+(U+i)*D+j] += mu;
      // SUPG
      double supg = uold[j]*taup;
      f1[(U+i)*D+j] += supg*f0[U+i];
      f10[((U+i)*D+j)*n_fields+U+i] += supg*f00[(U+i)*n_fields+U+i];
      f11[((U+i)*D+j)*(n_fields*D)+P*D+i] += supg*f01[(U+i)*n_fields*D+P*D+i];
      for (int k = 0; k < D; ++k)
        f11[((U+i)*D+j)*n_fields*D+(U+i)*D+k] += supg*f01[(U+i)*n_fields*D+(U+i)*D+k];
    }
    double lsic = tauc*rho;
    f1[(U+i)*D+i] += (divus+divu)*lsic;
    for (int j = 0; j < D; ++j) {
      f11[((U+i)*D+i)*(n_fields*D)+(U+j)*D+j] += lsic;
    }
    f1[P*D+i] = -u[i];
    f10[(P*D+i)*n_fields+U+i] += -1;
    // PSPG
    double pspg = taup/rho;
    f1[P*D+i] += pspg*(f0[U+i]+(problem->drag_in_stab?0:(1-c)*dp[i])) + problem->stab_param*dp[i];
    f11[(P*D+i)*(n_fields*D)+P*D+i] += pspg*(f01[(U+i)*n_fields*D+P*D+i]+(problem->drag_in_stab?0:1-c))+problem->stab_param;
    f10[(P*D+i)*n_fields+U+i] += pspg*f00[(U+i)*n_fields+(U+i)];
    for (int j = 0; j < D; ++j) {
      f11[(P*D+i)*n_fields*D+(U+i)*D+j] = pspg*f01[(U+i)*n_fields*D+(U+i)*D+j];
    }
  }
}

static void compute_weak_boundary_conditions(FluidProblem *problem, const double *solution_old, 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 *dsold = malloc(sizeof(double)*n_fields*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *f00 = malloc(sizeof(double)*n_fields*n_fields);
  double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
  double i_bnd = 0;
  double iv_bnd_up = 0;
  double iv_bnd_down = 0;
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    int n_value = weak_boundary_n_values(wbnd);
    double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
    weak_boundary_values(mesh, bnd, wbnd, data);

    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      grad_shape_functions(dxidx, dphi);
      double dc[D] = {0};
      for (int j = 0; j < n_fields*D; ++j) {
        ds[j] = 0;
      }
      for (int i = 0; i < N_SF; ++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];
          for (int k = 0; k < D; ++k){
            ds[j*D+k] += dphi[i][k]*dof;
            dsold[j*D+k] += dphi[i][k]*dofold;
          }
        }
        for (int k=0; k<D; ++k){
          dc[k] += problem->porosity[el[i]]*dphi[i][k];
        }
      }
      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
        double dp[D] = {0};
        double umesh[D] = {0};
        double *local_vector = &all_local_vector[local_size*iel];
        double *local_matrix = &all_local_matrix[local_size*local_size*iel];
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
          sold[i] = 0;
          f0[i] = 0;
        }
        for (int i = 0; i < n_fields*n_fields; ++i) {
          f00[i] = 0;
        }
        for (int i = 0; i < n_fields*n_fields*D; ++i) {
          f01[i] = 0;
        }
        double c = 0;
        double a = 0;
        for (int i = 0; i < DIMENSION; ++i) {
          c += problem->porosity[nodes[i]]*phi[i];
          if (problem->n_fluids == 2) {
            a += problem->concentration[iel*N_N+cl[i]]*phi[i];
          }
          for (int j = 0; j < n_fields; ++j) {
            double dof = solution[nodes[i]*n_fields+j];
            double dofold = solution_old[nodes[i]*n_fields+j];
            s[j] += phi[i]*dof;
            sold[j] += phi[i]*dofold;
          }
          for(int j = 0; j < D; ++j){
            double dofmesh = problem->mesh_velocity[nodes[i]*D+j];
            umesh[j] += phi[i]*dofmesh;
          }
        }
        double rho, mu;
        fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
        const double jw = LQW[iqp]*detbnd;
        if (wbnd->type != BND_WALL){
          for (int i = 0; i < D; ++i){
            if (wbnd->vid<0) {
              i_bnd -= a*s[U+i]*n[i]*jw*dt;
              iv_bnd_up -= s[U+i]*n[i]*jw;
            }
            else {
              iv_bnd_down -= dataqp[wbnd->vid+i]*n[i]*jw;
              i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
            }
          }
        }
        f_boundary(wbnd,problem,n,f0,s,ds,sold,dsold,umesh,c,dc,rho,mu,dt,iel,dataqp,f00,f01);
        for (int ifield = 0; ifield < n_fields; ++ifield) {
          for (int iphi = 0; iphi < D; ++iphi){
            if (ifield<D){
              problem->boundary_force[bndid*D + ifield] += phi[iphi]*f0[ifield]*jw;
            }
            local_vector[cl[iphi]+N_SF*ifield] += phi[iphi]*f0[ifield]*jw;
          }
        }
        for (int jfield = 0; jfield < n_fields; ++jfield) {
          for (int ifield = 0; ifield < n_fields; ++ifield){
            for (int iphi = 0; iphi < N_LSF; ++iphi){
              for (int jphi = 0; jphi < N_LSF; ++jphi){
                double d = phi[jphi]*f00[ifield*n_fields+jfield];
                LOCAL_MATRIX(cl[iphi],cl[jphi],ifield,jfield) += jw*phi[iphi]*d;
              }
              for (int jphi = 0; jphi < N_SF; ++jphi){
                double d = 0;
                for (int jD = 0; jD < D; ++jD) {
                  d += dphi[jphi][jD]*f01[ifield*n_fields*D+jfield*D+jD];
                }
                LOCAL_MATRIX(cl[iphi],jphi,ifield,jfield) += jw*phi[iphi]*d;
              }
            }
          }
        }
      }
    }
    free(data);
  }
  free(s);
  free(ds);
  free(sold);
  free(dsold);
  free(f0);
  free(f00);
}

double fluid_problem_a_integ_volume(FluidProblem *problem)
{
  const Mesh *mesh = problem->mesh;
  double I_a = 0;
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D];
    const double detj = mesh_dxidx(mesh,iel,dxidx);
    for (int i = 0; i < N_N; ++i)
      for (int j = 0; j< N_N; ++j)
        I_a += detj * problem->porosity[el[i]]*problem->concentration[iel*N_N+j]*mass_matrix[i][j];
  }
  return I_a;
}

double fluid_problem_volume_flux(FluidProblem *problem, const char *tagname)
{
  const Mesh *mesh = problem->mesh;
  const double *solution = problem->solution;
  double Q = 0;
  double wtotal = 0;
  int n_fields = D+1;
  int found = 0;
  for (int i = 0; i < N_LQP; ++i) wtotal += LQW[i];
  for (int ibnd = 0; ibnd < mesh->n_boundaries; ++ibnd) {
    if (strcmp(mesh->boundary_names[ibnd],tagname) != 0)
      continue;
    found = 1;
    MeshBoundary *bnd = &problem->boundaries[ibnd];
    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      double vnmean = 0;
      for (int i = 0; i < D; ++i) {
        for (int iD = 0; iD < D; ++iD) {
          vnmean += n[iD]*problem->solution[nodes[i]*n_fields+U+iD];
        }
      }
      Q += vnmean*detbnd*wtotal/D;
    }
  }
  if(found == 0) {
    printf("boundary '%s' not found\n", tagname);
    exit(1);
  }
  printf("Q = %g wtotal = %g\n", Q,wtotal);
  return Q;
}

double fluid_problem_a_integ_bnd(FluidProblem *problem, double dt)
{
  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 i_bnd = 0;
  for (int ibnd = 0; ibnd < problem->n_weak_boundaries; ++ibnd){
    WeakBoundary *wbnd = problem->weak_boundaries + ibnd;
    int bndid= -1;
    for (int i = 0; i < mesh->n_boundaries; ++i) {
      if (strcmp(mesh->boundary_names[i],wbnd->tag) == 0){
        bndid = i;
      }
    }
    if (bndid == -1) {
      printf("Mesh has no boundary with name \"%s\".", wbnd->tag);
    }
    MeshBoundary *bnd = &problem->boundaries[bndid];
    int n_value = weak_boundary_n_values(wbnd);
    double *data = malloc(n_value*bnd->n_interfaces*N_LQP*sizeof(double));
    weak_boundary_values(mesh, bnd, wbnd, data);

    for (int iinterface = 0; iinterface < bnd->n_interfaces; ++iinterface) {
      const int *interface = &mesh->interfaces[bnd->interfaces[iinterface]*4];
      const int iel = interface[0];
      const int icl = interface[1];
      const int *cl = elbnd[icl];
      int nodes[N_LSF];
      double *x[N_LSF];
      const int *el = &mesh->elements[iel*N_N];
      for (int i = 0; i < D; ++i){
        nodes[i] = el[cl[i]];
        x[i] = &mesh->x[nodes[i]*3];
      }
      double dxidx[D][D], dphi[N_N][D];
      mesh_dxidx(mesh,iel,dxidx);
      double n[D],detbnd;
      get_normal_and_det(x,n,&detbnd);
      grad_shape_functions(dxidx, dphi);

      for (int iqp = 0; iqp < N_LQP; ++iqp) {
        double *dataqp = &data[n_value*(N_LQP*iinterface+iqp)];
        double phi[DIMENSION];
        l_shape_functions(LQP[iqp],phi);
        for (int i = 0; i < n_fields; ++i) {
          s[i] = 0;
        }
        double a = 0;
        for (int i = 0; i < DIMENSION; ++i) {
          if (problem->n_fluids == 2) {
            a += problem->concentration[iel*N_N+cl[i]]*phi[i];
          }
          for (int j = 0; j < n_fields; ++j) {
            double dof = solution[nodes[i]*n_fields+j];
            s[j] += phi[i]*dof;
          }
        }
        const double jw = LQW[iqp]*detbnd;
        if (wbnd->type != BND_WALL){
          for (int i = 0; i < D; ++i){
            if (wbnd->vid<0) {
              i_bnd -= a*s[U+i]*n[i]*jw*dt;
            }
            else {
              i_bnd -= a*dataqp[wbnd->vid+i]*n[i]*jw*dt;
            }
          }
        }
      }
    }
  }
  return i_bnd;
}

static void fluid_problem_volume(FluidProblem *problem, const double *solution_old, const double *u_solid, double dt, double *all_local_vector, double *all_local_matrix)
{
  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 *dsold = malloc(sizeof(double)*n_fields*D);
  double *mesh_velocity = malloc(sizeof(double)*D);
  double *f0 = malloc(sizeof(double)*n_fields);
  double *f1 = malloc(sizeof(double)*n_fields*D);
  double *f00 = malloc(sizeof(double)*n_fields*n_fields);
  double *f10 = malloc(sizeof(double)*n_fields*n_fields*D);
  double *f01 = malloc(sizeof(double)*n_fields*n_fields*D);
  double *f11 = malloc(sizeof(double)*n_fields*n_fields*D*D);
  problem->kinetic_energy = 0;
  for (int iel=0; iel < mesh->n_elements; ++iel) {
    const int *el = &mesh->elements[iel*N_N];
    double dxidx[D][D], dphi[N_N][D];
    double dus[D][D] = {0};
    const double detj = mesh_dxidx(mesh,iel,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) {
      for (int j = 0; j < D; ++j) {
        ds[i*D+j] = 0;
        dsold[i*D+j] = 0;
      }
    }
    double dc[D] = {0}, da[D]={0};
    for (int i = 0; i < N_SF; ++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];
        for (int k = 0; k < D; ++k) {
          ds[j*D+k] += dphi[i][k]*dof;
          dsold[j*D+k] += dphi[i][k]*dofold;
        }
      }
      for (int k=0; k<D; ++k){
        dc[k] += problem->porosity[el[i]]*dphi[i][k];
        for (int j = 0; j < D; ++j) {
          dus[j][k] += u_solid[el[i]*DIMENSION+j]*dphi[i][k];
        }
      }
      if (problem->n_fluids==2) {
        for (int k=0; k<D; ++k){
          da[k] += problem->concentration[iel*N_N+i]*dphi[i][k];
        }
      }
    }
    for (int iqp = 0; iqp< N_QP; ++iqp) {
      double phi[N_SF];
      shape_functions(QP[iqp],phi);
      double bf[D];
      for (int j = 0; j < D; ++j) {
        bf[j] = 0;
        mesh_velocity[j] = 0;
      }
      for (int i = 0; i < n_fields; ++i) {
        s[i] = 0;
        sold[i] = 0;
      }
      for (int i = 0; i < n_fields*n_fields; ++i) {
        f00[i] = 0;
      }
      for (int i = 0; i < D*n_fields*n_fields; ++i) {
        f10[i] = 0;
        f01[i] = 0;
      }
      for (int i = 0; i < D*D*n_fields*n_fields; ++i) {
        f11[i] = 0;
      }
      double c = 0, a = 0;
      double cold = 0, aold = 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;
        }
        if (problem->n_fluids==2) {
          a += problem->concentration[iel*N_N+i]*phi[i];
        }
        for (int j=0; j<D; ++j) {
          mesh_velocity[j] += phi[i]*problem->mesh_velocity[el[i]*D+j];
          bf[j] += phi[i]*problem->bulk_force[el[i]*D+j];
        }
      }
      double mu;
      double rho;
      fluid_problem_interpolate_rho_and_nu(problem,a, &rho, &mu);
      const double jw = QW[iqp]*detj;
      problem->kinetic_energy += rho*sqrt(s[U]*s[U]+s[U+1]*s[U+1])*sqrt(s[U]*s[U]+s[U+1]*s[U+1])/2*jw;

      fluid_problem_f(problem,s,ds,sold,dsold,mesh_velocity,c,dc,bf,cold,dus,rho,mu,dt,iel,f0,f1,f00,f10,f01,f11);
      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) {
        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 m = jw*phi[iphi]*phi[jphi]*f00[ifield*n_fields+jfield];
              for (int id = 0; id < D; ++id) {
                m += jw*dphi[iphi][id]*phi[jphi]*f10[(ifield*D+id)*n_fields+jfield];
                for (int jd = 0; jd < D; ++jd) {
                  m += jw*dphi[iphi][id]*dphi[jphi][jd]*f11[(ifield*D+id)*n_fields*D+jfield*D+jd];
                }
              }
              for (int jd = 0; jd < D; ++jd) {