Carmen.h File Reference

The .h that includes all functions headers. More...

#include "carmen.pre"
#include "PreProcessor.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include "Vector.h"
#include "Matrix.h"
#include "Timer.h"
#include "Parameters.h"
#include "PrintGrid.h"
#include "TimeAverageGrid.h"
#include "Cell.h"
#include "Node.h"
#include "FineMesh.h"

Include dependency graph for Carmen.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros
#define	Max(x, y)   (((x) > (y)) ? (x):(y))
#define	Max3(x, y, z)   (Max((x),Max((y),(z))))
#define	Min(x, y)   (((x) < (y)) ? (x):(y))
#define	Min3(x, y, z)   (Min((x),Min((y),(z))))
#define	power2(x)   ((x)*(x))
#define	power3(x)   ((x)*(x)*(x))
#define	Abs(x)   ( ((x) < 0)? -(x):(x) )

Functions
void	AdaptTimeStep ()
	Adapts time step when required. More...
int	BC (int i, int AxisNo, int N=(1<< ScaleNb))
	Returns the position of i taking into account the boundary conditions in the direction AxisNo. The number of points in this direction is N. Example: for AxisNo=1 and for N=256, i must be between 0 and 255. If i=-1, the function returns 255 for periodic boundary conditions and 0 for Neumann boundary conditions. More...
int	BoundaryRegion (const Vector &X)
	Returns the boundary region type at the position X=(x,y,z). The returned value correspond to: 0 = Fluid (not in the boundary) 1 = Inflow 2 = Outflow 3 = Solid with free-slip walls 4 = Solid with adiabatic walls 5 = Solid with isothermal walls. More...
void	Backup (Node *Root)
	Stores the tree structure and data in order to restart a multiresolution computation. More...
void	Backup (FineMesh *Root)
	Stores the data countained in a regular mesh Root in order to restart a finite volume computation. More...
double	CPUTimeRef (int iterations, int scales)
	Returns the time required by a finite volume computation using iterations iterations and scales scales. It is use to estimate the CPU time compression. More...
real	ComputedTolerance (const int ScaleNo)
	Returns the computed tolerance at the scale ScaleNo, either using Harten or Donoho thresholding (if CVS=true). More...
int	DigitNumber (int arg)
	Returns the number of digits of the integer arg. More...
int	FileWrite (FILE *f, const char *format, real arg)
	Writes in binary or ASCII mode the real number arg into the file f with the format format. The global parameter DataIsBinary determines this choice. More...
Vector	Flux (Cell &Cell1, Cell &Cell2, Cell &Cell3, Cell &Cell4, int AxisNo)
	Returns the flux at the interface between Cell2 and Cell3. Here a 4-point space scheme is used. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides. More...
void	GetBoundaryCells (Cell &Cell1, Cell &Cell2, Cell &Cell3, Cell &Cell4, Cell &C1, Cell &C2, Cell &C3, Cell &C4, int AxisNo)
	Transform the 4 cells of the flux Cell1, Cell2, Cell3, Cell4 into C1, C2, C3, C4, to take into account boundary conditions (used in Flux.cpp). More...
Vector	InitAverage (real x, real y=0., real z=0.)
	Returns the initial condition in (x, y, z) form the one defined in carmen.ini. More...
real	InitResistivity (real x, real y=0., real z=0.)
	Returns the initial resistivity condition in (x, y, z) form the one defined in carmen.eta. More...
void	InitParameters ()
	Inits parameters from file carmen.par. If a parameter is not mentioned in this file, the default value is used. More...
void	InitTimeStep ()
	Inits time step and all the parameters which depend on it. More...
void	InitTree (Node *Root)
	Inits tree structure from initial condition, starting form the node Root. Only for multiresolution computations. More...
Vector	Limiter (const Vector u, const Vector v)
	Returns the value of the slope limiter between the slopes u and v. More...
real	Limiter (const real x)
	Returns the valur of slope limiter from a real value x. More...
real	MinAbs (const real a, const real b)
	Returns the minimum in module of a and b. More...
real	NormMaxQuantities (const Vector &V)
	Returns the Max-norm of the vector where every quantity is divided by its characteristic value. More...
void	Performance (const char *FileName)
	Computes the performance of the computation and, for cluster computations, write it into file FileName. More...
void	PrintIntegral (const char *FileName)
	Writes the integral values, like e.g flame velocity, global error, into file FileName. More...
void	RefreshTree (Node *Root)
	Refresh the tree structure, i.e. compute the cell-averages of the internal nodes by projection and those of the virtual leaves by prediction. The root node is Root. Only for multiresolution computations. More...
void	Remesh (Node *Root)
	Remesh the tree structure after a time evolution. The root node is Root. Only for multiresolution computations. More...
Vector	FluxX (const Vector &Avg)
	Returns the physical flux of MHD equations in X direction. More...
Vector	FluxY (const Vector &Avg)
	Returns the physical flux of MHD equations in Y direction. More...
Vector	FluxZ (const Vector &Avg)
	Returns the physical flux of MHD equations in Z direction. More...
Vector	SchemeHLL (const Cell &Cell1, const Cell &Cell2, const Cell &Cell3, const Cell &Cell4, const int AxisNo)
	Returns the HLL numerical flux for MHD equations. The scheme uses four cells to estimate the flux at the interface. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides. More...
Vector	SchemeHLLD (const Cell &Cell1, const Cell &Cell2, const Cell &Cell3, const Cell &Cell4, const int AxisNo)
	Returns the HLLD numerical flux for MHD equations. The scheme uses four cells to estimate the flux at the interface. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides. More...
Matrix	stateUstar (const Vector &AvgL, const Vector &AvgR, const real prel, const real prer, real &slopeLeft, real &slopeRight, real &slopeM, real &slopeLeftStar, real &slopeRightStar, int AxisNo)
	Returns the intermediary states of HLLD numerical flux for MHD equations. More...
void	fluxCorrection (Vector &Flux, const Vector &AvgL, const Vector &AvgR, int AxisNo)
	This function apply the divergence-free correction to the numerical flux. More...
void	ShowTime (Timer arg)
	Writes on screen the estimation of total and remaining CPU times. These informations are stored in the timer arg. More...
int	Sign (const real a)
	Returns 1 if a is non-negative, -1 elsewhere. More...
Vector	ArtificialViscosity (const Vector &Cell1, const Vector &Cell2, real dx, int AxisNo)
	Returns the artificial diffusion source terms in the cell UserCell. More...
Vector	ResistiveTerms (Cell &Cell1, Cell &Cell2, Cell &Cell3, Cell &Cell4, int AxisNo)
	Returns the resistive source terms in the cell UserCell. More...
Vector	Source (Cell &UserCell)
	Returns the source term in the cell UserCell. More...
real	Step (real x)
	Returns a step (1 if x <0, 0 if x >0, 0.5 if x=0) More...
void	TimeEvolution (FineMesh *Root)
	Computes a time evolution on the regular fine mesh Root. Only for finite volume computations. More...
void	TimeEvolution (Node *Root)
	Computes a time evolution on the tree structure, the root node being Root. Only for multiresolution computations. More...
void	View (FineMesh *Root, const char *AverageFileName)
	Writes the current cel–averages of the fine mesh Root into file AverageFileName. Only for finite volume computations. More...
void	View (Node *Root, const char *TreeFileName, const char *MeshFileName, const char *AverageFileName)
	Writes the data of the tree structure into files TreeFileName (tree structure), MeshFileName (mesh) and AverageFileName (cell-averages). The root node is Root. Only for multiresolution computations. More...
void	ViewEvery (FineMesh *Root, int arg)
	Same as previous for a fine mesh Root. Only for finite volume computations. More...
void	ViewEvery (Node *Root, int arg)
	Writes into file the data of the tree structure at iteration arg. The output file names are AverageNNN.dat and MeshNNN.dat, NNN being the iteration in an accurate format. The root node is Root. Only for multiresolution computations. More...
void	ViewIteration (FineMesh *Root)
	Same as previous for a fine mesh Root. Only for finite volume computations. More...
void	ViewIteration (Node *Root)
	Writes into file the data of the tree structure from physical time PrintTime1 to physical time PrintTime6. The output file names are Average_N.dat and Mesh_N.dat, N being between 1 and 6. The root node is Root. Only for multiresolution computations. More...
void	CreateMPIType (FineMesh *Root)
void	CPUExchange (FineMesh *Root, int)
	Parallel function DOES NOT WORK! More...
void	FreeMPIType ()
	Parallel function DOES NOT WORK! More...
void	CreateMPITopology ()
	Parallel function DOES NOT WORK! More...
void	CreateMPILinks ()
	Parallel function DOES NOT WORK! More...
void	ReduceIntegralValues ()
	Parallel function DOES NOT WORK! More...

Detailed Description

The .h that includes all functions headers.

Macro Definition Documentation

#define Abs

(

x

)

( ((x) < 0)? -(x):(x) )

Returns the absolute value of x.

#define Max	(		x,
			y
	)		(((x) > (y)) ? (x):(y))

Returns the Maximum value between x and y.

#define Max3	(		x,
			y,
			z
	)		(Max((x),Max((y),(z))))

Returns the Maximum value between x, y and z.

#define Min	(		x,
			y
	)		(((x) < (y)) ? (x):(y))

Returns the mininum value between x and y.

#define Min3	(		x,
			y,
			z
	)		(Min((x),Min((y),(z))))

Returns the minimum value between x, y and z.

#define power2

(

x

)

((x)*(x))

Returns the square of x.

#define power3

(

x

)

((x)*(x)*(x))

Returns the cube of x.

Function Documentation

void AdaptTimeStep

(

)

Adapts time step when required.

Returns

void

{
    int     RemainingIterations;
    real    RemainingTime;


    // Security : do nothing if ConstantTimeStep is true

    if (ConstantTimeStep)
        return;

    // Compute remaining time

    RemainingTime = PhysicalTime-ElapsedTime;


    // In this case, use time adaptivity based on CFL
    if(Resistivity)
        TimeStep = CFL*min(SpaceStep/Eigenvalue,SpaceStep*SpaceStep/(4*eta));
    else
        TimeStep = CFL*SpaceStep/Eigenvalue;

    // Recompute IterationNb

    if (RemainingTime <= 0.)
    {
        IterationNb = IterationNo;
    }
    else if (RemainingTime < TimeStep)
    {
        TimeStep = RemainingTime;
        IterationNb = IterationNo + 1;
    }
    else
    {
        RemainingIterations =  (int)(RemainingTime/TimeStep);
        IterationNb = IterationNo + RemainingIterations;
    }

    return;

}

Here is the caller graph for this function:

Vector ArtificialViscosity	(	const Vector &	Cell1,
		const Vector &	Cell2,
		real	dx,
		int	AxisNo
	)

Returns the artificial diffusion source terms in the cell UserCell.

Parameters

Cell1	Left cell value
Cell2	Right cell value
AxisNo	Axis of interest

Returns

Vector

X - direction

Y - direction

Z - direction

{
    // ---  Local variables ---
    Vector  ML(3), MR(3);
    Vector  Result(QuantityNb);
    real    EL,ER,RL,RR;
    real    viscR, viscX, viscY, viscZ, viscE;


    for(int i=1; i <= 3; i++){
        ML.setValue(i, Cell1.value(i+1));
        MR.setValue(i, Cell2.value(i+1));
    }

    RL = Cell1.value(1);
    RR = Cell2.value(1);
    EL = Cell1.value(5);
    ER = Cell2.value(5);


    if(AxisNo == 1){
        viscR = (RR - RL)/dx;
        viscE = (ER - EL)/dx;

        viscX = (MR.value(1) - ML.value(1))/dx;
        viscY = (MR.value(2) - ML.value(2))/dx;
        viscZ = (MR.value(3) - ML.value(3))/dx;


    }else if(AxisNo == 2){
        viscR = (RR - RL)/dx;
        viscE = (ER - EL)/dx;

        viscX = (MR.value(1) - ML.value(1))/dx;
        viscY = (MR.value(2) - ML.value(2))/dx;
        viscZ = (MR.value(3) - ML.value(3))/dx;

    }else{
        viscR = (RR - RL)/dx;
        viscE = (ER - EL)/dx;

        viscX = (MR.value(1) - ML.value(1))/dx;
        viscY = (MR.value(2) - ML.value(2))/dx;
        viscZ = (MR.value(3) - ML.value(3))/dx;

    }


    Result.setZero();

    // These values will be added to the numerical flux
    Result.setValue(1, chi*viscR);
    Result.setValue(2, chi*viscX);
    Result.setValue(3, chi*viscY);
    Result.setValue(4, chi*viscZ);
    Result.setValue(5, chi*viscE);

    return Result;

}

Here is the caller graph for this function:

void Backup

(

Node *

Root

)

Stores the tree structure and data in order to restart a multiresolution computation.

• Root denotes the pointer to the first node of the tree structure.

  Parameters

  Root

  Root

  Returns

  void

{
    Root->backup();
}

Here is the caller graph for this function:

void Backup

(

FineMesh *

Root

)

Stores the data countained in a regular mesh Root in order to restart a finite volume computation.

Parameters

Root

Root

Returns

void

{
    Root->backup();
}

int BC	(	int	i,
		int	AxisNo,
		int	N = (1<< ScaleNb)
	)

Returns the position of i taking into account the boundary conditions in the direction AxisNo. The number of points in this direction is N.
Example: for AxisNo=1 and for N=256, i must be between 0 and 255. If i=-1, the function returns 255 for periodic boundary conditions and 0 for Neumann boundary conditions.

Parameters

i	Position
AxisNo	Axis of interest
N	Defaults to (1<<ScaleNb).

Returns

int

{
  int result=-999999;

    if (AxisNo > Dimension)
        return 0;

#if defined PARMPI
    if (CMin[AxisNo] == 3) result=i;    //Periodic
    else
    {
        if (i<0) if ((coords[0]==0 && AxisNo==1) || (coords[1]==0 && AxisNo==2) || (coords[2]==0 && AxisNo==3)) result=-i-1;  //Neumann
        if (i>=N)       if ((coords[0]==CartDims[0]-1 && AxisNo==1) ||  (coords[1]==CartDims[1]-1 && AxisNo==2)
                                         || (coords[2]==CartDims[2]-1 && AxisNo==3)) result=2*N-i-1;
      if (result==-999999) result=i;    //Not the boundary, simple cell from anouther CPU
    }


#else
    if (CMin[AxisNo] == 3)
        result = (i+N)%N;                        // Periodic
    else if(CMin[AxisNo] == 2)
      result = ((i+N)/N==1)? i : (2*N-i-1)%N;  // Neumann

#endif

    return result;
}

Here is the caller graph for this function:

int BoundaryRegion

(

const Vector &

X

)

Returns the boundary region type at the position X=(x,y,z).
The returned value correspond to: 0 = Fluid (not in the boundary) 1 = Inflow 2 = Outflow 3 = Solid with free-slip walls 4 = Solid with adiabatic walls 5 = Solid with isothermal walls.

Parameters

X

Vector

Returns

int

{
    real x=0., y=0., z=0.;

    int Fluid;
    int Inflow;
    int Outflow;
    int FreeSlipSolid;
    int AdiabaticSolid;
    int IsothermalSolid;

    int Region;

    // Only in UseBoundaryRegions = true

    if (!UseBoundaryRegions) return 0;

    // --- Init values --

    Fluid       = 0;
    Inflow      = 1;
    Outflow     = 2;
    FreeSlipSolid   = 3;
    AdiabaticSolid  = 4;
    IsothermalSolid = 5;

    Region = Fluid;

    x = X.value(1);
    y = (Dimension > 1)? X.value(2):0.;
    z = (Dimension > 2)? X.value(3):0.;

    #include "carmen.bc"

    return Region;
}

Here is the caller graph for this function:

real ComputedTolerance

(

const int

ScaleNo

)

Returns the computed tolerance at the scale ScaleNo, either using Harten or Donoho thresholding (if CVS=true).

Parameters

ScaleNo

Level of interest.

Returns

double

{
//if ThresholdNorm==0 const Tolerance, else L1 Harten norm

    if(ThresholdNorm)
        return((Tolerance/GlobalVolume)*exp(Dimension*(ScaleNo-ScaleNb+1)*log(2.)));
    else
        return(Tolerance);

}

void CPUExchange	(	FineMesh *	Root,
		int
	)

Parallel function DOES NOT WORK!

Parameters

Root

Fine mesh

Returns

