implicit.h

Go to the documentation of this file.

/** @file implicit.h
 */
/**
# Ohmic conduction
 
This function approximates implicitly the EHD equation set given by
the electric Poisson equation and the ohmic conduction term of the
charge density conservation (the advection term is computed elsewhere
using a tracer advection scheme),
\f[
\nabla \cdot( \epsilon \nabla  \phi^{n+1}) = -\rho_e^{n+1} \quad \mathrm{and}
\quad (\rho_e^{n+1}-\rho_e^n) = \Delta t\nabla \cdot (K\nabla \phi^{n+1})
\f] 
where \f$\rho_e\f$ is the charge density, and \f$\phi\f$ is the electric potential.
\f$K\f$ and \f$\epsilon\f$ are the conductivity and permittivity respectively.
 
Substituting the Poisson equation into the conservation equation gives
the following time-implicit scheme,
\f[
\nabla \cdot [(K \, \Delta t + \epsilon) \nabla \phi^{n+1}] = -\rho_e^{n}
\f]
 
We need the Poisson solver and the timestep *dt*. */
 
#include "poisson.h"
extern double dt;
 
/**
The electric potential and the volume charge density are scalars while
the permittivity and conductivity are face vectors. In *mgphi* we will
store the statistics for the multigrid resolution of the electric
poisson equation. */
 
scalar phi[], rhoe[];
vector epsilon[], K[];
mgstats mgphi;
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) 
{
 
  /**
  The restriction/refine attributes of the charge density are those of a tracer
  otherwise the conservation is not guaranteed. */
 
#if TREE
  rhoe.restriction = restriction_volume_average;
  rhoe.refine = refine_linear;
#endif
 
  /**
  By default the permittivity is unity and other quantities are zero. */
 
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    epsilon.x[] = K.x[] = fm.x[];
}
 
/** @brief Event: tracer_diffusion (i++) */
void event_tracer_diffusion (void) 
{
  scalar rhs[];
  
  /**
  The r.h.s of the poisson equation is set to \f$-\rho_e^n\f$, */ 
  
  for (int _i = 0; _i < _N; _i++) /* foreach */
    rhs[] = - rhoe[]*cm[];
 
  if (K.x.i) {
    
    /**
    We compute the coefficients of the Laplacian: $(K dt +\epsilon)$. */
 
    vector f[];
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      f.x[] = K.x[]*dt + epsilon.x[];
 
    /**
    The poisson equation is solved to obtain \f$\phi^{n+1}\f$. */
 
    mgphi = poisson (phi, rhs, f);
 
    /**
    Finally \f$\rho_e^{n+1}\f$ is computed as
    \f$\rho_e^{n+1} = -\nabla \cdot (\epsilon \nabla \phi^{n+1})\f$. */
 
#if TREE
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      f.x[] = epsilon.x[]*(phi[] - phi[-1])/Delta;
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      rhoe[] = 0.;
      for (int _d = 0; _d < dimension; _d++)
  rhoe[] -= f.x[1] - f.x[];
      rhoe[] /= cm[]*Delta;
    }
#else // Cartesian
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      rhoe[] = 0.;
      for (int _d = 0; _d < dimension; _d++)
  rhoe[] -= (epsilon.x[1]*(phi[1] - phi[]) -
       epsilon.x[]*(phi[] - phi[-1]));
      rhoe[] /= cm[]*sq(Delta);
    }
#endif
  }
 
  /**
  In the absence of conductivity, the electric potential only varies 
  if the electrical permittivity varies, */
  
  else
    mgphi = poisson (phi, rhs, epsilon);
}