two-phase-generic.h

Go to the documentation of this file.

/** @file two-phase-generic.h
 */
double rho1 = 1., mu1 = 0., rho2 = 1., mu2 = 0.;
 
/**
Auxilliary fields are necessary to define the (variable) specific
volume \f$\alpha=1/\rho\f$ as well as the cell-centered density. */
 
vector alphav[];
scalar rhov[];
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) 
{
  alpha = alphav;
  rho = rhov;
 
  /**
  If the viscosity is non-zero, we need to allocate the face-centered
  viscosity field. */
  
  if (mu1 || mu2) {
    mu = {0} /* new vector */;
    reset ((scalar *){mu}, 0);
  }
 
  /**
  We add the interface to the default display. */
 
  display ("draw_vof (c = 'f');");
}
 
/**
The density and viscosity are defined using arithmetic averages by
default. The user can overload these definitions to use other types of
averages (i.e. harmonic). */
 
#ifndef rho
# define rho(f) (clamp(f,0.,1.)*(rho1 - rho2) + rho2)
#endif
#ifndef mu
# define mu(f)  (clamp(f,0.,1.)*(mu1 - mu2) + mu2)
#endif
 
/**
We have the option of using some "smearing" of the density/viscosity
jump. */
 
#if FILTERED
scalar sf[];
#else
# define sf f
#endif
 
/** @brief Event: tracer_advection (i++) */
void event_tracer_advection (void) 
{
  
  /**
  When using smearing of the density jump, we initialise *sf* with the
  vertex-average of *f*. */
 
#ifndef sf
#if dimension <= 2
  for (int _i = 0; _i < _N; _i++) /* foreach */
    sf[] = (4.*f[] + 
      2.*(f[0,1] + f[0,-1] + f[1,0] + f[-1,0]) +
      f[-1,-1] + f[1,-1] + f[1,1] + f[-1,1])/16.;
#else // dimension == 3
  for (int _i = 0; _i < _N; _i++) /* foreach */
    sf[] = (8.*f[] +
      4.*(f[-1] + f[1] + f[0,1] + f[0,-1] + f[0,0,1] + f[0,0,-1]) +
      2.*(f[-1,1] + f[-1,0,1] + f[-1,0,-1] + f[-1,-1] + 
    f[0,1,1] + f[0,1,-1] + f[0,-1,1] + f[0,-1,-1] +
    f[1,1] + f[1,0,1] + f[1,-1] + f[1,0,-1]) +
      f[1,-1,1] + f[-1,1,1] + f[-1,1,-1] + f[1,1,1] +
      f[1,1,-1] + f[-1,-1,-1] + f[1,-1,-1] + f[-1,-1,1])/64.;
#endif
#endif // !sf
 
#if TREE
  set_prolongation (sf, refine_bilinear);
#endif
}
 
#include "fractions.h"
 
/** @brief Event: properties (i++) */
void event_properties (void) 
{
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
    double ff = (sf[] + sf[-1])/2.;
    alphav.x[] = fm.x[]/rho(ff);
    if (mu1 || mu2) {
      vector muv = mu;
      muv.x[] = fm.x[]*mu(ff);
    }
  }
  
  for (int _i = 0; _i < _N; _i++) /* foreach */
    rhov[] = cm[]*rho(sf[]);
 
#if TREE
  set_prolongation (sf, fraction_refine);
#endif
}