void

                                        {
  CommTimer.start();
#if defined PARMPI
  int i,k;
    int exNb=0;

  WhatSend=WS;
  CellElementsNb=0;

  for (i=0;i<16;i++) {
    k=1<<i;
    if ((WS & k) != 0) CellElementsNb++;
  }

  static bool ft=true;
//  if (ft==true) {
    CreateMPIType(Root);
//    CreateMPILinks();
//    ft=false;
//  }

//  MPI_Startall(4*Dimension,req);


//Send
    switch (MPISendType) {
    case 0:
      MPI_Ibsend(MPI_BOTTOM, 1, MPItypeSiL, rank_il, 100, comm_cart ,&req[exNb++]);
      MPI_Ibsend(MPI_BOTTOM, 1, MPItypeSiU, rank_iu, 200, comm_cart, &req[exNb++] );
        break;

    case 10:
      MPI_Isend(MPI_BOTTOM, 1, MPItypeSiL, rank_il, 100, comm_cart,&req[exNb++]);
      MPI_Isend(MPI_BOTTOM, 1, MPItypeSiU, rank_iu, 200, comm_cart,&req[exNb++]);
        break;

    case 20:
      MPI_Issend(MPI_BOTTOM, 1, MPItypeSiL, rank_il, 100, comm_cart,&req[exNb++]);
      MPI_Issend(MPI_BOTTOM, 1, MPItypeSiU, rank_iu, 200, comm_cart,&req[exNb++]);
        break;
    }

    if (Dimension >= 2) {
        switch (MPISendType) {
        case 0:
          MPI_Ibsend(MPI_BOTTOM, 1, MPItypeSjL, rank_jl, 300, comm_cart,&req[exNb++]);
        MPI_Ibsend(MPI_BOTTOM, 1, MPItypeSjU, rank_ju, 400, comm_cart,&req[exNb++]);
          break;

        case 10:
          MPI_Isend(MPI_BOTTOM, 1, MPItypeSjL, rank_jl, 300, comm_cart,&req[exNb++]);
          MPI_Isend(MPI_BOTTOM, 1, MPItypeSjU, rank_ju, 400, comm_cart,&req[exNb++]);
            break;

        case 20:
          MPI_Issend(MPI_BOTTOM, 1, MPItypeSjL, rank_jl, 300, comm_cart,&req[exNb++]);
          MPI_Issend(MPI_BOTTOM, 1, MPItypeSjU, rank_ju, 400, comm_cart,&req[exNb++]);
            break;
        }
    }

    if (Dimension == 3) {
        switch (MPISendType) {
        case 0:
          MPI_Ibsend(MPI_BOTTOM, 1, MPItypeSkL, rank_kl, 500, comm_cart,&req[exNb++]);
          MPI_Ibsend(MPI_BOTTOM, 1, MPItypeSkU, rank_ku, 600, comm_cart,&req[exNb++]);
            break;

        case 10:
          MPI_Isend(MPI_BOTTOM, 1, MPItypeSkL, rank_kl, 500, comm_cart,&req[exNb++]);
          MPI_Isend(MPI_BOTTOM, 1, MPItypeSkU, rank_ku, 600, comm_cart,&req[exNb++]);
            break;

        case 20:
        MPI_Issend(MPI_BOTTOM, 1, MPItypeSkL, rank_kl, 500, comm_cart,&req[exNb++]);
          MPI_Issend(MPI_BOTTOM, 1, MPItypeSkU, rank_ku, 600, comm_cart,&req[exNb++]);
            break;
        }
    }

//Recv

    if (MPIRecvType==0) {
        MPI_Recv(MPI_BOTTOM, 1, MPItypeRiL, rank_il, 200, comm_cart, &st[6]);
        MPI_Recv(MPI_BOTTOM, 1, MPItypeRiU, rank_iu, 100, comm_cart, &st[7]);
    } else
    {
        MPI_Irecv(MPI_BOTTOM, 1, MPItypeRiL, rank_il, 200, comm_cart, &req[exNb++]);
        MPI_Irecv(MPI_BOTTOM, 1, MPItypeRiU, rank_iu, 100, comm_cart, &req[exNb++]);
    }

    if (Dimension >= 2) {
        if (MPIRecvType==0) {
            MPI_Recv(MPI_BOTTOM, 1, MPItypeRjL, rank_jl, 400, comm_cart, &st[8]);
            MPI_Recv(MPI_BOTTOM, 1, MPItypeRjU, rank_ju, 300, comm_cart, &st[9]);
        } else
        {
            MPI_Irecv(MPI_BOTTOM, 1, MPItypeRjL, rank_jl, 400, comm_cart, &req[exNb++]);
            MPI_Irecv(MPI_BOTTOM, 1, MPItypeRjU, rank_ju, 300, comm_cart, &req[exNb++]);
        }
    }

    if (Dimension == 3) {
        if (MPIRecvType==0) {
            MPI_Recv(MPI_BOTTOM, 1, MPItypeRkL, rank_kl, 600, comm_cart, &st[10]);
            MPI_Recv(MPI_BOTTOM, 1, MPItypeRkU, rank_ku, 500, comm_cart, &st[11]);
        } else
        {
            MPI_Irecv(MPI_BOTTOM, 1, MPItypeRkL, rank_kl, 600, comm_cart, &req[exNb++]);
            MPI_Irecv(MPI_BOTTOM, 1, MPItypeRkU, rank_ku, 500, comm_cart, &req[exNb++]);
        }
    }

  FreeMPIType();
#endif
  CommTimer.stop();
}

Here is the caller graph for this function:

double CPUTimeRef	(	int	iterations,
		int	scales
	)

Returns the time required by a finite volume computation using iterations iterations and scales scales. It is use to estimate the CPU time compression.

Parameters

iterations	Number of iterations.
scales	Scales

Returns

double

{
    // --- Local variables ---------------------------------------------------

    int OldIterationNb=0;
    int OldScaleNb=0;
    real OldTimeStep=0.;
    bool ConstantTimeStepOld=ConstantTimeStep;

    double result=0.;

    Timer CPURef;
    FineMesh* MeshRef;

    // --- Execution ---------------------------------------------------------

    // Toggle on : Compute reference CPU time

    ComputeCPUTimeRef = true;

    // Toggle off : Constant time step

    ConstantTimeStep = true;

    // backup values of IterationNb and ScaleNb

    OldIterationNb  = IterationNb;
    OldScaleNb  = ScaleNb;
    OldTimeStep     = TimeStep;

    // use reference values
    IterationNb     = iterations;
    ScaleNb         = scales;
    TimeStep    = 0.;

    one_D=1; two_D=1;
    if (Dimension >= 2) one_D=1<<ScaleNb;
    if (Dimension == 3) two_D=1<<ScaleNb;

    // init mesh
    MeshRef = new FineMesh;

    // Iterate on time

    for (IterationNo = 1; IterationNo <= IterationNb; IterationNo++)
    {
        // start timer
        CPURef.start();

        // Compute time evolution
        TimeEvolution(MeshRef);

        // check CPU Time
        CPURef.check();

        // stop timer
        CPURef.stop();
    }

    // Compute CPUTimeRef
    result = CPURef.CPUTime();
    result *= 1./IterationNb;
    result *= 1<<(Dimension*(OldScaleNb-ScaleNb));

    // delete MeshRef
    delete MeshRef;

    // restore values of IterationNb and ScaleNb
    IterationNb = OldIterationNb;
    ScaleNb = OldScaleNb;
    TimeStep = OldTimeStep;
    IterationNo = 0;

    one_D=1; two_D=1;
    if (Dimension >= 2) one_D=1<<ScaleNb;
    if (Dimension == 3) two_D=1<<ScaleNb;

    // Toggle off : Compute reference CPU time

    ComputeCPUTimeRef = false;

    // Restore the value of ConstantTimeStep

    ConstantTimeStep = ConstantTimeStepOld;

    return result;
}

Here is the caller graph for this function:

void CreateMPILinks

(

)

Parallel function DOES NOT WORK!

Returns

void

                      {
    int exNb;
    exNb=0;
#if defined PARMPI

//Send

    switch (MPISendType) {
    case 0:
      MPI_Bsend_init(MPI_BOTTOM, 1, MPItypeSiL, rank_il, 100, comm_cart ,&req[exNb++]);
      MPI_Bsend_init(MPI_BOTTOM, 1, MPItypeSiU, rank_iu, 200, comm_cart, &req[exNb++]);
        break;

    case 10:
      MPI_Send_init(MPI_BOTTOM, 1, MPItypeSiL, rank_il, 100, comm_cart,&req[exNb++]);
      MPI_Send_init(MPI_BOTTOM, 1, MPItypeSiU, rank_iu, 200, comm_cart,&req[exNb++]);
        break;

    case 20:
      MPI_Ssend_init(MPI_BOTTOM, 1, MPItypeSiL, rank_il, 100, comm_cart,&req[exNb++]);
      MPI_Ssend_init(MPI_BOTTOM, 1, MPItypeSiU, rank_iu, 200, comm_cart,&req[exNb++]);
        break;
    }

    if (Dimension >= 2) {
        switch (MPISendType) {
        case 0:
          MPI_Bsend_init(MPI_BOTTOM, 1, MPItypeSjL, rank_jl, 300, comm_cart,&req[exNb++]);
        MPI_Bsend_init(MPI_BOTTOM, 1, MPItypeSjU, rank_ju, 400, comm_cart,&req[exNb++]);
          break;

        case 10:
          MPI_Send_init(MPI_BOTTOM, 1, MPItypeSjL, rank_jl, 300, comm_cart,&req[exNb++]);
          MPI_Send_init(MPI_BOTTOM, 1, MPItypeSjU, rank_ju, 400, comm_cart,&req[exNb++]);
            break;

        case 20:
          MPI_Ssend_init(MPI_BOTTOM, 1, MPItypeSjL, rank_jl, 300, comm_cart,&req[exNb++]);
          MPI_Ssend_init(MPI_BOTTOM, 1, MPItypeSjU, rank_ju, 400, comm_cart,&req[exNb++]);
            break;
        }
    }

    if (Dimension == 3) {
        switch (MPISendType) {
        case 0:
          MPI_Bsend_init(MPI_BOTTOM, 1, MPItypeSkL, rank_kl, 500, comm_cart,&req[exNb++]);
          MPI_Bsend_init(MPI_BOTTOM, 1, MPItypeSkU, rank_ku, 600, comm_cart,&req[exNb++]);
            break;

        case 10:
          MPI_Send_init(MPI_BOTTOM, 1, MPItypeSkL, rank_kl, 500, comm_cart,&req[exNb++]);
          MPI_Send_init(MPI_BOTTOM, 1, MPItypeSkU, rank_ku, 600, comm_cart,&req[exNb++]);
            break;

        case 20:
        MPI_Ssend_init(MPI_BOTTOM, 1, MPItypeSkL, rank_kl, 500, comm_cart,&req[exNb++]);
          MPI_Ssend_init(MPI_BOTTOM, 1, MPItypeSkU, rank_ku, 600, comm_cart,&req[exNb++]);
            break;
        }
    }

//Recv

    MPI_Recv_init(MPI_BOTTOM, 1, MPItypeRiL, rank_il, 200, comm_cart, &req[exNb++]);
    MPI_Recv_init(MPI_BOTTOM, 1, MPItypeRiU, rank_iu, 100, comm_cart, &req[exNb++]);

    if (Dimension >= 2) {
        MPI_Recv_init(MPI_BOTTOM, 1, MPItypeRjL, rank_jl, 400, comm_cart, &req[exNb++]);
        MPI_Recv_init(MPI_BOTTOM, 1, MPItypeRjU, rank_ju, 300, comm_cart, &req[exNb++]);
    }

    if (Dimension == 3) {
        MPI_Recv_init(MPI_BOTTOM, 1, MPItypeRkL, rank_kl, 600, comm_cart, &req[exNb++]);
        MPI_Recv_init(MPI_BOTTOM, 1, MPItypeRkU, rank_ku, 500, comm_cart, &req[exNb++]);
    }
#endif
}

void CreateMPITopology

(

)

Parallel function DOES NOT WORK!

Returns

void

                         {
#if defined PARMPI
 int src;
 int periods[]={1,1,1};
 CartDims[0]=CartDims[1]=CartDims[2]=0;

 MPI_Dims_create(size,Dimension,CartDims);
 MPI_Cart_create(MPI_COMM_WORLD,Dimension,CartDims,periods,1,&comm_cart);
 MPI_Comm_rank(comm_cart, &rank);
 MPI_Cart_coords(comm_cart,rank,Dimension,coords);

 MPI_Cart_shift(comm_cart, 0, -1, &src, &rank_il);
 MPI_Cart_shift(comm_cart, 0, 1, &src, &rank_iu);

 if (Dimension >= 2) {
   MPI_Cart_shift(comm_cart, 1, -1, &src, &rank_jl);
   MPI_Cart_shift(comm_cart, 1, 1, &src, &rank_ju);
 }

 if (Dimension == 3) {
  MPI_Cart_shift(comm_cart, 2, -1, &src, &rank_kl);
  MPI_Cart_shift(comm_cart, 2, 1, &src, &rank_ku);
 }
#endif
}

Here is the caller graph for this function:

void CreateMPIType

(

FineMesh *

Root

)

                                   {
#if defined PARMPI
  int i,j,k;
    int n,d,l;

    Cell *MeshCell;
    MeshCell=Root->MeshCell;

    n=0;
  for (l=0;l<NeighbourNb;l++)
        for (j=0;j<one_D;j++)
            for (k=0;k<two_D;k++)   FillNbAddr(Root->Neighbour_iL,l,j,k,n);
  MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeRiL);
    MPI_Type_commit(&MPItypeRiL);

    n=0;
  for (l=0;l<NeighbourNb;l++)
        for (j=0;j<one_D;j++)
            for (k=0;k<two_D;k++)   FillNbAddr(Root->Neighbour_iU,l,j,k,n);
  MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeRiU);
    MPI_Type_commit(&MPItypeRiU);

    n=0;
    for (l=0;l<NeighbourNb;l++)
        for (j=0;j<one_D;j++)
            for (k=0;k<two_D;k++)   {
                i=l;
                d=i + (1<<ScaleNb)*(j + (1<<ScaleNb)*k);
                FillCellAddr(MeshCell,d,n);
            }
  MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeSiL);
    MPI_Type_commit(&MPItypeSiL);

    n=0;
    for (l=0;l<NeighbourNb;l++)
        for (j=0;j<one_D;j++)
            for (k=0;k<two_D;k++)   {
                i=(1<<ScaleNb)-NeighbourNb+l;
                d=i + (1<<ScaleNb)*(j + (1<<ScaleNb)*k);
                FillCellAddr(MeshCell,d,n);
            }
  MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeSiU);
    MPI_Type_commit(&MPItypeSiU);

    if (Dimension >= 2) {
        n=0;
      for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (k=0;k<two_D;k++)   FillNbAddr(Root->Neighbour_jL,l,i,k,n);
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeRjL);
        MPI_Type_commit(&MPItypeRjL);

        n=0;
      for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (k=0;k<two_D;k++)   FillNbAddr(Root->Neighbour_jU,l,i,k,n);
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeRjU);
        MPI_Type_commit(&MPItypeRjU);

        n=0;
        for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (k=0;k<two_D;k++)   {
                    j=l;
                    d=i + (1<<ScaleNb)*(j + (1<<ScaleNb)*k);
                    FillCellAddr(MeshCell,d,n);
                }
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeSjL);
        MPI_Type_commit(&MPItypeSjL);

        n=0;
        for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (k=0;k<two_D;k++)   {
                    j=(1<<ScaleNb)-NeighbourNb+l;
                    d=i + (1<<ScaleNb)*(j + (1<<ScaleNb)*k);
                    FillCellAddr(MeshCell,d,n);
                }
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeSjU);
        MPI_Type_commit(&MPItypeSjU);
  }


    if (Dimension == 3) {
        n=0;
      for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (j=0;j<two_D;j++)   FillNbAddr(Root->Neighbour_kL,l,i,j,n);
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeRkL);
        MPI_Type_commit(&MPItypeRkL);

        n=0;
      for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (j=0;j<two_D;j++)   FillNbAddr(Root->Neighbour_kU,l,i,j,n);
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeRkU);
        MPI_Type_commit(&MPItypeRkU);

        n=0;
        for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (j=0;j<two_D;j++)   {
                    k=l;
                    d=i + (1<<ScaleNb)*(j + (1<<ScaleNb)*k);
                    FillCellAddr(MeshCell,d,n);
                }
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeSkL);
        MPI_Type_commit(&MPItypeSkL);

        n=0;
        for (l=0;l<NeighbourNb;l++)
            for (i=0;i<one_D;i++)
                for (j=0;j<two_D;j++)   {
                    k=(1<<ScaleNb)-NeighbourNb+l;
                    d=i + (1<<ScaleNb)*(j + (1<<ScaleNb)*k);
                    FillCellAddr(MeshCell,d,n);
                }
      MPI_Type_hindexed(CellElementsNb*NeighbourNb*one_D*two_D,blocklen,disp,MPI_Type,&MPItypeSkU);
        MPI_Type_commit(&MPItypeSkU);
    }

#endif
}

Here is the caller graph for this function:

int DigitNumber

(

int

arg

)

Returns the number of digits of the integer arg.

Parameters

arg

Argument

Returns

int

{
    int result;
    int i;

    result = 0;
    i = arg;

    while(i != 0)
    {
        i/=10;
        result++;
    }

    return result;
}

Here is the caller graph for this function:

int FileWrite	(	FILE *	f,
		const char *	format,
		real	arg
	)

Writes in binary or ASCII mode the real number arg into the file f with the format format. The global parameter DataIsBinary determines this choice.

Parameters

f	File name
format	Format
arg	Argument

Returns

int

{
    int result;
    real x;

    x = arg;

    if (DataIsBinary)
        result = fwrite(&x,sizeof(real),1,f);
    else
        result = fprintf(f, format, x);


  return result;
}

Here is the caller graph for this function:

Vector Flux	(	Cell &	Cell1,
		Cell &	Cell2,
		Cell &	Cell3,
		Cell &	Cell4,
		int	AxisNo
	)

Returns the flux at the interface between Cell2 and Cell3. Here a 4-point space scheme is used. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides.

Parameters

Cell1	second neighbour on the left side
Cell2	first neighbour on the left side
Cell3	first neighbour on the right side
Cell4	second neighbour on the right side
AxisNo	Axis of interest

