ParabolicFunctionNC1.5Stage.c File Reference

Evaluate the parabolic term using a "1.5"-stage discretization. More...

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
int	ParabolicFunctionNC1_5Stage (double *par, double *u, void *s, void *m, double t)

Detailed Description

Evaluate the parabolic term using a "1.5"-stage discretization.

Author

Debojyoti Ghosh

Definition in file ParabolicFunctionNC1.5Stage.c.

Function Documentation

ParabolicFunctionNC1_5Stage()

int ParabolicFunctionNC1_5Stage	(	double *	par,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Evaluate the parabolic term using a "1.5"-stage finite-difference spatial discretization: The parabolic term is assumed to be of the form:

\begin{equation} {\bf P}\left({\bf u}\right) = \sum_{d1=0}^{D-1}\sum_{d2=0}^{D-1} \frac {\partial^2 h_{d1,d2}\left(\bf u\right)} {\partial x_{d1} \partial x_{d2}}, \end{equation}

where \(d1\) and \(d2\) are spatial dimension indices, and \(D\) is the total number of spatial dimensions (HyPar::ndims). This term is discretized at a grid point as:

\begin{equation} \left.{\bf P}\left({\bf u}\right)\right|_j = \sum_{d1=0}^{D-1} \sum_{d2=0}^{D-1} \frac { \mathcal{D}_{d1}\mathcal{D}_{d2} \left[ {\bf h}_{d1,d2} \right] } {\Delta x_{d1} \Delta x_{d2}}, \end{equation}

where \(\mathcal{D}\) denotes the finite-difference approximation to the first derivative.

• When \(d = d1 = d2\), the operator \(\mathcal{D}^2_d\left[\cdot\right] = \mathcal{D}_{d}\mathcal{D}_{d}\left[\cdot\right]\) is computed in one step as the finite-difference approximation to the second derivative (Laplacian) \(\mathcal{L}_d\left[\cdot\right]\) using HyPar::SecondDerivativePar.
• When \(d1 \ne d2 \), each of the first derivative approximations are \(\mathcal{D}_{d1}\) and \(\mathcal{D}_{d2}\) are computed separately, and thus the cross-derivative is evaluated in two steps using HyPar::FirstDerivativePar.

Note: this form of the parabolic term does allow for cross-derivatives ( \( d1 \ne d2 \)).

To use this form of the parabolic term:

• specify "par_space_type" in solver.inp as "nonconservative-1.5stage" (HyPar::spatial_type_par).
• the physical model must specify \({\bf h}_{d1,d2}\left({\bf u}\right)\) through HyPar::HFunction.

See also

ParabolicFunctionNC2Stage()

Parameters

par	array to hold the computed parabolic term
u	solution
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current simulation time

Definition at line 35 of file ParabolicFunctionNC1.5Stage.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  double        *Func   = solver->fluxC;
  double        *Deriv1 = solver->Deriv1;
  double        *Deriv2 = solver->Deriv2;
  int           d, d1, d2, v, p, done;
  double        dxinv1, dxinv2;
  _DECLARE_IERR_;

  int     ndims   = solver->ndims;
  int     nvars   = solver->nvars;
  int     ghosts  = solver->ghosts;
  int     *dim    = solver->dim_local;
  double  *dxinv  = solver->dxinv;
  int     size   = solver->npoints_local_wghosts;

  if (!solver->HFunction) return(0); /* zero parabolic terms */
  solver->count_par++;

  int index[ndims];
  _ArraySetValue_(par,size*nvars,0.0);

  for (d1 = 0; d1 < ndims; d1++) {
    for (d2 = 0; d2 < ndims; d2++) {

      /* calculate the diffusion function */
      IERR solver->HFunction(Func,u,d1,d2,solver,t);                      CHECKERR(ierr);
      if (d1 == d2) {
        /* second derivative with respect to the same independent variable */
        IERR solver->SecondDerivativePar(Deriv2,Func,d1,solver,mpi);      CHECKERR(ierr);
      } else {
        /* cross derivative */
        IERR solver->FirstDerivativePar(Deriv1,Func  ,d1, 1,solver,mpi);  CHECKERR(ierr);
        IERR MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,Deriv1);  CHECKERR(ierr);
        IERR solver->FirstDerivativePar(Deriv2,Deriv1,d2,-1,solver,mpi);  CHECKERR(ierr);
      }

      /* calculate the final term - second derivative of the diffusion function */
      done = 0; _ArraySetValue_(index,ndims,0);
      while (!done) {
        _ArrayIndex1D_(ndims,dim,index,ghosts,p);
        _GetCoordinate_(d1,index[d1],dim,ghosts,dxinv,dxinv1);
        _GetCoordinate_(d2,index[d2],dim,ghosts,dxinv,dxinv2);
        for (v=0; v<nvars; v++) par[nvars*p+v] += (dxinv1*dxinv2 * Deriv2[nvars*p+v]);
        _ArrayIncrementIndex_(ndims,dim,index,done);
      }

    }
  }

  if (solver->flag_ib) _ArrayBlockMultiply_(par,solver->iblank,size,nvars);
  return(0);
}