Returns

Vector

{
    // --- Local variables ---

    Vector Result(QuantityNb);

    Cell C1, C2, C3, C4;

    int BoundaryCell1 = BoundaryRegion(Cell1.center());
    int BoundaryCell2 = BoundaryRegion(Cell2.center());
    int BoundaryCell3 = BoundaryRegion(Cell3.center());
    int BoundaryCell4 = BoundaryRegion(Cell4.center());

    bool UseBoundaryCells = (UseBoundaryRegions && (BoundaryCell1!=0 || BoundaryCell2!=0 || BoundaryCell3!=0 || BoundaryCell4!=0));
    
    // --- Take into account boundary conditions ---

            if (UseBoundaryCells)
                GetBoundaryCells(Cell1, Cell2, Cell3, Cell4, C1, C2, C3, C4, AxisNo);

    switch(SchemeNb)
    {
        case 1:
        default:
            if (UseBoundaryCells)
                Result = SchemeHLL(C1, C2, C3, C4, AxisNo);
            else
                Result = SchemeHLL(Cell1, Cell2, Cell3, Cell4, AxisNo);
            break;

        case 2:
            if (UseBoundaryCells)
                Result = SchemeHLLD(C1, C2, C3, C4, AxisNo);
            else
                Result = SchemeHLLD(Cell1, Cell2, Cell3, Cell4, AxisNo);
            break;
        break;

    }


    return Result;
}

Here is the caller graph for this function:

void fluxCorrection	(	Vector &	Flux,
		const Vector &	AvgL,
		const Vector &	AvgR,
		int	AxisNo
	)

This function apply the divergence-free correction to the numerical flux.

Parameters

Flux	Numerical flux vector
AvgL	Left average vector
AvgR	Right average vector
AxisNo	Axis of interest

Returns

void

{
    auxvar = Flux.value(AxisNo+6);

    Flux.setValue(AxisNo+6, Flux.value(AxisNo+6) + (AvgL.value(6) +    .5*(AvgR.value(6) - AvgL.value(6))
                                          - ch*.5*(AvgR.value(AxisNo+6) - AvgL.value(AxisNo+6))));

    Flux.setValue(6, ch*ch*(AvgL.value(AxisNo+6) + .5*(AvgR.value(AxisNo+6) - AvgL.value(AxisNo+6))
                                          - .5*(AvgR.value(6) - AvgL.value(6))/ch));

}

Here is the caller graph for this function:

Vector FluxX

(

const Vector &

Avg

)

Returns the physical flux of MHD equations in X direction.

Parameters

Avg	Average vector.

Returns

Vector

{
    real rho;
    real vx, vy, vz;
    real pre, e;
    real Bx, By, Bz;
    real Bx2, By2, Bz2, B2;
    real vx2, vy2, vz2, v2;
    real half = 0.5;
    Vector F(QuantityNb);

    //Variables
    rho = Avg.value(1);
    vx  = Avg.value(2)/rho;
    vy  = Avg.value(3)/rho;
    vz  = Avg.value(4)/rho;
    e   = Avg.value(5);
    Bx  = Avg.value(7);
    By  = Avg.value(8);
    Bz  = Avg.value(9);

    Bx2 = Bx*Bx;
    By2 = By*By;
    Bz2 = Bz*Bz;
    B2  = half*(Bz2+Bx2+By2);

    vx2 = vx*vx;
    vy2 = vy*vy;
    vz2 = vz*vz;
    v2  = half*(vz2+vx2+vy2);

    //pressure
    pre = (Gamma -1.)*(e - rho*v2 - B2);

    //Physical flux - x-direction
    F.setValue(1,rho*vx);
    F.setValue(2,rho*vx2 + pre + half*(Bz2+By2-Bx2));
    F.setValue(3,rho*vx*vy - Bx*By);
    F.setValue(4,rho*vx*vz - Bx*Bz);
    F.setValue(5,(e + pre + B2)*vx - Bx*(vx*Bx + vy*By + vz*Bz) );
    F.setValue(6,0.0);
    F.setValue(7,0.0);
    F.setValue(8,vx*By - vy*Bx);
    F.setValue(9,vx*Bz - vz*Bx);


    return F;
}

Here is the caller graph for this function:

Vector FluxY

(

const Vector &

Avg

)

Returns the physical flux of MHD equations in Y direction.

Parameters

Avg	Average vector.

Returns

Vector

{
    real rho;
    real vx, vy, vz;
    real pre, e;
    real Bx, By, Bz;
    real Bx2, By2, Bz2, B2;
    real vx2, vy2, vz2, v2;
    real half = 0.5;

    Vector G(QuantityNb);

    //Variables
    rho = Avg.value(1);
    vx  = Avg.value(2)/rho;
    vy  = Avg.value(3)/rho;
    vz  = Avg.value(4)/rho;
    e   = Avg.value(5);
    Bx  = Avg.value(7);
    By  = Avg.value(8);
    Bz  = Avg.value(9);

    Bx2 = Bx*Bx;
    By2 = By*By;
    Bz2 = Bz*Bz;
    B2  = half*(Bz2+Bx2+By2);

    vx2 = vx*vx;
    vy2 = vy*vy;
    vz2 = vz*vz;
    v2  = half*(vz2+vx2+vy2);

    //pressure
    pre = (Gamma -1.)*(e - rho*v2 - B2);

    //Physical flux - y-direction
    G.setValue(1,rho*vy);
    G.setValue(2,rho*vx*vy - Bx*By);
    G.setValue(3,rho*vy2 + pre + half*(Bx2+Bz2-By2));
    G.setValue(4,rho*vy*vz - By*Bz);
    G.setValue(5,(e + pre + B2)*vy - By*(vx*Bx + vy*By + vz*Bz));
    G.setValue(6,0.0);
    G.setValue(7,vy*Bx - vx*By);
    G.setValue(8,0.0);
    G.setValue(9,vy*Bz - vz*By);
    return G;
}

Here is the caller graph for this function:

Vector FluxZ

(

const Vector &

Avg

)

Returns the physical flux of MHD equations in Z direction.

Parameters

Avg	Average vector.

Returns

Vector

{
    real rho;
    real vx, vy, vz;
    real pre, e;
    real Bx, By, Bz;
    real Bx2, By2, Bz2, B2;
    real vx2, vy2, vz2, v2;
    real half = 0.5;

    Vector H(QuantityNb);

    //Variables
    rho = Avg.value(1);
    vx  = Avg.value(2)/rho;
    vy  = Avg.value(3)/rho;
    vz  = Avg.value(4)/rho;
    e   = Avg.value(5);
    Bx  = Avg.value(7);
    By  = Avg.value(8);
    Bz  = Avg.value(9);

    Bx2 = Bx*Bx;
    By2 = By*By;
    Bz2 = Bz*Bz;
    B2  = half*(Bz2+Bx2+By2);

    vx2 = vx*vx;
    vy2 = vy*vy;
    vz2 = vz*vz;
    v2  = half*(vz2+vx2+vy2);

    //pressure
    pre = (Gamma -1.)*(e - rho*v2 - B2);

    //Physical flux - y-direction
    H.setValue(1,rho*vz);
    H.setValue(2,rho*vz*vx - Bz*Bx);
    H.setValue(3,rho*vz*vy - Bz*By);
    H.setValue(4,rho*vz2 + pre + half*(Bx2+By2-Bz2));
    H.setValue(5,(e + pre + B2)*vz - Bz*(vx*Bx + vy*By + vz*Bz));
    H.setValue(6,0.0);
    H.setValue(7,vz*Bx - vx*Bz);
    H.setValue(8,vz*By - vy*Bz);
    H.setValue(9,0.0);

    return H;
}

Here is the caller graph for this function:

void FreeMPIType

(

)

Parallel function DOES NOT WORK!

Returns

void

                   {
#if defined PARMPI
    MPI_Type_free(&MPItypeSiL);
    MPI_Type_free(&MPItypeSiU);
    MPI_Type_free(&MPItypeRiL);
    MPI_Type_free(&MPItypeRiU);

    if (Dimension >= 2) {
        MPI_Type_free(&MPItypeSjL);
        MPI_Type_free(&MPItypeSjU);
        MPI_Type_free(&MPItypeRjL);
        MPI_Type_free(&MPItypeRjU);
    }

    if (Dimension == 3) {
        MPI_Type_free(&MPItypeSkL);
        MPI_Type_free(&MPItypeSkU);
        MPI_Type_free(&MPItypeRkL);
        MPI_Type_free(&MPItypeRkU);
    }
#endif
}

Here is the caller graph for this function:

void GetBoundaryCells	(	Cell &	Cell1,
		Cell &	Cell2,
		Cell &	Cell3,
		Cell &	Cell4,
		Cell &	C1,
		Cell &	C2,
		Cell &	C3,
		Cell &	C4,
		int	AxisNo
	)

Transform the 4 cells of the flux Cell1, Cell2, Cell3, Cell4 into C1, C2, C3, C4, to take into account boundary conditions (used in Flux.cpp).

Parameters

Cell1	second neighbour on the left side
Cell2	first neighbour on the left side
Cell3	first neighbour on the right side
Cell4	second neighbour on the right side
C1	Auxiliar cell1
C2	Auxiliar cell2
C3	Auxiliar cell3
C4	Auxiliar cell4
AxisNo	...

Returns

void

{
  // --- Local variables ---

    int InCell1, InCell2, InCell3, InCell4;                     // Boundary conditions in cells 1, 2, 3, 4
    real P1, P2, P3, P4;                            // Pressures in cells 1, 2, 3, 4
    real T1, T2, T3, T4;                            // Temperatures in cells 1, 2, 3, 4
    real rho1, rho2, rho3, rho4;                        // Densities in cells 1, 2, 3, 4
    Vector V1(Dimension), V2(Dimension), V3(Dimension), V4(Dimension);  // Velocities in cells 1, 2, 3, 4
    real e1, e2, e3, e4;                            // Energies in cell 1, 2, 3, 4
    real Y1=0., Y2=0., Y3=0., Y4=0.;                    // Partial mass in cell 1, 2, 3, 4

    int i;          // Counter

    // --- Init C1, C2, C3, C4 ---

    C1 = Cell1;
    C2 = Cell2;
    C3 = Cell3;
    C4 = Cell4;

    // --- Depending on the boundary region type, transform C1, C2, C3, C4 ---

    InCell1 = BoundaryRegion(Cell1.center());
    InCell2 = BoundaryRegion(Cell2.center());
    InCell3 = BoundaryRegion(Cell3.center());
    InCell4 = BoundaryRegion(Cell4.center());

    // --- Cell2 IN THE BOUNDARY, Cell3 IN THE FLUID -----------------------------------------

    if (InCell2 !=0 && InCell3 == 0)
    {
        switch(InCell2)
        {
            // INFLOW
            case 1:
                // Dirichlet on temperature
                T2 = Cell2.temperature();
                T1 = Cell1.temperature();

                // Extrapolate pressure
                P2 = Cell3.oldPressure();
                P1 = P2;
                // P1 = 2*P2 - Cell3.pressure();

                // Compute density
                rho2 = Gamma*Ma*Ma*P2/T2;
                rho1 = Gamma*Ma*Ma*P1/T1;

                // Dirichlet on momentum
                V2 = (Cell2.density()/rho2)*Cell2.velocity();
                V1 = (Cell1.density()/rho1)*Cell1.velocity();

                // Dirichlet on partial mass
                if (ScalarEqNb == 1)
                {
                    Y2 = Cell2.average(Dimension+3)/Cell2.density();
                    Y1 = Cell1.average(Dimension+3)/Cell1.density();
                }

                // Compute energies
                e2 = P2/((Gamma-1.)*rho2) + 0.5*N2(V2);
                e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);

                // Correct C1 and C2
                C2.setAverage(1, rho2);
                C1.setAverage(1, rho1);
                for (i=1; i<=Dimension; i++)
                {
                    C2.setAverage(i+1, rho2*V2.value(i));
                    C1.setAverage(i+1, rho1*V1.value(i));
                }
                C2.setAverage(Dimension+2,rho2*e2);
                C1.setAverage(Dimension+2,rho1*e1);

                if (ScalarEqNb == 1)
                {
                    C2.setAverage(Dimension+3, rho2*Y2);
                    C1.setAverage(Dimension+3, rho1*Y1);
                }
                break;

            // OUTFLOW : use the old value of the neighbour
            case 2:
                C2.setAverage(Cell3.oldAverage());
                C1.setAverage(Cell3.oldAverage());

                // Also change the values in the boundary
                Cell2.setAverage(C2.average());
                Cell1.setAverage(C1.average());
                break;

            // FREE-SLIP WALL : Neuman on all quantities
            case 3:

                C2 = Cell3;
                C1 = Cell4;
                break;

            // ADIABATIC WALL
            case 4:

                // Dirichlet on velocity
                V2 = Cell2.velocity();
                V1 = Cell1.velocity();

                // Neuman on temperature
                T2 = Cell3.temperature();
                T1 = Cell4.temperature();

                // Neuman on pressure
                P2 = Cell3.pressure();
                P1 = Cell4.pressure();

                // Extrapolate pressure
                //P2 = 2*Cell3.pressure()-Cell4.pressure();
                //P1 = P2;

                // Compute densities
                rho2 = Gamma*Ma*Ma*P2/T2;
                rho1 = Gamma*Ma*Ma*P1/T1;

                // Compute energies
                e2 = P2/((Gamma-1.)*rho2) + 0.5*N2(V2);
                e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);

                // Correct C1 and C2
                C2.setAverage(1, rho2);
                C1.setAverage(1, rho1);
                for (i=1; i<=Dimension; i++)
                {
                    C2.setAverage(i+1, rho2*V2.value(i));
                    C1.setAverage(i+1, rho1*V1.value(i));
                }
                C2.setAverage(Dimension+2,rho2*e2);
                C1.setAverage(Dimension+2,rho1*e1);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                {
                    C2.setAverage(Dimension+3, Cell3.average(Dimension+3));
                    C1.setAverage(Dimension+3, Cell4.average(Dimension+3));
                }

                break;

            // ISOTHERMAL WALL
            case 5:

                // Dirichlet on velocity
                V2 = Cell2.velocity();
                V1 = Cell1.velocity();

                // Dirichlet on temperature
                T2 = Cell2.temperature();
                T1 = Cell1.temperature();

                // Neuman on pressure
                P2 = Cell3.pressure();
                P1 = Cell4.pressure();

                // Extrapolate pressure
                //P2 = 2*Cell3.pressure()-Cell4.pressure();
                //P1 = P2;

                // Compute densities
                rho2 = Gamma*Ma*Ma*P2/T2;
                rho1 = Gamma*Ma*Ma*P1/T1;

                // Compute energies
                e2 = P2/((Gamma-1.)*rho2) + 0.5*N2(V2);
                e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);

                // Correct C1 and C2
                C2.setAverage(1, rho2);
                C1.setAverage(1, rho1);
                for (i=1; i<=Dimension; i++)
                {
                    C2.setAverage(i+1, rho2*V2.value(i));
                    C1.setAverage(i+1, rho1*V1.value(i));
                }
                C2.setAverage(Dimension+2,rho2*e2);
                C1.setAverage(Dimension+2,rho1*e1);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                {
                    C2.setAverage(Dimension+3, Cell3.average(Dimension+3));
                    C1.setAverage(Dimension+3, Cell4.average(Dimension+3));
                }

                break;
        };
    return;
    }

  // --- Cell1 IN THE BOUNDARY, Cell2 IN THE FLUID ------------------------------------------------

    if (InCell1 != 0 && InCell2 == 0)
    {
        switch(InCell1)
        {
            // INFLOW
            case 1:
                // Dirichlet on temperature
                T1 = Cell1.temperature();

                // Extrapolate pressure from old value
                P1 = Cell2.oldPressure();

                // Compute density
                rho1 = Gamma*Ma*Ma*P1/T1;

                // Dirichlet on momentum
                V1 = (Cell1.density()/rho1)*Cell1.velocity();

                // Dirichlet on partial mass
                if (ScalarEqNb == 1)
                    Y1 = Cell1.average(Dimension+3)/Cell1.density();

                // Compute energies
                e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);

                // Correct C1
                C1.setAverage(1, rho1);
                for (i=1; i<=Dimension; i++)
                    C1.setAverage(i+1, rho1*V1.value(i));
                C1.setAverage(Dimension+2,rho1*e1);

                if (ScalarEqNb == 1)
                    C1.setAverage(Dimension+3, rho1*Y1);
                break;

            // OUTFLOW : Get old value of the neighbour
            case 2:

                C1.setAverage(Cell2.oldAverage());
                break;

            // FREE-SLIP WALL : Neuman on all quantities
            case 3:

                C1 = Cell2;
                break;

            // ADIABATIC WALL
            case 4:

                // Dirichlet on velocity
                V1 = Cell1.velocity();

                // Neuman on temperature
                T1 = Cell2.temperature();

                // Neuman on pressure
                P1 = Cell2.pressure();

                // Extrapolate pressure
                //P1 = 2*Cell2.pressure()-Cell3.pressure();

                // Compute density
                rho1 = Gamma*Ma*Ma*P1/T1;

                // Compute energy
                e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);

                // Correct C1
                C1.setAverage(1, rho1);
                for (i=1; i<=Dimension; i++)
                    C1.setAverage(i+1, rho1*V1.value(i));
                C1.setAverage(Dimension+2,rho1*e1);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                    C1.setAverage(Dimension+3, Cell2.average(Dimension+3));

                break;

            // ISOTHERMAL WALL
            case 5:

                // Dirichlet on velocity
                V1 = Cell1.velocity();

                // Dirichlet on temperature
                T1 = Cell1.temperature();

                // Neuman on pressure
                P1 = Cell2.pressure();

                // Extrapolate pressure
                //P1 = 2*Cell2.pressure()-Cell3.pressure();

                // Compute density
                rho1 = Gamma*Ma*Ma*P1/T1;

                // Compute energies
                e1 = P1/((Gamma-1.)*rho1) + 0.5*N2(V1);

                // Correct C1
                C1.setAverage(1, rho1);
                for (i=1; i<=Dimension; i++)
                    C1.setAverage(i+1, rho1*V1.value(i));
                C1.setAverage(Dimension+2,rho1*e1);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                    C1.setAverage(Dimension+3, Cell2.average(Dimension+3));

                break;
        };
    return;
    }
    // --- Cell3 IN THE BOUNDARY, Cell2 IN THE FLUID -----------------------------------------

  if (InCell3 !=0 && InCell2 == 0)
    {
        switch(InCell3)
        {
            // INFLOW
            case 1:
                // Dirichlet on temperature
                T3 = Cell3.temperature();
                T4 = Cell4.temperature();

                // Extrapolate pressure from old value
                P3 = Cell2.oldPressure();
                P4 = P3;
                //P4 = 2*P3 - Cell2.pressure();

                // Compute densities
                rho3 = Gamma*Ma*Ma*P3/T3;
                rho4 = Gamma*Ma*Ma*P4/T4;

                // Dirichlet on momentum
                V3 = (Cell3.density()/rho3)*Cell3.velocity();
                V4 = (Cell4.density()/rho4)*Cell4.velocity();

                // Dirichlet on partial mass
                if (ScalarEqNb == 1)
                {
                    Y3 = Cell3.average(Dimension+3)/Cell3.density();
                    Y4 = Cell4.average(Dimension+3)/Cell4.density();
                }

                // Compute energies
                e3 = P3/((Gamma-1.)*rho3) + 0.5*N2(V3);
                e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);

                // Correct C1 and C2
                C3.setAverage(1, rho3);
                C4.setAverage(1, rho4);
                for (i=1; i<=Dimension; i++)
                {
                    C3.setAverage(i+1, rho3*V3.value(i));
                    C4.setAverage(i+1, rho4*V4.value(i));
                }
                C3.setAverage(Dimension+2,rho3*e3);
                C4.setAverage(Dimension+2,rho4*e4);

                if (ScalarEqNb == 1)
                {
                    C3.setAverage(Dimension+3,rho3*Y3);
                    C4.setAverage(Dimension+3,rho4*Y4);
                }
                break;

            // OUTFLOW
            case 2:

                C3.setAverage(Cell2.oldAverage());
                C4.setAverage(Cell2.oldAverage());
                //C4.setAverage(2*C3.average()-Cell2.average());

                // Also change the values in the boundary
                Cell3.setAverage(C3.average());
                Cell4.setAverage(C4.average());
                break;

            // FREE-SLIP WALL : Neuman on all quantities
            case 3:

                C3 = Cell2;
                C4 = Cell1;
                break;

            // ADIABATIC WALL
            case 4:

                // Dirichlet on velocity
                V3 = Cell3.velocity();
                V4 = Cell4.velocity();

                // Neuman on temperature
                T3 = Cell2.temperature();
                T4 = Cell1.temperature();

                // Neuman on pressure
                P3 = Cell2.pressure();
                P4 = Cell1.pressure();

                // Extrapolate pressure
                //P3 = 2*Cell2.pressure()-Cell1.pressure();
                //P4 = P3;

                // Compute densities
                rho3 = Gamma*Ma*Ma*P3/T3;
                rho4 = Gamma*Ma*Ma*P4/T4;

                // Compute energies
                e3 = P3/((Gamma-1.)*rho3) + 0.5*N2(V3);
                e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);

                // Correct C3 and C4
                C3.setAverage(1, rho3);
                C4.setAverage(1, rho4);
                for (i=1; i<=Dimension; i++)
                {
                    C3.setAverage(i+1, rho3*V3.value(i));
                    C4.setAverage(i+1, rho4*V4.value(i));
                }
                C3.setAverage(Dimension+2,rho3*e3);
                C4.setAverage(Dimension+2,rho4*e4);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                {
                    C3.setAverage(Dimension+3, Cell2.average(Dimension+3));
                    C4.setAverage(Dimension+3, Cell1.average(Dimension+3));
                }

                break;

            // ISOTHERMAL WALL
            case 5:

                // Dirichlet on velocity
                V3 = Cell3.velocity();
                V4 = Cell4.velocity();

                // Dirichlet on temperature
                T3 = Cell3.temperature();
                T4 = Cell4.temperature();

                // Neuman on pressure
                P3 = Cell2.pressure();
                P4 = Cell1.pressure();

                // Extrapolate pressure
                //P3 = 2*Cell2.pressure()-Cell1.pressure();
                //P4 = P3;

                // Compute densities
                rho3 = Gamma*Ma*Ma*P3/T3;
                rho4 = Gamma*Ma*Ma*P4/T4;

                // Compute energies
                e3 = P3/((Gamma-1.)*rho3) + 0.5*N2(V3);
                e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);

                // Correct C3 and C4
                C3.setAverage(1, rho3);
                C4.setAverage(1, rho4);
                for (i=1; i<=Dimension; i++)
                {
                    C3.setAverage(i+1, rho3*V3.value(i));
                    C4.setAverage(i+1, rho4*V4.value(i));
                }
                C3.setAverage(Dimension+2,rho3*e3);
                C4.setAverage(Dimension+2,rho4*e4);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                {
                    C3.setAverage(Dimension+3, Cell2.average(Dimension+3));
                    C4.setAverage(Dimension+3, Cell1.average(Dimension+3));
                }

                break;
        };
    return;
    }

  // --- Cell4 IN THE BOUNDARY, Cell3 IN THE FLUID ------------------------------------------------

    if (InCell4 != 0 && InCell3 == 0)
    {
        switch(InCell4)
        {
            // INFLOW
            case 1:
                // Dirichlet on temperature
                T4 = Cell4.temperature();

                // Extrapolate pressure from old value
                P4 = Cell3.oldPressure();

                // Compute density
                rho4 = Gamma*Ma*Ma*P4/T4;

                // Dirichlet on momentum
                V4 = (Cell4.density()/rho4)*Cell4.velocity();

                // Dirichlet on partial mass
                if (ScalarEqNb == 1)
                    Y4 = Cell4.average(Dimension+3)/Cell4.density();

                // Compute energies
                e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);

                // Correct C4
                C4.setAverage(1, rho4);
                for (i=1; i<=Dimension; i++)
                    C4.setAverage(i+1, rho4*V4.value(i));
                C4.setAverage(Dimension+2,rho4*e4);

                if (ScalarEqNb == 1)
                    C4.setAverage(Dimension+3,rho4*Y4);
                break;

            // OUTFLOW : Use old cell-average values of the neighbour
            case 2:

                C4.setAverage(Cell3.oldAverage());
                break;

            // FREE-SLIP WALL : Neuman on all quantities
            case 3:

                C4 = Cell3;
                break;

            // ADIABATIC WALL
            case 4:

                // Dirichlet on velocity
                V4 = Cell4.velocity();

                // Neuman on temperature
                T4 = Cell3.temperature();

                // Neuman on pressure
                P4 = Cell3.pressure();

                // Extrapolate pressure
                //P4 = 2*Cell3.pressure()-Cell2.pressure();

                // Compute density
                rho4 = Gamma*Ma*Ma*P4/T4;

                // Compute energy
                e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);

                // Correct C4
                C4.setAverage(1, rho4);
                for (i=1; i<=Dimension; i++)
                    C4.setAverage(i+1, rho4*V4.value(i));
                C4.setAverage(Dimension+2,rho4*e4);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                    C4.setAverage(Dimension+3, Cell3.average(Dimension+3));

                break;

            // ISOTHERMAL WALL
            case 5:

                // Dirichlet on velocity
                V4 = Cell4.velocity();

                // Dirichlet on temperature
                T4 = Cell4.temperature();

                // Neuman on pressure
                P4 = Cell3.pressure();

                // Extrapolate pressure
                //P4 = 2*Cell3.pressure()-Cell2.pressure();

                // Compute density
                rho4 = Gamma*Ma*Ma*P4/T4;

                // Compute energies
                e4 = P4/((Gamma-1.)*rho4) + 0.5*N2(V4);

                // Correct C4
                C4.setAverage(1, rho4);
                for (i=1; i<=Dimension; i++)
                    C4.setAverage(i+1, rho4*V4.value(i));
                C4.setAverage(Dimension+2,rho4*e4);

                // Neuman on partial mass
                if (ScalarEqNb == 1)
                    C4.setAverage(Dimension+3, Cell3.average(Dimension+3));

                break;
        };
    return;
    }

}

Here is the caller graph for this function:

Vector InitAverage	(	real	x,
		real	y = 0.,
		real	z = 0.
	)

Returns the initial condition in (x, y, z) form the one defined in carmen.ini.

Parameters

x	Position x
y	Position y. Defaults to 0..
z	Position z. Defaults to 0..

Returns

Vector

{
    // --- Local variables ---

    Vector Result(QuantityNb);
    real     *Q;
  Q = new real [QuantityNb+1];
    int      n;

    // --- Init Q ---

    for (n = 1; n <= QuantityNb; n++)
        Q[n]=0.;

    // --- Use definition of initial Q contained in file 'initial' ---

    #include "carmen.ini"

    // --- Fill vector Result and return it ---

    for (n = 1; n <= QuantityNb; n++)
        Result.setValue(n, Q[n]);

  delete[] Q;

    return Result;
}

Here is the caller graph for this function:

void InitParameters

(

)

Inits parameters from file carmen.par. If a parameter is not mentioned in this file, the default value is used.

Returns

void

— Compute ch --------------------------------------------------------—

{
    // --- Local variables ---------------------------------------------------------------------

    int   i;        // Counter

    // --- Set global variables from file "carmen.par" -----------------------------------------

    #include "carmen.par"

    // --- Adapt IterationNbRef to the dimension -----------------------------------------------

    IterationNbRef=(int)(exp((4.-Dimension)*log(10.)));

    // --- Compute the number of children of a given parent cell ---

    ChildNb = (1<<Dimension);

#if defined PARMPI

    AllTaskScaleNb=ScaleNb;
    for (i=0;i<4;i++)
    {
        AllXMax[i]=XMax[i];
        AllXMin[i]=XMin[i];
    }

    //some combinations give deadlock...
    MPISendType = 10;   //0 - Ibsend; 10 - Isend; 20 - Issend;
    MPIRecvType = 1;    //0 - Recv;  1 - Irecv;

    CPUScales=0;
    int tmp=size;
    while ((tmp=(tmp>>1))>0) CPUScales++;
    ScaleNb-=CPUScales/Dimension;

    one_D=1; two_D=1;
    if (Dimension >= 2) one_D=1<<ScaleNb;
    if (Dimension == 3) two_D=1<<ScaleNb;

//#if defined PARMPI

    NeighbourNb=2;
    MaxCellElementsNb=6;

    // -- Create memory arrays thats are needs for the MPI Type creation ---
    disp = new MPI_Aint[NeighbourNb*MaxCellElementsNb*one_D*two_D];
    blocklen = new int [NeighbourNb*MaxCellElementsNb*one_D*two_D];

    // --- Allocate additional memory for MPI buffer send---
    Cell tc;
    int CellElNb,bufsize;
    CellElNb=tc.size().dimension()+tc.center().dimension()+tc.average().dimension()+tc.tempAverage().dimension()+tc.divergence().dimension();

    if (EquationType==6)
        CellElNb += tc.gradient().lines()*tc.gradient().columns();

    bufsize=(CellElNb*one_D*two_D*NeighbourNb+MPI_BSEND_OVERHEAD)*2*Dimension+1024;
    MPIbuffer=new real[bufsize];
    MPI_Buffer_attach(MPIbuffer,bufsize*sizeof(real));

#else
    NeighbourNb=0;
#endif

#if defined PARMPI
    CreateMPITopology();

  // --- Compute domain coordinates for the processors ---
    XMin[1] = AllXMin[1] + coords[0]*(AllXMax[1]-AllXMin[1])/CartDims[0];
    XMax[1] = AllXMin[1] + (coords[0]+1)*(AllXMax[1]-AllXMin[1])/CartDims[0];

    if (Dimension >= 2)
    {
        XMin[2] = AllXMin[2] + coords[1]*(AllXMax[2]-AllXMin[2])/CartDims[1];
        XMax[2] = AllXMin[2] + (coords[1]+1)*(AllXMax[2]-AllXMin[2])/CartDims[1];
    }

    if (Dimension == 3)
    {
        XMin[3] = AllXMin[3] + coords[2]*(AllXMax[3]-AllXMin[3])/CartDims[2];
        XMax[3] = AllXMin[3] + (coords[2]+1)*(AllXMax[3]-AllXMin[3])/CartDims[2];
    }

    // --- Set the backup file name for the current processor
    sprintf(BackupName,"%d_%d_%d_%s",coords[0],coords[1],coords[2],"carmen.bak");
#else
    sprintf(BackupName,"%s","carmen.bak");
#endif


    // --- Use CVS only if Dimension > 1 -------------------------------------------------------

    if (Dimension == 1)
        CVS = false;

    // --- TimeAveraging always false if not Navier-Stokes -------------------------------------

    if (EquationType != 6)
        TimeAveraging = false;

    // --- If there is no file "carmen.bak", set Recovery=false --------------------------------

    if (!fopen(BackupName,"r"))
        Recovery = false;

    // --- If PrintMoreScales != 0 or 1 with Multiresolution = false, print error and stop ---

    if (!Multiresolution && (!(PrintMoreScales == 0 || PrintMoreScales == -1)) )
    {
        cout << "Parameters.cpp: In method `void InitParameters()':\n";
        cout << "Parameters.cpp: value of PrintMoreScales incompatible with FV computations\n";
        cout << "Parameters.cpp: must be 0 or -1\n";
        cout << "carmen: *** [Parameters.o] Execution error\n";
        cout << "carmen: abort execution.\n";
        exit(1);
    }

    // --- Compute global volume ---------------------------------------------------------------

    GlobalVolume = fabs(XMax[1]-XMin[1]);

    if (Dimension > 1)
        GlobalVolume *= fabs(XMax[2]-XMin[2]);

    if (Dimension > 2)
        GlobalVolume *= fabs(XMax[3]-XMin[3]);

    // --- Compute PostProcessing and DataIsBinary ---------------------------------------------

    // In 1D, use Gnuplot instead of Data Explorer

    if (Dimension == 1 && PostProcessing == 2)
        PostProcessing = 1;

    // In 2D-3D, use Data Explorer instead of Gnuplot

    if (Dimension != 1 && PostProcessing == 1)
        PostProcessing = 2;

    // --- Compute number of conservative quantities -------------------------------------------

    QuantityNb = 9;

    // --- Set the dimension of QuantityMax to QuantityNb ---------------------------------------

    QuantityMax.setDimension(QuantityNb);

    // --- Set the dimension of QuantityAverage to 4 (pressure, vorticity, entropy, volume)

    QuantityAverage.setDimension(4);

    // --- Set the dimension of IntMomentum to dimension ----------------------------------------

    IntMomentum.setDimension(Dimension);

    // --- Compute minimal space step -----------------------------------------------------------

    SpaceStep = fabs(XMax[1]-XMin[1]);

    for (i = 2; i <= Dimension; i++)
        SpaceStep = Min(SpaceStep, fabs(XMax[i]-XMin[i]));

    SpaceStep /= (1<<ScaleNb);

    // --- Compute time step from CFL if TimeStep = 0 -------------------------------------------

    if (TimeStep == 0.)
    {
        if (fabs(Eigenvalue)>0.0e-20)
            TimeStep = CFL*SpaceStep/Eigenvalue;
        else
            TimeStep = 0.0001;
    }
    else
        ConstantTimeStep = true;

    ch = CFL*SpaceStep/TimeStep;

}

Here is the caller graph for this function:

real InitResistivity	(	real	x,
		real	y = 0.,
		real	z = 0.
	)

Returns the initial resistivity condition in (x, y, z) form the one defined in carmen.eta.

Parameters

x	Position x
y	Position y. Defaults to 0..
z	Position z. Defaults to 0..

Returns

double

{
    // --- Local variables ---

    real Result=0.;
    real Res = 0.;

    #include "carmen.eta"

    Result = Res;

    return Result;
}

Here is the caller graph for this function:

void InitTimeStep

(

)

Inits time step and all the parameters which depend on it.

Returns

void

{

    // --- Init TimeStep --------------------------------------------------------------

    if (TimeStep == 0)
    {
        if(Resistivity)TimeStep = CFL*SpaceStep/(Eigenvalue + eta);
        else TimeStep = CFL*SpaceStep/Eigenvalue;
    }

    // --- Compute number of iterations -----------------------------------------------

    if (PhysicalTime != 0. && IterationNb == 0)
        IterationNb =  (int)(ceil(PhysicalTime/TimeStep));

    // --- Compute Refresh ------------------------------------------------------------

    if (Refresh == 0)
        Refresh = (int)(ceil(IterationNb/(RefreshNb*1.)));

    // --- Compute PrintEvery ---------------------------------------------------------

    if ((PrintEvery == 0)&&(ImageNb != 0))
        PrintEvery = (int)(ceil(IterationNb/(ImageNb*1.)));

    // --- Compute iterations for print  ----------------------------------------------

    if (PrintTime1 != 0.)
        PrintIt1 = (int)(ceil(PrintTime1/TimeStep));

    if (PrintTime2 != 0.)
        PrintIt2 = (int)(ceil(PrintTime2/TimeStep));

    if (PrintTime3 != 0.)
        PrintIt3 = (int)(ceil(PrintTime3/TimeStep));

    if (PrintTime4 != 0.)
        PrintIt4 = (int)(ceil(PrintTime4/TimeStep));

    if (PrintTime5 != 0.)
        PrintIt5 = (int)(ceil(PrintTime5/TimeStep));

    if (PrintTime6 != 0.)
        PrintIt6 = (int)(ceil(PrintTime6/TimeStep));

    // --- Compute FV reference time -----------------------------------------------------------

    if (Multiresolution)
        FVTimeRef = CPUTimeRef(IterationNbRef,ScaleNbRef);

}

Here is the caller graph for this function:

void InitTree

(

Node *

Root

)

Inits tree structure from initial condition, starting form the node Root. Only for multiresolution computations.

Parameters

Root

Root

Returns

void

{
    // --- Local variables ---

    int l;  // Counter on levels

    // --- Init cell-average value in root and split it ---

    if (Recovery && UseBackup)

         Root->restore();

    else
    {
        Root->initValue();

        // --- Add and init nodes in different levels, when necessary ---

        for (l=1; l <= ScaleNb; l++)
            Root->addLevel();
    }

    // -- Check if tree is graded ---

    if (debug) Root->checkGradedTree();

}

Here is the caller graph for this function:

Vector Limiter	(	const Vector	u,
		const Vector	v
	)

Returns the value of the slope limiter between the slopes u and v.

Parameters

u	Vector
v	Vector

Returns

Vector

{
        // Min Mod limiter

        int LimiterNo = 3;

        Vector Result(u.dimension());
        int i;
        real x, y;  // slopes

    for (i=1; i<=u.dimension(); i++)
        {
            x = u.value(i);
            y = v.value(i);

            switch(LimiterNo)
            {
                // MIN-MOD
                case 1:
                    if (x == y)
                        Result.setValue(i, 0.);
                    else
                        Result.setValue(i,Min(1., fabs(x)/fabs(x-y)));
                    break;

        // VAN LEER
                case 3:
                default:
                    if ((fabs(x) + fabs(y)) == 0.)
                        Result.setValue(i, 0.);
                    else
                        Result.setValue(i,fabs(x)/(fabs(x)+fabs(y)));
                    break;
            };
        }

        return Result;
}

Here is the caller graph for this function:

real Limiter

(

const real

x

)

Returns the valur of slope limiter from a real value x.

Parameters

x

...

Returns

double

{
  real Result = 0.;

  switch (LimiterNo)
  {
    case 1: // Min-Mod
        Result = Max(0., Min(1.,r));
        break;

    case 2: // Van Albada
        Result = (r<=0.)? 0.: (r*r+r)/(r*r+1.);
        break;

    case 3: // Van Leer
        Result = (r<=0.) ? 0.: (r+Abs(r))/(1.+Abs(r));
        break;

    case 4: // Superbee
        Result = (r<=0.) ? 0.: Max(0.,Max(Min(2.*r,1.),Min(r,2.)));
        break;
    case 5: // Monotonized Central
        Result = max(0.0,min(min(2*r,0.5*(1+r)),2.0));
        break;

  };

        return Result;
}

real MinAbs	(	const real	a,
		const real	b
	)

Returns the minimum in module of a and b.

Parameters

a	Real value
b	Real value

Returns

double

{
    return (fabs(a) <= fabs(b))? a:b;
}

real NormMaxQuantities

(

const Vector &

V

)

Returns the Max-norm of the vector where every quantity is divided by its characteristic value.

Parameters

V

Vector

Returns

double

{
    Vector W(QuantityNb);
    int AxisNo=1;
    real MomentumMax=0.;
    real MagMax=0.;



    W.setZero();

/*
    // Density
    W.setValue(1, V.value(1)/QuantityMax.value(1));

    // Momentum
    W.setValue(2 , V.value(2)/QuantityMax.value(2) );
    W.setValue(3 , V.value(3)/QuantityMax.value(3) );
    W.setValue(4 , V.value(4)/QuantityMax.value(4) );

    // Energy
    W.setValue(5 , V.value(5)/QuantityMax.value(5) );

    // psi
   // W.setValue(6, V.value(6)/QuantityMax.value(6));

    // Magnetic Field
    W.setValue(7 , V.value(7)/QuantityMax.value(7) );
    W.setValue(8 , V.value(8)/QuantityMax.value(8) );
    W.setValue(9 , V.value(9)/QuantityMax.value(9) );
*/

  // --- Compute Linf norm --

    W.setValue(1, (V.value(1))/QuantityMax.value(1));
    W.setValue(5, (V.value(5))/QuantityMax.value(5));
    //W.setValue(6, (V.value(6))/QuantityMax.value(6));


    for (AxisNo = 1; AxisNo <= Dimension; AxisNo++)
    {
        MomentumMax = Max( MomentumMax, QuantityMax.value(AxisNo+1) );
        MagMax      = Max( MagMax,      QuantityMax.value(AxisNo+6) );
        W.setValue(2, W.value(2) + V.value(AxisNo+1)*V.value(AxisNo+1));
        W.setValue(7, W.value(7) + V.value(AxisNo+6)*V.value(AxisNo+6));
    }

    if(Dimension==2){
        W.setValue(4, (V.value(4))/QuantityMax.value(4));
        W.setValue(9, (V.value(9))/QuantityMax.value(9));
    }

    W.setValue(2, sqrt(W.value(2))/MomentumMax);
    W.setValue(7, sqrt(W.value(7))/MagMax );

    if(IterationNo==0)return NMax(V);
    return NMax(W);
}

void Performance

(

const char *

FileName

)

Computes the performance of the computation and, for cluster computations, write it into file FileName.

Parameters

FileName

Name of the file.

Returns

void

{
    // --- Local variables ---

    bool EndComputation;            // True if end computation
    int FineCellNb;             // Number of cells on fine grid
    int CellPVirt;
    FILE    *output;            // Pointer to output file

    double realtimefull;            //full real time
    double ftime;               // real time
    double ctime;               // CPU time
    unsigned int    ttime, rtime;       // total and remaining  real time (in seconds)
    unsigned int  tctime=0, rctime=0;   // total and remaining CPU time (in seconds)
    unsigned int    rest;
    int day, hour, min, sec;

    // --- Init EndComputation

    EndComputation = (IterationNo > IterationNb);

    // --- Compute FineCellNb ---

    FineCellNb = 1<<(ScaleNb*Dimension);
    CellPVirt = 1<<(ScaleNb*(Dimension-1));
    CellPVirt = CellPVirt*2*Dimension;
    // --- Write in file ---

/*
    char CPUFileName[255];
#if defined PARMPI
    sprintf(CPUFileName,"%d_%d_%d_%s",coords[0],coords[1],coords[2],FileName);
//  strcpy(CPUFileName, FileName);
#else
    strcpy(CPUFileName, FileName);
#endif
*/

    if ((output = fopen(FileName,"w")))
    {

        realtimefull=ftime = CPUTime.realTime();
        ctime = CPUTime.CPUTime();

        if (!EndComputation)
        {
            ttime  = (unsigned int)((ftime*IterationNb)/IterationNo);
            rtime  = (unsigned int)((ftime*(IterationNb-IterationNo))/IterationNo);
            tctime = (unsigned int)((ctime*IterationNb)/IterationNo);
            rctime = (unsigned int)((ctime*(IterationNb-IterationNo))/IterationNo);
    }

        fprintf(output, "Dimension             : %12i\n", Dimension);

        if (EndComputation)
            fprintf(output, "Iterations            : %12i\n", IterationNb);
        else
        {
            fprintf(output, "Iterations (total)    : %12i\n", IterationNb);
            fprintf(output, "Iterations (elapsed)  : %12i\n", IterationNo);
            fprintf(output, "In progress           :%13.6f %%\n", 100.*IterationNo/(1.*IterationNb));
    }

        fprintf(output, "Scales (max)          : %12i\n", ScaleNb);
        fprintf(output, "Cells (max)           : %12i\n", (1<<(ScaleNb*Dimension)));

        if (Multiresolution)
          fprintf(output, "Solver                :           MR\n");
        else
          fprintf(output, "Solver                :           FV\n");

        //fprintf(output, "Time integration      :     explicit\n");
        fprintf(output, "Time accuracy order   : %12i\n", StepNb);
        fprintf(output, "Time step             :%13.6e s\n", TimeStep);
        fprintf(output, "Threshold parameter   :%13.6e \n", Tolerance);
        fprintf(output, "Threshold norm        : %12i\n", ThresholdNorm);
        fprintf(output, "CFL                   :%13.6e \n", CFL);

        if(Resistivity)
            fprintf(output, "Eta                   :%13.6e \n", eta);
        if(Diffusivity)
            fprintf(output, "Chi                   :%13.6e \n", chi);

        if (EndComputation)
        {
            fprintf(output, "Physical time         :%13.6e s\n", ElapsedTime); //TimeStep * IterationNb);
            fprintf(output, "CPU time (s)          :%13.6e s\n", ctime);

            if (Multiresolution)
                fprintf(output, "CPU time / it. x pt   :%13.6e s\n", ctime/TotalLeafNb);
            else
                    fprintf(output, "CPU time / it. x pt   :%13.6e s\n", ctime/((1<<(ScaleNb*Dimension))*IterationNb));

            if (Multiresolution)
            {

                fprintf(output, "Leaf compression      :%13.6f %% \n", (100.*TotalLeafNb)/(1.0*IterationNb*FineCellNb));
                //fprintf(output, "Memory compression    :%13.6f %% \n", (100.*TotalCellNb)/(1.0*IterationNb*(FineCellNb)));
                fprintf(output, "Memory compression    :%13.6f %% \n", (100.*TotalCellNb)/(1.0*IterationNb*(FineCellNb + CellPVirt)));
                fprintf(output, "CPU compression       :%13.6f %% \n", (100.*ctime)/(IterationNb*FVTimeRef));
            }
            else
            {
                fprintf(output, "Leaf compression      :%13.6f %% \n", 100.);
                fprintf(output, "Memory compression    :%13.6f %% \n", 100.);
                fprintf(output, "CPU compression       :%13.6f %% \n", 100.);
            }
    }
        else
        {
            fprintf(output, "Total physical time   :%13.6e s\n", TimeStep * IterationNb);
            fprintf(output, "Elapsed physical time :%13.6e s\n", TimeStep * IterationNo);
            if (Multiresolution)
            {
                fprintf(output, "Leaf compression      :%13.6f %% \n", (100.*TotalLeafNb)/(1.0*IterationNb*FineCellNb));
                fprintf(output, "Memory compression    :%13.6f %% \n", (100.*TotalCellNb)/(1.0*IterationNo*FineCellNb));
                fprintf(output, "CPU compression       :%13.6f %% \n", (100.*ctime)/(IterationNo*FVTimeRef));
            }
            else
            {
                fprintf(output, "Leaf compression      :%13.6f %% \n", 100.);
                fprintf(output, "Memory compression    :%13.6f %% \n", 100.);
                fprintf(output, "CPU compression       :%13.6f %% \n", 100.);
            }
        }

        fprintf(output, "\n");

        if (EndComputation)
        {
            // --- Print final time ------------------------------------------------

            rest =  (unsigned int)(ctime);
        day  =  rest/86400;
            rest %= 86400;
            hour =  rest/3600;
            rest %= 3600;
            min  =  rest/60;
            rest %= 60;
            sec  =  rest;
            rest =  (unsigned int)(ctime);

            if (rest >= 86400)
                fprintf(output, "CPU time : %5d day %2d h %2d min %2d s\n", day, hour, min, sec);

            if ((rest < 86400)&&(rest >= 3600))
                fprintf(output, "CPU time : %2d h %2d min %2d s\n", hour, min, sec);

            if ((rest < 3600)&&(rest >= 60))
                fprintf(output, "CPU time : %2d min %2d s\n", min, sec);

            if (rest < 60)
                fprintf(output, "CPU time : %2d s\n", sec);
        }
        else
        {
            // --- Print total time ------------------------------------------------

            rest =  tctime;
        day  =  rest/86400;
            rest %= 86400;
            hour =  rest/3600;
            rest %= 3600;
            min  =  rest/60;
            rest %= 60;
            sec  =  rest;

            if (tctime >= 86400)
                fprintf(output, "Total CPU time     (estimation) : %5d day %2d h %2d min %2d s\n", day, hour, min, sec);

            if ((tctime < 86400)&&(tctime >= 3600))
                fprintf(output, "Total CPU time     (estimation) : %2d h %2d min %2d s\n", hour, min, sec);

            if ((tctime < 3600)&&(tctime >= 60))
                fprintf(output, "Total CPU time     (estimation) : %2d min %2d s\n", min, sec);

            if (tctime < 60)
                fprintf(output,"Total CPU time     (estimation) : %2d s\n", sec);

            // --- Print remaining time ------------------------------------------------

            rest =  rctime;
        day  =  rest/86400;
            rest %= 86400;
            hour =  rest/3600;
            rest %= 3600;
            min  =  rest/60;
            rest %= 60;
            sec  =  rest;

            if (rctime >= 86400)
                fprintf(output, "Remaining CPU time (estimation) : %5d day %2d h %2d min %2d s\n", day, hour, min, sec);

            if ((rctime < 86400)&&(rctime >= 3600))
                fprintf(output, "Remaining CPU time (estimation) : %2d h %2d min %2d s\n", hour, min, sec);

            if ((rctime < 3600)&&(rctime >= 60))
                fprintf(output, "Remaining CPU time (estimation) : %2d min %2d s\n", min, sec);

            if (rctime < 60)
                fprintf(output, "Remaining CPU time (estimation) : %2d s\n", sec);

        }


#if defined PARMPI
    fprintf(output,  "\n");
    fprintf(output,  "Real time (time() function) :%lf\n",  realtimefull);
        fprintf(output,  "clock() function :%lf\n",  ctime);
        fprintf(output,  "\nCommunications real timer: %lf\n",  CommTimer.realTime());
        fprintf(output,  "Communications clock():%lf\n",  CommTimer.CPUTime());
#endif

        fprintf(output, "\n");
        fclose (output);
    }
    else
    {
        cout << "Performance.cpp: In method `void Performance(Node*, char*)':\n";
        cout << "Performance.cpp: cannot open file " << FileName << '\n';
        cout << "carmen: *** [Performance.o] Execution error\n";
        cout << "carmen: abort execution.\n";
        exit(1);
    }
}

Here is the caller graph for this function:

void PrintIntegral

(

const char *

FileName

)

Writes the integral values, like e.g flame velocity, global error, into file FileName.

Parameters

FileName

Name of the file

Returns

void

{
    // --- Local variables ---

    real    t;              // time
    FILE    *output;  // output file
    int i;                  // counter
    real Volume=1.; // Total volume

    // --- Open file ---

    if ( (IterationNo == 0) ? (output = fopen(FileName,"w")) : (output = fopen(FileName,"a")) )
    {
    // HEADER

        if (IterationNo == 0)
        {

            fprintf(output, "#");
            fprintf(output, TXTFORMAT, " Time");
            fprintf(output, TXTFORMAT, "CFL");
            fprintf(output, TXTFORMAT, "Energy");
            fprintf(output, TXTFORMAT, "Div B Max");
            fprintf(output, TXTFORMAT, "ch");
            fprintf(output, TXTFORMAT, "Helicity");
            fprintf(output, TXTFORMAT, "DivB Max norm");
/*
            if (Multiresolution)
            {
                fprintf(output, TXTFORMAT, "Memory comp.");
                fprintf(output, TXTFORMAT, "CPU comp.");
                if (ExpectedCompression != 0. || CVS)
                    fprintf(output, TXTFORMAT, "Tolerance");
            //  if (CVS)
            //      fprintf(output, TXTFORMAT, "Av. Pressure");

            }
*/
            if (!ConstantTimeStep)
            {
                if (StepNb == 3) fprintf(output, TXTFORMAT, "RKF Error");
                fprintf(output, TXTFORMAT,"Next time step");
                fprintf(output, "%13s  ", "IterationNo");
                fprintf(output, "%13s  ", "IterationNb");
            }

            fprintf(output,"\n");

        }

        if (ConstantTimeStep)
            t=IterationNo*TimeStep;
        else
            t = ElapsedTime;

        fprintf(output, FORMAT, t);

        // --- Compute total volume ---

        for (i=1; i<= Dimension; i++)
            Volume *= fabs(XMax[i]-XMin[i]);

        // Print CFL
        fprintf(output, FORMAT,Eigenvalue*TimeStep/SpaceStep);

        // Print momentum and energy
        //fprintf(output, FORMAT, GlobalMomentum);
        fprintf(output, FORMAT, GlobalEnergy);
        fprintf(output, FORMAT, DIVBMax);
        fprintf(output, FORMAT, ch);
        fprintf(output, FORMAT, Helicity);
        fprintf(output, FORMAT, DIVB);


/*
        if (Multiresolution)
        {
            fprintf(output, FORMAT, (1.*CellNb)/(1<<(ScaleNb*Dimension)));
            fprintf(output, FORMAT, CPUTime.CPUTime()/(IterationNo*FVTimeRef));

            if (ExpectedCompression != 0.)
                fprintf(output, FORMAT, Tolerance);

            //  if (CVS)
            //{
        //      fprintf(output, FORMAT, ComputedTolerance(ScaleNb));
            //  fprintf(output, FORMAT, QuantityAverage.value(1));
            //}
        }
*/
        if (!ConstantTimeStep)
        {
            if (StepNb == 3) fprintf(output, FORMAT, RKFError);
            fprintf(output, FORMAT, TimeStep);
            fprintf(output, "%13i  ", IterationNo);
            fprintf(output, "%13i  ", IterationNb);
        }

        fprintf(output, "\n");
        fclose(output);
    }
    else
    {
        cout << "PrintIntegral.cpp: In method `void PrintIntegral(Node*, char*)´:\n";
        cout << "PrintIntegral.cpp: cannot open file " << FileName << '\n';
        cout << "carmen: *** [PrintIntegral.o] Execution error\n";
        cout << "carmen: abort execution.\n";
        exit(1);
    }
}

Here is the caller graph for this function:

void ReduceIntegralValues

(

)

Parallel function DOES NOT WORK!

Returns

void

                            {
real rb; //Recieve Buffer
rb=0.0;
  CommTimer.start();
#if defined PARMPI
  MPI_Reduce(&ErrorMax,&rb,1,MPI_Type,MPI_MAX,0,MPI_COMM_WORLD);
  ErrorMax=rb;

  MPI_Reduce(&ErrorMid,&rb,1,MPI_Type,MPI_SUM,0,MPI_COMM_WORLD);
  ErrorMid=rb/size;

  MPI_Reduce(&ErrorL2,&rb,1,MPI_Type,MPI_SUM,0,MPI_COMM_WORLD);
  ErrorL2=rb/size;

  MPI_Reduce(&ErrorNb,&rb,1,MPI_Type,MPI_SUM,0,MPI_COMM_WORLD);
  ErrorNb=rb;

  MPI_Allreduce(&FlameVelocity,&rb,1,MPI_Type,MPI_SUM,MPI_COMM_WORLD);
  FlameVelocity=rb;

  MPI_Allreduce(&GlobalMomentum,&rb,1,MPI_Type,MPI_SUM,MPI_COMM_WORLD);
  GlobalMomentum=rb;

  MPI_Allreduce(&GlobalEnergy,&rb,1,MPI_Type,MPI_SUM,MPI_COMM_WORLD);
  GlobalEnergy=rb;

  MPI_Reduce(&ExactMomentum,&rb,1,MPI_Type,MPI_SUM,0,MPI_COMM_WORLD);
  ExactMomentum=rb;

  MPI_Reduce(&ExactEnergy,&rb,1,MPI_Type,MPI_SUM,0,MPI_COMM_WORLD);
  ExactEnergy=rb;

  MPI_Allreduce(&GlobalReactionRate,&rb,1,MPI_Type,MPI_SUM,MPI_COMM_WORLD);
  GlobalReactionRate=rb;

  MPI_Allreduce(&EigenvalueMax, &rb,1,MPI_Type,MPI_MAX,MPI_COMM_WORLD);
  EigenvalueMax=rb;

#endif
  CommTimer.stop();
}

Here is the caller graph for this function:

void RefreshTree

(

Node *

Root

)

Refresh the tree structure, i.e. compute the cell-averages of the internal nodes by projection and those of the virtual leaves by prediction. The root node is Root. Only for multiresolution computations.

Parameters

Root

Root

Returns

void

{
    // --- Project : compute cell-average values in all nodes ---
    Root->project();

    // --- Fill virtual children with predicted values ---
    Root->fillVirtualChildren();
}

Here is the caller graph for this function:

void Remesh

(

Node *

Root

)

Remesh the tree structure after a time evolution. The root node is Root. Only for multiresolution computations.

Parameters

Root

Root

Returns

void

{
    // --- Refresh tree structure ---
//  RefreshTree(Root);

    // --- Check if tree is graded ---
    if (debug) Root->checkGradedTree();

    // --- Adapt : depending on details, refine or combine ---
    Root->adapt();

    // --- Check if tree is graded ---
    if (debug) Root->checkGradedTree();
}

Here is the caller graph for this function:

Vector ResistiveTerms	(	Cell &	Cell1,
		Cell &	Cell2,
		Cell &	Cell3,
		Cell &	Cell4,
		int	AxisNo
	)

Returns the resistive source terms in the cell UserCell.

Parameters

Cell1	Cell 1
Cell2	Cell 2
Cell3	Cell 3
Cell4	Cell 4
AxisNo	Axis of interest

Returns

Vector

X - direction

2D

Y - direction

2D

Z - direction

3D

{
    // ---  Local variables ---
    Vector  B(3), Bi(3), Bj(3), Bk(3);
    Vector  Result(QuantityNb);
    Vector  Bavg(3);
    real    Jx = 0., Jy = 0., Jz = 0.;
    real    dx, dy, dz;
    real    ResX= 0., ResY= 0., ResZ= 0., ResE= 0.;
    real    eta0=0., etai=0., etaj=0., etak=0., etaR=0.;

    dx = Cell2.size(1);
    dy = Cell2.size(2);
    dz = Cell2.size(3);

    eta0 = Cell1.Res;
    etai = Cell2.Res;
    etaj = Cell3.Res;
    etak = Cell4.Res;

    for(int i=1; i <= 3; i++){
        B.setValue (i, Cell1.average(i+6));
        Bi.setValue(i, Cell2.average(i+6));
        Bj.setValue(i, Cell3.average(i+6));
        Bk.setValue(i, Cell4.average(i+6));
    }


    if(AxisNo == 1){
        etaR = (eta0 + etai)/2.;

        Bavg.setValue(2, 0.5*(B.value(2) + Bi.value(2)));
        Bavg.setValue(3, 0.5*(B.value(3) + Bi.value(3)));

        Jy = -(B.value(3) - Bi.value(3))/dx;
        Jz =  (B.value(2) - Bi.value(2))/dx;

        Jz = Jz - (B.value(1) - Bj.value(1))/dy;

        if(Dimension==3){
            Jy   = Jy   + (B.value(1) - Bk.value(1))/dz;
        }

        ResE =  etaR*(Bavg.value(2)*Jz - Bavg.value(3)*Jy);
        ResX = 0.;
        ResY =  etaR*Jz;
        ResZ = -etaR*Jy;


    }else if(AxisNo == 2){
        etaR = (eta0 + etaj)/2.;

        Bavg.setValue(1, 0.5*(B.value(1) + Bj.value(1)));
        Bavg.setValue(3, 0.5*(B.value(3) + Bj.value(3)));

        Jx =  (B.value(3) - Bj.value(3))/dy;
        Jz = -(B.value(1) - Bj.value(1))/dy;

        Jz = Jz + (B.value(2) - Bi.value(2))/dx;

        if(Dimension==3){
            Jx   = Jx   + (B.value(2) - Bk.value(2))/dz;
        }

        ResE =  etaR*(Bavg.value(3)*Jx - Bavg.value(1)*Jz);
        ResX = -etaR*Jz;
        ResY = 0.;
        ResZ =  etaR*Jx;

    }else{
        etaR = (eta0 + etak)/2.;

        Bavg.setValue(1, 0.5*(B.value(1) + Bk.value(1)));
        Bavg.setValue(2, 0.5*(B.value(2) + Bk.value(2)));

        Jx = -(B.value(2) - Bk.value(2))/dz;
        Jy =  (B.value(1) - Bk.value(1))/dz;

        Jx = Jx + (B.value(3) - Bj.value(3))/dy;
        Jy = Jy - (B.value(3) - Bi.value(3))/dx;

        ResE =  etaR*(Bavg.value(1)*Jy - Bavg.value(2)*Jx);
        ResX =  etaR*Jy;
        ResY = -etaR*Jx;
        ResZ = 0.;
    }


    Result.setZero();

    // These values will be added to the numerical flux
    Result.setValue(5, ResE);
    Result.setValue(7, ResX);
    Result.setValue(8, ResY);
    Result.setValue(9, ResZ);

    return Result;

}

Here is the caller graph for this function:

Vector SchemeHLL	(	const Cell &	Cell1,
		const Cell &	Cell2,
		const Cell &	Cell3,
		const Cell &	Cell4,
		const int	AxisNo
	)

Returns the HLL numerical flux for MHD equations. The scheme uses four cells to estimate the flux at the interface. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides.

Parameters

Cell1	second neighbour on the left side
Cell2	first neighbour on the left side
Cell3	first neighbour on the right side
Cell4	second neighbour on the right side
AxisNo	Axis of interest.

Returns

Vector

{

    // General variables

    Vector  LeftAverage(QuantityNb);    //
    Vector  RightAverage(QuantityNb);   // Conservative quantities
    Vector  Result(QuantityNb);         // MHD numerical flux
    int aux=0;
    // Variables for the HLL scheme
    Vector  FL(QuantityNb), FR(QuantityNb); //Left and right physical fluxes
    Vector  VL(3), VR(3);   // Left and right velocities
    Vector  BL(3), BR(3);   // Left and right velocities
    real    rhoL=0., rhoR=0.;       // Left and right densities
    real    eL=0., eR=0.;           // Left and right energies
    real    preL=0., preR=0.;
    real    bkL=0., bkR=0.;
    real    aL=0., aR=0.;
    real    bL=0., bR=0.;
    real    cfL=0., cfR=0.;
    real    SL=0., SR=0.;
    real    dx=0.;
    dx = Cell2.size(AxisNo);
    real   r, Limit, LeftSlope = 0., RightSlope = 0.; // Left and right slopes
    int i;

//  --- Limiter function ---------------------------------------------------------

    for (i=1; i<=QuantityNb; i++)
    {
        // --- Compute left cell-average value ---

        if (Cell2.average(i) != Cell1.average(i))
        {
            RightSlope  = Cell3.average(i)-Cell2.average(i);
            LeftSlope   = Cell2.average(i)-Cell1.average(i);
            r           = RightSlope/LeftSlope;
            Limit       = Limiter(r);
            LeftAverage.setValue(i, Cell2.average(i) + 0.5*Limit*LeftSlope);
            aux = 1;
        }
        else
            LeftAverage.setValue(i, Cell2.average(i));

        // --- Compute right cell-average value ---

        if (Cell3.average(i) != Cell2.average(i))
        {
            RightSlope  = Cell4.average(i)-Cell3.average(i);
            LeftSlope   = Cell3.average(i)-Cell2.average(i);
            r           = RightSlope/LeftSlope;
            Limit       = Limiter(r);
            RightAverage.setValue(i, Cell3.average(i) - 0.5*Limit*LeftSlope);
            aux = 1;
        }
        else
            RightAverage.setValue(i, Cell3.average(i));
    }


    // --- HLL scheme -------------------------------------------------------------

    // --- Conservative variables ---

    // Left and right densities
    rhoL = LeftAverage.value(1);
    rhoR = RightAverage.value(1);

    // Left and right momentum and magnetic field
    for (int i=1;i<=3;i++)
    {
        VL.setValue( i, LeftAverage.value(i+1));
        VR.setValue( i, RightAverage.value(i+1));
        BL.setValue( i, LeftAverage.value(i+6));
        BR.setValue( i, RightAverage.value(i+6));
    }

    // Left and right energies
    eL = LeftAverage.value(5);
    eR = RightAverage.value(5);

    // Left and right pressures
    preL = (Gamma -1.)*(eL - 0.5*(VL*VL)/rhoL - 0.5*(BL*BL));
    preR = (Gamma -1.)*(eR - 0.5*(VR*VR)/rhoR - 0.5*(BR*BR));

    // --- Magnetoacoustic waves calculations --

    bkL = power2(BL.value(AxisNo))/rhoL;
    bkR = power2(BR.value(AxisNo))/rhoR;

    aL = Gamma*preL/rhoL;
    aR = Gamma*preR/rhoR;

    bL = (BL*BL)/rhoL;
    bR = (BR*BR)/rhoR;

    // Left and Right fast speeds
    cfL = sqrt(0.5*(aL + bL + sqrt(power2(aL + bL) - 4.0*aL*bkL)));
    cfR = sqrt(0.5*(aR + bR + sqrt(power2(aR + bR) - 4.0*aR*bkR)));

    // Left and right slopes
    SL = Min(Min(VL.value(AxisNo)/rhoL - cfL, VR.value(AxisNo)/rhoR - cfR),0.0);
    SR = Max(Max(VL.value(AxisNo)/rhoL + cfL, VR.value(AxisNo)/rhoR + cfR),0.0);

    // --- Physical flux ---
     if(AxisNo ==1){
        EigenvalueX = Max(Max(Abs(SL),Abs(SR)),EigenvalueX);
        FL = FluxX(LeftAverage);
        FR = FluxX(RightAverage);
     }else if(AxisNo ==2){
        EigenvalueY = Max(Max(Abs(SL),Abs(SR)),EigenvalueY);
        FL = FluxY(LeftAverage);
        FR = FluxY(RightAverage);
    }else{
        EigenvalueZ = Max(Max(Abs(SL),Abs(SR)),EigenvalueZ);
        FL = FluxZ(LeftAverage);
        FR = FluxZ(RightAverage);
    }


    // --- HLL Riemann Solver ---

    for(int i=1;i<=QuantityNb;i++)
    {
        Result.setValue(i, (SR*FL.value(i) - SL*FR.value(i) + SR*SL*(RightAverage.value(i) - LeftAverage.value(i)))/(SR-SL));
    }

    // Parabolic-Hyperbolic divergence Cleaning (Dedner, 2002)
    fluxCorrection(Result, LeftAverage, RightAverage, AxisNo);

    // Artificial diffusion terms
    if(Diffusivity && aux==1) Result = Result - ArtificialViscosity(LeftAverage,RightAverage,dx,AxisNo);

    return Result;

}

Here is the caller graph for this function:

Vector SchemeHLLD	(	const Cell &	Cell1,
		const Cell &	Cell2,
		const Cell &	Cell3,
		const Cell &	Cell4,
		const int	AxisNo
	)

Returns the HLLD numerical flux for MHD equations. The scheme uses four cells to estimate the flux at the interface. Cell2 and Cell3 are the first neighbours on the left and right sides. Cell1 and Cell4 are the second neighbours on the left and right sides.

Parameters

Cell1	second neighbour on the left side
Cell2	first neighbour on the left side
Cell3	first neighbour on the right side
Cell4	second neighbour on the right side
AxisNo	Axis of interest.

Returns

Vector

{

    // General variables

    Vector  LeftAverage(QuantityNb);    //
    Vector  RightAverage(QuantityNb);   // Conservative quantities
    Vector  Result(QuantityNb);         // MHD numerical flux
    int aux=0;

    // Variables for the HLL scheme
    Vector  FL(QuantityNb), FR(QuantityNb); //Left and right physical fluxes
    Vector  VL(3), VR(3);   // Left and right velocities
    Vector  BL(3), BR(3);   // Left and right velocities
    real    rhoL=0., rhoR=0.;       // Left and right densities
    real    eL=0., eR=0.;           // Left and right energies
    real    preL=0., preR=0.;
    real    bkL=0., bkR=0.;
    real    aL=0., aR=0.;
    real    bL=0., bR=0.;
    real    cfL=0., cfR=0.;
    real    SL=0., SR=0.;
    real    SLS=0., SRS=0.;
    real    SM=0.;
    Matrix U(QuantityNb,4);
    Matrix F(QuantityNb,2);
    real    dx=0.;
    dx = Cell2.size(AxisNo);
    real   r, Limit, LeftSlope = 0., RightSlope = 0.; // Left and right slopes
    int i;

//  --- Limiter function ---------------------------------------------------------

    for (i=1; i<=QuantityNb; i++)
    {
        // --- Compute left cell-average value ---

        if (Cell2.average(i) != Cell1.average(i))
        {
            RightSlope  = Cell3.average(i)-Cell2.average(i);
            LeftSlope   = Cell2.average(i)-Cell1.average(i);
            r           = RightSlope/LeftSlope;
            Limit       = Limiter(r);
            LeftAverage.setValue(i, Cell2.average(i) + 0.5*Limit*LeftSlope);
            aux = 1;
        }
        else
            LeftAverage.setValue(i, Cell2.average(i));

        // --- Compute right cell-average value ---

        if (Cell3.average(i) != Cell2.average(i))
        {
            RightSlope  = Cell4.average(i)-Cell3.average(i);
            LeftSlope   = Cell3.average(i)-Cell2.average(i);
            r           = RightSlope/LeftSlope;
            Limit       = Limiter(r);
            RightAverage.setValue(i, Cell3.average(i) - 0.5*Limit*LeftSlope);
            aux = 1;
        }
        else
            RightAverage.setValue(i, Cell3.average(i));
    }

    // --- HLLD scheme -------------------------------------------------------------

    // --- Conservative variables ---

    // Left and right densities
    rhoL = LeftAverage.value(1);
    rhoR = RightAverage.value(1);

    // Left and right momentum and magnetic field
    for (int i=1;i<=3;i++)
    {
        VL.setValue( i, LeftAverage.value (i+1));
        VR.setValue( i, RightAverage.value(i+1));
        BL.setValue( i, LeftAverage.value (i+6));
        BR.setValue( i, RightAverage.value(i+6));
    }

    // Left and right energies
    eL = LeftAverage.value(5);
    eR = RightAverage.value(5);

    // Left and right pressures
    preL = (Gamma -1.)*(eL - 0.5*(VL*VL)/rhoL - 0.5*(BL*BL));
    preR = (Gamma -1.)*(eR - 0.5*(VR*VR)/rhoR - 0.5*(BR*BR));

    // --- Magnetoacoustic waves computation --
    bkL = power2(BL.value(AxisNo))/rhoL;
    bkR = power2(BR.value(AxisNo))/rhoR;

    aL = Gamma*preL/rhoL;
    aR = Gamma*preR/rhoR;

    bL = (BL*BL)/rhoL;
    bR = (BR*BR)/rhoR;

    // Left and Right fast speeds
    cfL = sqrt(0.5*(aL + bL + sqrt(power2(aL + bL) - 4.0*aL*bkL)));
    cfR = sqrt(0.5*(aR + bR + sqrt(power2(aR + bR) - 4.0*aR*bkR)));

    // Left and Right slopes
    SL = Min((VL.value(AxisNo))/rhoL, (VR.value(AxisNo))/rhoR) - Max(cfL,cfR);
    SR = Max((VL.value(AxisNo))/rhoL, (VR.value(AxisNo))/rhoR) + Max(cfL,cfR);

    // --- Physical flux ---
    if(AxisNo ==1){
        EigenvalueX = Max(Max(Abs(SL),Abs(SR)),EigenvalueX);
        FL = FluxX(LeftAverage);
        FR = FluxX(RightAverage);
    }else if(AxisNo ==2){
        EigenvalueY = Max(Max(Abs(SL),Abs(SR)),EigenvalueY);
        FL = FluxY(LeftAverage);
        FR = FluxY(RightAverage);
    }else{
        EigenvalueZ = Max(Max(Abs(SL),Abs(SR)),EigenvalueZ);
        FL = FluxZ(LeftAverage);
        FR = FluxZ(RightAverage);
    }

    // Intermediary states U* and U**
    U = stateUstar(LeftAverage, RightAverage, preL, preR, SL, SR, SM, SLS, SRS,AxisNo);

    // --- HLLD Riemann Solver ---

    for(int i=1;i<=QuantityNb;i++)
    {
                //Flux Function - Equation 66
                //F_L
                if(SL>=0.)
                    Result.setValue(i, FL.value(i));
                //F-star left //  FL=FLstar
                else if(SLS>=0. && SL<0.)
                    Result.setValue(i, FL.value(i) + SL*(U.value(i,1) - LeftAverage.value(i)));
                //F-star-star left
                else if(SM>=0. && SLS<0.)
                    Result.setValue(i, FL.value(i) + SLS*U.value(i,3) - (SLS - SL)*U.value(i,1) - SL*LeftAverage.value(i));
                //F-star-star right
                else if(SRS>=0. && SM<0.)
                    Result.setValue(i, FR.value(i) + SRS*U.value(i,4) - (SRS - SR)*U.value(i,2) - SR*RightAverage.value(i));
                //F-star right
                else if(SR>=0. && SRS<0.)
                    Result.setValue(i, FR.value(i) + SR*(U.value(i,2) - RightAverage.value(i)));
                //F_R
                else
                    Result.setValue(i, FR.value(i));
    }
    // Parabolic-Hyperbolic divergence Cleaning (Dedner, 2002)
    //fluxCorrection(Result, Cell2.average(), Cell3.average(), AxisNo);
    fluxCorrection(Result, LeftAverage, RightAverage, AxisNo);

    // Artificial diffusion terms
    if(Diffusivity && aux==1) Result = Result - ArtificialViscosity(LeftAverage,RightAverage,dx,AxisNo);

    return Result;

}

Here is the caller graph for this function:

void ShowTime

(

Timer

arg

)

Writes on screen the estimation of total and remaining CPU times. These informations are stored in the timer arg.

Parameters

arg

Argument

Returns

void

{
//  double ftime;                                   // real time
    double ctime;                                   // CPU time
//  unsigned int    ttime, rtime;   // total and remaining  real time (in seconds)
    unsigned int  tctime, rctime; // total and remaining CPU time (in seconds)

    int day, hour, min, sec;
    unsigned int rest;

    // --- Write total and remaining estimated time -----------------------

//  ftime  = arg.GetRealTime();
    ctime  = arg.CPUTime();
//  ttime  = (unsigned int)((ftime*IterationNb)/IterationNo);
//  rtime  = (unsigned int)((ftime*(IterationNb-IterationNo))/IterationNo);
    tctime = (unsigned int)((ctime*IterationNb)/IterationNo);
    rctime = (unsigned int)((ctime*(IterationNb-IterationNo))/IterationNo);

    // --- Show total time ------------------------------------------------

    rest =  tctime;
  day  =  rest/86400;
    rest %= 86400;
    hour =  rest/3600;
    rest %= 3600;
    min  =  rest/60;
    rest %= 60;
    sec  =  rest;

    printf("\033[1A\033[1A");

    if (tctime >= 86400)
        printf("Total CPU time     (estimation) : %5d day %2d h %2d min %2d s\n", day, hour, min, sec);

    if ((tctime < 86400)&&(tctime >= 3600))
        printf("Total CPU time     (estimation) : %2d h %2d min %2d s        \n", hour, min, sec);

    if ((tctime < 3600)&&(tctime >= 60))
        printf("Total CPU time     (estimation) : %2d min %2d s              \n", min, sec);

    if (tctime < 60)
        printf("Total CPU time     (estimation) : %2d s                      \n", sec);

    // --- Show remaining time ------------------------------------------------

    rest =  rctime;
  day  =  rest/86400;
    rest %= 86400;
    hour =  rest/3600;
    rest %= 3600;
    min  =  rest/60;
    rest %= 60;
    sec  =  rest;

    if (rctime >= 86400)
        printf("Remaining CPU time (estimation) : %5d day %2d h %2d min %2d s\n", day, hour, min, sec);

    if ((rctime < 86400)&&(rctime >= 3600))
        printf("Remaining CPU time (estimation) : %2d h %2d min %2d s        \n", hour, min, sec);

    if ((rctime < 3600)&&(rctime >= 60))
        printf("Remaining CPU time (estimation) : %2d min %2d s              \n", min, sec);

    if (rctime < 60)
        printf("Remaining CPU time (estimation) : %2d s                      \n", sec);

}

Here is the caller graph for this function:

int Sign

(

const real

a

)

Returns 1 if a is non-negative, -1 elsewhere.

Parameters

a	Real value

Returns

int

{
    if (a >= 0)
        return 1;
    else
        return -1;
}

Vector Source

(

Cell &

UserCell

)

Returns the source term in the cell UserCell.

Parameters

UserCell

Cell value

Returns

Vector

Gravity vector

{
    // ---  Local variables ---

    Vector Force(Dimension);
    Vector Result(QuantityNb);
    Result.setZero();

    Vector V(3);
    real Gx=0., Gy=0., Gz=0.,rho=0.;
    for(int i=1;i<=3;i++)
        V.setValue(i,UserCell.average(i+1));
    rho = UserCell.density();
    Gz = 0.2;
    Result.setValue(2,rho*Gx);
    Result.setValue(3,rho*Gy);
    Result.setValue(4,rho*Gz);
    Result.setValue(5,rho*(Gx*V.value(1) + Gy*V.value(2) + Gz*V.value(3)));

    Result.setZero();
    return Result;

}

Here is the caller graph for this function:

Matrix stateUstar	(	const Vector &	AvgL,
		const Vector &	AvgR,
		const real	prel,
		const real	prer,
		real &	slopeLeft,
		real &	slopeRight,
		real &	slopeM,
		real &	slopeLeftStar,
		real &	slopeRightStar,
		int	AxisNo
	)

Returns the intermediary states of HLLD numerical flux for MHD equations.

Parameters

AvgL	Left average vector
AvgR	Right average vector
prel	Left pressure
prer	Right pressure
slopeLeft	Left slope
slopeRight	Right slope
slopeM	Slope value computed for HLLD
slopeLeftStar	Left star slope (intermediary states)
slopeRightStar	Right star slope (intermediary states)
AxisNo	Axis of interest

Returns

Matrix

variables U-star

variables U-star-star

{
    real rhol=0.,rhor=0.;
    real psil=0.,psir=0.;
    real vxl=0.,vxr=0.,vyl=0.,vyr=0.,vzl=0.,vzr=0.;
    real Bxr=0.,Bxl=0.,Byl=0.,Byr=0.,Bzl=0.,Bzr=0.;
    real vl=0.,vr=0.,mf=0.;
    real Bl=0.,Br=0.;
    real pTl=0.,pTr=0.;
    real vBl=0.,vBr=0.;
    real el=0., er=0.;
    real rhols=0.,rhors=0.;
    real vxls=0.,vxrs=0.,vyls=0.,vyrs=0.,vzls=0.,vzrs=0.;
    real Bxls=0.,Bxrs=0.,Byls=0.,Byrs=0.,Bzls=0.,Bzrs=0.;
    real pTs=0.;
    real vBls=0.,vBrs=0.;
    real Els=0.,Ers=0.;
    real rholss=0.,rhorss=0.;
    real vxlss=0.,vxrss=0.,vylss=0.,vyrss=0.,vzlss=0.,vzrss=0.;
    real Bxlss=0.,Bxrss=0.,Bylss=0.,Byrss=0.,Bzlss=0.,Bzrss=0.;
    real vBlss=0.,vBrss=0.;
    real Elss=0.,Erss=0.;
    int Bsign=0;
    real half=0.5;
    real epsilon2=1e-16;
    real auxl=0.,auxr=0., Sqrtrhols=0.;

    Matrix U(QuantityNb,4);

    //Variables
    rhol = AvgL.value(1);
    vxl  = AvgL.value(2)/rhol;
    vyl  = AvgL.value(3)/rhol;
    vzl  = AvgL.value(4)/rhol;
    el   = AvgL.value(5);
    psil = AvgL.value(6);
    Bxl  = AvgL.value(7);
    Byl  = AvgL.value(8);
    Bzl  = AvgL.value(9);

    rhor = AvgR.value(1);
    vxr  = AvgR.value(2)/rhor;
    vyr  = AvgR.value(3)/rhor;
    vzr  = AvgR.value(4)/rhor;
    er   = AvgR.value(5);
    psir = AvgR.value(6);
    Bxr  = AvgR.value(7);
    Byr  = AvgR.value(8);
    Bzr  = AvgR.value(9);

    //v_x and v_y
    vl = AvgL.value(AxisNo+1)/rhol;
    vr = AvgR.value(AxisNo+1)/rhor;

    // Average B value
    mf = (AvgL.value(AxisNo+6)+AvgR.value(AxisNo+6))/(double (2.0));

        if(AxisNo==1)
        {
            //|B|
            Bl = sqrt( mf*mf + Byl*Byl + Bzl*Bzl );
            Br = sqrt( mf*mf + Byr*Byr + Bzr*Bzr );
                //Inner product v.B
            vBl = vxl*mf + vyl*Byl + vzl*Bzl;
            vBr = vxr*mf + vyr*Byr + vzr*Bzr;
        }else if(AxisNo==2){
            //|B|
            Bl = sqrt( Bxl*Bxl + mf*mf + Bzl*Bzl );
            Br = sqrt( Bxr*Bxr + mf*mf + Bzr*Bzr );
                //Inner product v.B
            vBl = vxl*Bxl + vyl*mf + vzl*Bzl;
            vBr = vxr*Bxr + vyr*mf + vzr*Bzr;
        }else{
            //|B|
            Bl = sqrt( Bxl*Bxl + Byl*Byl + mf*mf );
            Br = sqrt( Bxr*Bxr + Byr*Byr + mf*mf );
                //Inner product v.B
            vBl = vxl*Bxl + vyl*Byl + vzl*mf;
            vBr = vxr*Bxr + vyr*Byr + vzr*mf;
        }

    //Total Pressure
    pTl = prel + half*Bl*Bl;
    pTr = prer + half*Br*Br;

    // S_M - Equation 38
    slopeM = ((slopeRight - vr)*rhor*vr - (slopeLeft - vl)*rhol*vl - pTr +pTl)/
             ((slopeRight - vr)*rhor    - (slopeLeft - vl)*rhol);

    //Sign function
    if(mf > 0)  Bsign =  1;
    else        Bsign = -1;


            //density - Equation 43
            rhols = rhol*(slopeLeft-vl)/(slopeLeft-slopeM);
            rhors = rhor*(slopeRight-vr)/(slopeRight-slopeM);

            if(AxisNo==1){
                //velocities
                //Equation 39
                vxls = slopeM;
                vxrs = vxls;
                //B_x
                Bxls = mf;
                Bxrs = Bxls;

                auxl= (rhol*(slopeLeft  - vl)*(slopeLeft  - slopeM)-mf*mf);
                auxr= (rhor*(slopeRight - vr)*(slopeRight - slopeM)-mf*mf);

                    if( fabs(auxl)<epsilon2 ||
                    ( ( ( fabs(slopeM -vl) ) < epsilon2 )   &&
                    ( ( fabs(Byl) + fabs(Bzl) ) < epsilon2 )  &&
                    ( ( mf*mf ) >  (Gamma*prel) ) ) ) {
                        vyls = vyl;
                        vzls = vzl;
                        Byls = Byl;
                        Bzls = Bzl;
                    }else{
                        //Equation 44
                        vyls = vyl - mf*Byl*(slopeM-vl)/auxl;
                        //Equation 46
                        vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
                        //Equation 45
                        Byls = Byl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
                        //Equation 47
                        Bzls = Bzl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
                    }
                    if( fabs(auxr)<epsilon2 ||
                    ( ( ( fabs(slopeM -vr) ) < epsilon2 )   &&
                    ( ( fabs(Byr) + fabs(Bzr) ) < epsilon2 )  &&
                    ( ( mf*mf ) >  (Gamma*prer) ) ) ) {
                        vyrs = vyr;
                        vzrs = vzr;
                        Byrs = Byr;
                        Bzrs = Bzr;
                    }else{
                        //Equation 44
                        vyrs = vyr - mf*Byr*(slopeM-vr)/auxr;
                        //Equation 46
                        vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
                        //Equation 45;
                        Byrs = Byr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
                        //Equation 47
                        Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
                    }

            }else if(AxisNo==2){
                //velocities
                //Equation 39
                vyls = slopeM;
                vyrs = vyls;
                //B_y
                Byls = mf;
                Byrs = Byls;

                auxl= (rhol*(slopeLeft  - vl)*(slopeLeft  - slopeM)-mf*mf);
                auxr= (rhor*(slopeRight - vr)*(slopeRight - slopeM)-mf*mf);

                if( fabs(auxl)<epsilon2 ||
                    ( ( ( fabs(slopeM -vl) ) < epsilon2 )   &&
                    ( ( fabs(Bxl) + fabs(Bzl) ) < epsilon2 )  &&
                    ( ( mf*mf ) >  (Gamma*prel) ) ) ) {
                    vxls = vxl;
                    vzls = vzl;
                    Bxls = Bxl;
                    Bzls = Bzl;
                }else{
                    //Equation 44
                    vxls = vxl - mf*Bxl*(slopeM-vl)/auxl;
                    //Equation 46
                    vzls = vzl - mf*Bzl*(slopeM-vl)/auxl;
                    //Equation 45
                    Bxls = Bxl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
                    //Equation 47
                    Bzls = Bzl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
                }

                if( fabs(auxr)<epsilon2 ||
                    ( ( ( fabs(slopeM -vr) ) < epsilon2 )   &&
                    ( ( fabs(Bxr) + fabs(Bzr) ) < epsilon2 )  &&
                    ( ( mf*mf ) >  (Gamma*prer) ) ) ) {
                    vxrs = vxr;
                    vzrs = vzr;
                    Bxrs = Bxr;
                    Bzrs = Bzr;
                }else{
                    //Equation 44
                    vxrs = vxr - mf*Bxr*(slopeM-vr)/auxr;
                    //Equation 46
                    vzrs = vzr - mf*Bzr*(slopeM-vr)/auxr;
                    //Equation 45
                    Bxrs = Bxr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
                    //Equation 47
                    Bzrs = Bzr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
                }
            }else{
                //velocities
                //Equation 39
                vzls = slopeM;
                vzrs = vzls;
                //B_z
                Bzls = mf;
                Bzrs = Bzls;

                auxl= (rhol*(slopeLeft  - vl)*(slopeLeft  - slopeM)-mf*mf);
                auxr= (rhor*(slopeRight - vr)*(slopeRight - slopeM)-mf*mf);

                if( fabs(auxl)<epsilon2 ||
                    ( ( ( fabs(slopeM -vl) ) < epsilon2 )   &&
                    ( ( fabs(Bxl) + fabs(Byl) ) < epsilon2 )  &&
                    ( ( mf*mf ) >  (Gamma*prel) ) ) ) {
                    vxls = vxl;
                    vyls = vyl;
                    Bxls = Bxl;
                    Byls = Byl;
                }else{
                    //Equation 44
                    vxls = vxl - mf*Bxl*(slopeM-vl)/auxl;
                    //Equation 46
                    vyls = vyl - mf*Byl*(slopeM-vl)/auxl;
                    //Equation 45
                    Bxls = Bxl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
                    //Equation 47
                    Byls = Byl*(rhol*(slopeLeft-vl)*(slopeLeft-vl) - mf*mf)/auxl;
                }

                if( fabs(auxr)<epsilon2 ||
                    ( ( ( fabs(slopeM -vr) ) < epsilon2 )   &&
                    ( ( fabs(Bxr) + fabs(Byr) ) < epsilon2 )  &&
                    ( ( mf*mf ) >  (Gamma*prer) ) ) ) {
                    vxrs = vxr;
                    vyrs = vyr;
                    Bxrs = Bxr;
                    Byrs = Byr;
                }else{
                    //Equation 44
                    vxrs = vxr - mf*Bxr*(slopeM-vr)/auxr;
                    //Equation 46
                    vyrs = vyr - mf*Byr*(slopeM-vr)/auxr;
                    //Equation 45
                    Bxrs = Bxr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
                    //Equation 47
                    Byrs = Byr*(rhor*(slopeRight-vr)*(slopeRight-vr) - mf*mf)/auxr;
                }
            }

            //total pressure - Equation 41
            pTs = pTl + rhol*(slopeLeft - vl)*(slopeM - vl);

            //inner product vs*Bs
            vBls = vxls*Bxls + vyls*Byls + vzls*Bzls;
            vBrs = vxrs*Bxrs + vyrs*Byrs + vzrs*Bzrs;

            //energy - Equation 48
            Els = ((slopeLeft -vl)*el - pTl*vl+pTs*slopeM+mf*(vBl - vBls))/(slopeLeft  - slopeM);
            Ers = ((slopeRight-vr)*er - pTr*vr+pTs*slopeM+mf*(vBr - vBrs))/(slopeRight - slopeM);

            //U-star variables
            //Left
            U.setValue(1,1,rhols);
            U.setValue(2,1,rhols*vxls);
            U.setValue(3,1,rhols*vyls);
            U.setValue(4,1,rhols*vzls);
            U.setValue(5,1,Els);
            U.setValue(6,1,psil);
            U.setValue(7,1,Bxls);
            U.setValue(8,1,Byls);
            U.setValue(9,1,Bzls);
            //Right
            U.setValue(1,2,rhors);
            U.setValue(2,2,rhors*vxrs);
            U.setValue(3,2,rhors*vyrs);
            U.setValue(4,2,rhors*vzrs);
            U.setValue(5,2,Ers);
            U.setValue(6,2,psir);
            U.setValue(7,2,Bxrs);
            U.setValue(8,2,Byrs);
            U.setValue(9,2,Bzrs);




            Sqrtrhols = (sqrt(rhols)+sqrt(rhors));

            if((mf*mf/min((Bl*Bl),(Br*Br))) < epsilon2){
            for(int k=1;k<=QuantityNb;k++){
                    U.setValue(k,3, U.value(k,1));
                    U.setValue(k,4, U.value(k,2));
                }
            }else{
                //density - Equation 49
                rholss = rhols;
                rhorss = rhors;

                if(AxisNo==1){
                    //velocities
                    //Equation 39
                    vxlss = slopeM;
                    vxrss = slopeM;
                    //Equation 59
                    vylss = (sqrt(rhols)*vyls + sqrt(rhors)*vyrs + (Byrs - Byls)*Bsign)/Sqrtrhols;
                    vyrss = vylss;
                    //Equation 60
                    vzlss = (sqrt(rhols)*vzls + sqrt(rhors)*vzrs + (Bzrs - Bzls)*Bsign)/Sqrtrhols;
                    vzrss = vzlss;

                    //Magnetic Field components
                    //B_x
                    Bxlss = mf;
                    Bxrss = mf;
                    //Equation 61
                    Bylss = (sqrt(rhols)*Byrs + sqrt(rhors)*Byls + sqrt(rhols*rhors)*(vyrs - vyls)*Bsign)/Sqrtrhols;
                    Byrss = Bylss;
                    //Equation 62
                    Bzlss = (sqrt(rhols)*Bzrs + sqrt(rhors)*Bzls + sqrt(rhols*rhors)*(vzrs - vzls)*Bsign)/Sqrtrhols;
                    Bzrss = Bzlss;
                }
                else if(AxisNo==2){
                    //velocities
                    //Equation 59
                    vxlss = (sqrt(rhols)*vxls + sqrt(rhors)*vxrs + (Bxrs - Bxls)*Bsign)/Sqrtrhols;
                    vxrss = vxlss;
                    //Equation 39
                    vylss = slopeM;
                    vyrss = slopeM;
                    //Equation 60
                    vzlss = (sqrt(rhols)*vzls + sqrt(rhors)*vzrs + (Bzrs - Bzls)*Bsign)/Sqrtrhols;
                    vzrss = vzlss;

                    //Magnetic Field components
                    //Equation 61
                    Bxlss = (sqrt(rhols)*Bxrs + sqrt(rhors)*Bxls + sqrt(rhols*rhors)*(vxrs - vxls)*Bsign)/Sqrtrhols;
                    Bxrss = Bxlss;
                    //B_y
                    Bylss = mf;
                    Byrss = mf;
                    //Equation 62
                    Bzlss = (sqrt(rhols)*Bzrs + sqrt(rhors)*Bzls + sqrt(rhols*rhors)*(vzrs - vzls)*Bsign)/Sqrtrhols;
                    Bzrss = Bzlss;

                }else{
                    //velocities
                    //Equation 59
                    vxlss = (sqrt(rhols)*vxls + sqrt(rhors)*vxrs + (Bxrs - Bxls)*Bsign)/Sqrtrhols;
                    vxrss = vxlss;
                    //Equation 60
                    vylss = (sqrt(rhols)*vyls + sqrt(rhors)*vyrs + (Byrs - Byls)*Bsign)/Sqrtrhols;
                    vyrss = vylss;
                    //Equation 39
                    vzlss = slopeM;
                    vzrss = slopeM;

                    //Magnetic Field components
                    //Equation 61
                    Bxlss = (sqrt(rhols)*Bxrs + sqrt(rhors)*Bxls + sqrt(rhols*rhors)*(vxrs - vxls)*Bsign)/Sqrtrhols;
                    Bxrss = Bxlss;
                    //Equation 62
                    Bylss = (sqrt(rhols)*Byrs + sqrt(rhors)*Byls + sqrt(rhols*rhors)*(vyrs - vyls)*Bsign)/Sqrtrhols;
                    Byrss = Bylss;
                    //B_y
                    Bzlss = mf;
                    Bzrss = mf;

                }


                //inner product vss*Bss
                vBlss = vxlss*Bxlss + vylss*Bylss + vzlss*Bzlss;
                vBrss = vxrss*Bxrss + vyrss*Byrss + vzrss*Bzrss;

                //Energy - Equation 63
                Elss = Els - sqrt(rhols)*(vBls - vBlss)*Bsign;
                Erss = Ers + sqrt(rhors)*(vBrs - vBrss)*Bsign;




                //U-star-star variables
                //Left
                U.setValue(1,3,rholss);
                U.setValue(2,3,rholss*vxlss);
                U.setValue(3,3,rholss*vylss);
                U.setValue(4,3,rholss*vzlss);
                U.setValue(5,3,Elss);
                U.setValue(6,3,psil);
                U.setValue(7,3,Bxlss);
                U.setValue(8,3,Bylss);
                U.setValue(9,3,Bzlss);
                //Right
                U.setValue(1,4,rhorss);
                U.setValue(2,4,rhorss*vxrss);
                U.setValue(3,4,rhorss*vyrss);
                U.setValue(4,4,rhorss*vzrss);
                U.setValue(5,4,Erss);
                U.setValue(6,4,psir);
                U.setValue(7,4,Bxrss);
                U.setValue(8,4,Byrss);
                U.setValue(9,4,Bzrss);
            }
            //Equation 61 - S-star left and S-star right
            slopeLeftStar  = slopeM - Abs(mf)/sqrt(rhols);
            slopeRightStar = slopeM + Abs(mf)/sqrt(rhors);

    return U;
}

Here is the caller graph for this function:

real Step

(

real

x

)

Returns a step (1 if x <0, 0 if x >0, 0.5 if x=0)

Parameters

x	Real value

Returns

double

{
    if (x < 0.)
        return 1.;
    else if (x > 0.)
        return 0.;
    else
        return .5;
}

void TimeEvolution

(

FineMesh *

Root

)

Computes a time evolution on the regular fine mesh Root. Only for finite volume computations.

Parameters

Root

Fine mesh

Returns

void

{

        // --- Store cell-average values into temporary ---
        Root->store();

        for (StepNo = 1; StepNo <= StepNb; StepNo++)
        {
            // --- Compute divergence for neighbour cells ---
            //The same conception with background computations, see upper...
            Root->computeDivergence(1);
            // --- Runge-Kutta step for neighbour cells ---
            Root->RungeKutta(1);
            // --- Divergence cleaning source term
            Root->computeCorrection(1);
            // --- Start inter-CPU exchanges ---
            CPUExchange(Root, SendQ);
            // --- Compute divergence for internal cells ---
            Root->computeDivergence(0);
            // --- Runge-Kutta step for internal cells ---
            Root->RungeKutta(0);
            // --- Divergence cleaning source term
            Root->computeCorrection(0);

#if defined PARMPI
            CommTimer.start();    //Communication Timer Start
            //Waiting while inter-CPU exchanges are finished
            if (MPIRecvType == 1)                     //for nonblocking recive...
                MPI_Waitall(4*Dimension,req,st);
            CommTimer.stop();
#endif

        }

        // --- Check stability ---
        Root->checkStability();

        // --- Compute integral values ---
        Root->computeIntegral();

        // --- Compute elapsed time and adapt time step ---

        if (!ComputeCPUTimeRef)
        {
            Eigenvalue = Max(EigenvalueX,Max(EigenvalueY,EigenvalueZ));
            ElapsedTime += TimeStep;
            if (!ConstantTimeStep) AdaptTimeStep();
            // --- Compute divergence-free correction constant
            //ch = CFL*SpaceStep/TimeStep;
            ch = Max(CFL*SpaceStep/TimeStep, Eigenvalue);
        }

        // --- Compute time-average values ---

        if (TimeAveraging)
            Root->computeTimeAverage();

}

Here is the caller graph for this function:

void TimeEvolution

(

Node *

Root

)

Computes a time evolution on the tree structure, the root node being Root. Only for multiresolution computations.

Parameters

Root

Root node

Returns

void

{

        // --- Smooth data ---

        if ( SmoothCoeff != 0.)
            Root->smooth();

        // --- Store cell-average values of leaves ---
        Root->store();

        // --- Refresh tree structure ---
        RefreshTree(Root);

        for (StepNo = 1; StepNo <= StepNb; StepNo++)
        {
            // --- Compute divergence ---
            Root->computeDivergence();
            // --- Runge-Kutta step ---
            Root->RungeKutta();
            // --- Divergence cleaning source-terms
            Root->computeCorrection();

        }

        // --- Refresh tree structure ---
        RefreshTree(Root);



        // --- Check stability ---
        Root->checkStability();

        // --- Compute integral values ---
        Root->computeIntegral();

        // --- Compute total number of cells and leaves ---

        TotalCellNb += CellNb;
        TotalLeafNb += LeafNb;
 //cout<<"eigen= "<<Eigenvalue<<endl;
        // --- Compute elapsed time and adapt time step ---
        Eigenvalue = Max(EigenvalueX,Max(EigenvalueY,EigenvalueZ));
        ElapsedTime += TimeStep;
        if (!ConstantTimeStep) AdaptTimeStep();

        // --- Compute divergence-free correction constant
        //ch = CFL*SpaceStep/TimeStep;
        ch = Max(CFL*SpaceStep/TimeStep, Eigenvalue);

}

void View	(	FineMesh *	Root,
		const char *	AverageFileName
	)

Writes the current cel–averages of the fine mesh Root into file AverageFileName. Only for finite volume computations.

Parameters

Root	Fine mesh
AverageFileName	File name

Returns

void

{
    char buf[256];
    int iaux;


  char CPUFileName[255];
#if defined PARMPI
    sprintf(CPUFileName,"%d_%d_%d_%s",coords[0],coords[1],coords[2],AverageFileName);
#else
    strcpy(CPUFileName, AverageFileName);
#endif

    // write header for graphic visualization
    Root->writeHeader(CPUFileName);

    // write cell-average values for graphic visualization
    Root->writeAverage(CPUFileName);

  // Compress data
    if (Dimension != 1)
    {
        if (ZipData)
        {
                sprintf(buf,"gzip %s",CPUFileName);
            iaux=system(buf);
        }
    }

        // --- Write time-average values into file ---

        if (TimeAveraging)
            Root->writeTimeAverage("TimeAverage.dat");

}

Here is the caller graph for this function:

void View	(	Node *	Root,
		const char *	TreeFileName,
		const char *	MeshFileName,
		const char *	AverageFileName
	)

Writes the data of the tree structure into files TreeFileName (tree structure), MeshFileName (mesh) and AverageFileName (cell-averages). The root node is Root. Only for multiresolution computations.

Parameters

Root	Root node
TreeFileName	Tree file name
MeshFileName	Mesh file name
AverageFileName	Average file name

Returns

void

{
    char buf[256];
    int iaux;

    // write tree (debugging only)
    if (debug) Root->writeTree(TreeFileName);

//  Root->computeCorrection();

    // write mesh for graphic visualisation
    if (Dimension != 1)
    {
        Root->writeHeader(MeshFileName);
        Root->writeAverage(MeshFileName);

        // Compress data (if parameter ZipData is true)
        if (ZipData)
        {
            sprintf(buf,"gzip %s",MeshFileName);
            iaux=system(buf);
    }
  }
    else
        Root->writeMesh(MeshFileName);


    // write cell-averages in multiresolution representation (1D) or on fine grid (2-3D)
    if (Dimension != 1)
    {
        Root->writeFineGrid(AverageFileName,ScaleNb+PrintMoreScales);

        // Compress data
        if (ZipData)
        {
            sprintf(buf,"gzip %s",AverageFileName);
            iaux=system(buf);
        }
    }
    else
    {
        Root->writeHeader(AverageFileName);
        Root->writeAverage(AverageFileName);
  }
}

void ViewEvery	(	FineMesh *	Root,
		int	arg
	)

Same as previous for a fine mesh Root. Only for finite volume computations.

Parameters

Root	Fine mesh
arg	argument

Returns

void

{
    char AverageName[256];      // File name for AverageNNN.dat
    char AverageFormat[256];    // File format for AverageNNN.dat

    sprintf(AverageFormat, "Average%s0%ii.vtk", "%", DigitNumber(IterationNb));
    sprintf(AverageName, AverageFormat, arg);

    View(Root, AverageName);
}

Here is the caller graph for this function:

void ViewEvery	(	Node *	Root,
		int	arg
	)

Writes into file the data of the tree structure at iteration arg. The output file names are AverageNNN.dat and MeshNNN.dat, NNN being the iteration in an accurate format. The root node is Root. Only for multiresolution computations.

Parameters

Root	Root node
arg	Argument

Returns

void

{
    char AverageName[256];      // File name for AverageNNN.dat
    char MeshName[256];       // File name for MeshNNN.dat
    char AverageFormat[256];    // File format for AverageNNN.dat
    char MeshFormat[256];           // File format for MeshNNN.dat

    sprintf(AverageFormat, "Average%s0%ii.vtk", "%", DigitNumber(IterationNb));
    sprintf(AverageName, AverageFormat, arg);
    sprintf(MeshFormat, "Mesh%s0%ii.dat", "%", DigitNumber(IterationNb));
    sprintf(MeshName, MeshFormat,arg);

    View(Root, "Tree.dat", MeshName, AverageName);
}

void ViewIteration

(

FineMesh *

Root

)

Same as previous for a fine mesh Root. Only for finite volume computations.

Parameters

Root

Fine mesh

Returns

void

{
        if (IterationNo == PrintIt1)
                View(Root, "Average_1.vtk");

        if (IterationNo == PrintIt2)
                View(Root, "Average_2.vtk");

        if (IterationNo == PrintIt3)
                View(Root, "Average_3.vtk");

        if (IterationNo == PrintIt4)
                View(Root, "Average_4.vtk");

        if (IterationNo == PrintIt5)
                View(Root, "Average_5.vtk");

        if (IterationNo == PrintIt6)
                View(Root, "Average_6.vtk");
}

Here is the caller graph for this function:

void ViewIteration

(

Node *

Root

)

Writes into file the data of the tree structure from physical time PrintTime1 to physical time PrintTime6. The output file names are Average_N.dat and Mesh_N.dat, N being between 1 and 6. The root node is Root. Only for multiresolution computations.

Parameters

Root

Root node

Returns

void

{
        if (IterationNo == PrintIt1)
                View(Root, "Tree.dat", "Mesh_1.dat", "Average_1.vtk");

        if (IterationNo == PrintIt2)
                View(Root, "Tree.dat", "Mesh_2.dat", "Average_2.vtk");

        if (IterationNo == PrintIt3)
                View(Root, "Tree.dat", "Mesh_3.dat", "Average_3.vtk");

        if (IterationNo == PrintIt4)
                View(Root, "Tree.dat", "Mesh_4.dat", "Average_4.vtkt");

        if (IterationNo == PrintIt5)
                View(Root, "Tree.dat", "Mesh_5.dat", "Average_5.vtk");

        if (IterationNo == PrintIt6)
                View(Root, "Tree.dat", "Mesh_6.dat", "Average_6.vtk");
}