VlasovInitialize.c File Reference

Initialize the Vlasov module. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <bandedmatrix.h>
#include <physicalmodels/vlasov.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
double	VlasovComputeCFL (void *, void *, double, double)
int	VlasovAdvection (double *, double *, int, void *, double)
int	VlasovUpwind (double *, double *, double *, double *, double *, double *, int, void *, double)
int	VlasovWriteEFieldAndPotential (void *, void *, double)
int	VlasovPreStep (double *, void *, void *, double)
int	VlasovPostStage (double *, void *, void *, double)
int	VlasovEField (double *, void *, double)
int	VlasovInitialize (void *s, void *m)

Detailed Description

Initialize the Vlasov module.

Author

John Loffeld

Definition in file VlasovInitialize.c.

Function Documentation

double VlasovComputeCFL	(	void *	s,
		void *	m,
		double	dt,
		double	t
	)

Compute the maximum CFL number

Computes the maximum CFL number over the domain. Note that the CFL is computed over the local domain on this processor only.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
dt	Time step size for which to compute the CFL
t	Time

Definition at line 19 of file VlasovComputeCFL.c.

{
  HyPar   *solver = (HyPar*)  s;
  Vlasov  *params = (Vlasov*) solver->physics;

  int     ndims  = solver->ndims;
  int     ghosts = solver->ghosts;
  int     *dim   = solver->dim_local;
  double  *u     = solver->u;

  double max_cfl = 0;
  int done = 0; 
  int index[ndims];
  _ArraySetValue_(index,ndims,0);

  while (!done) {

    for (int dir=0; dir<ndims; dir++) {
      double dxinv;
      _GetCoordinate_(dir,index[dir],dim,ghosts,solver->dxinv,dxinv);
      double eig = VlasovAdvectionCoeff(index, dir, solver);
      double local_cfl = eig*dt*dxinv; 
      if (local_cfl > max_cfl) max_cfl = local_cfl;
    }

    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(max_cfl);
}

int VlasovAdvection	(	double *	f,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the advection term

Compute the hyperbolic flux over the local domain.

\begin{equation} {\bf F}\left({\bf u}(x,v)\right) = c {\bf u}(x,v), \end{equation}

where the advection coefficient \(c\) is computed by VlasovAdvectionCoeff().

Parameters

f	Array to hold the computed flux (same size and layout as u)
u	Array containing the conserved solution
dir	Spatial dimension
s	Solver object of type HyPar
t	Current time

Definition at line 22 of file VlasovAdvection.c.

{

  HyPar  *solver = (HyPar*)  s;
  Vlasov *param  = (Vlasov*) solver->physics;

  int* dim    = solver->dim_local;
  int  ghosts = solver->ghosts;
  int  ndims  = solver->ndims;

  // set bounds for array index to include ghost points
  int bounds[ndims];
  _ArrayCopy1D_(dim,bounds,ndims);
  for (int i=0; i<ndims; i++) bounds[i] += 2*ghosts;

  // set offset such that index is compatible with ghost point arrangement
  int offset[ndims];
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; 
  int index_wog[ndims];
  int index[ndims]; _ArraySetValue_(index,ndims,0);
  while (!done) {

    _ArrayCopy1D_(index, index_wog, ndims);
    for (int i = 0; i < ndims; i++) index_wog[i] -= ghosts;
    double adv_coeff = VlasovAdvectionCoeff(index_wog, dir, solver);

    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    f[p] = adv_coeff * u[p];
    _ArrayIncrementIndex_(ndims,bounds,index,done);

  }

  return 0;
}

int VlasovUpwind	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the upwind flux at interfaces

Upwinding scheme for the Vlasov equations

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension
s	Solver object of type HyPar
t	Current solution time

Definition at line 17 of file VlasovUpwind.c.

{
  HyPar  *solver = (HyPar*)  s;
  Vlasov *param  = (Vlasov*) solver->physics;
  int    done;

  int ndims   = solver->ndims;
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;

  int index_outer[ndims], 
      index_inter[ndims], 
      bounds_outer[ndims], 
      bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {

      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int idxL[ndims], idxR[ndims];

      int neig = 2+4*(ndims-1);
      double eig[neig];
      int count = 0;
      eig[count] = VlasovAdvectionCoeff(indexL, dir, solver); count++;
      eig[count] = VlasovAdvectionCoeff(indexR, dir, solver); count++;
      for (int tdir = 0; tdir < ndims; tdir++ ) {
        if (tdir != dir) {
          _ArrayCopy1D_(indexL, idxL, ndims); idxL[tdir]--;
          _ArrayCopy1D_(indexR, idxR, ndims); idxR[tdir]--;
          eig[count] = VlasovAdvectionCoeff(idxL, dir, solver); count++;
          eig[count] = VlasovAdvectionCoeff(idxR, dir, solver); count++;
          _ArrayCopy1D_(indexL, idxL, ndims); idxL[tdir]++;
          _ArrayCopy1D_(indexR, idxR, ndims); idxR[tdir]++;
          eig[count] = VlasovAdvectionCoeff(idxL, dir, solver); count++;
          eig[count] = VlasovAdvectionCoeff(idxR, dir, solver); count++;
        }
      }
      if (count != neig) {
        fprintf(stderr, "Error in VlasovUpwind(): count != neig for dir=%d\n",dir);
        return 1;
      }

      int all_positive = 1, all_negative = 1;
      double maxabs_eig = 0.0;
      for (int n = 0; n<neig; n++) {
        if (eig[n] <= 0) all_positive = 0;
        if (eig[n] >= 0) all_negative = 0;
        if (absolute(eig[n]) > maxabs_eig) {
          maxabs_eig = absolute(eig[n]);
        }
      }

      if (all_positive) {
        fI[p] = fL[p];
      } else if (all_negative) {
        fI[p] = fR[p];
      } else { 
        fI[p] = 0.5 * (fL[p] + fR[p] - maxabs_eig * (uR[p] - uL[p]));
      }

    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);

  }

  return 0;
}

int VlasovWriteEFieldAndPotential	(	void *	s,
		void *	m,
		double	a_t
	)

Write E-field and potential to file

Write out the electric field and potential to file

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
a_t	Current simulation time

Definition at line 23 of file VlasovWriteEFieldAndPotential.c.

{
  HyPar         *solver = (HyPar*)        s;
  MPIVariables  *mpi    = (MPIVariables*) m;
  Vlasov        *param  = (Vlasov*)       solver->physics;

  {
    char fname_root[_MAX_STRING_SIZE_] = "efield";
    VlasovWriteSpatialField( solver, mpi, param->e_field, fname_root );
    if (!strcmp(solver->plot_solution,"yes")) {
      VlasovPlotSpatialField( solver, mpi, param->e_field, a_t, fname_root );
    }
  }

  if (param->self_consistent_electric_field) {
    char fname_root[_MAX_STRING_SIZE_] = "potential";
    VlasovWriteSpatialField( solver, mpi, param->potential, fname_root );
    if (!strcmp(solver->plot_solution,"yes")) {
      VlasovPlotSpatialField( solver, mpi, param->potential, a_t, fname_root );
    }
  }

  return 0;
}

int VlasovPreStep	(	double *	u,
		void *	s,
		void *	m,
		double	waqt
	)

Vlasov-specific function called at the beginning of each time-step: Calls the function to set the electric field

Parameters

u	Solution (conserved variables)
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current solution time

Definition at line 14 of file VlasovFunctions.c.

{
  HyPar  *solver = (HyPar*) s;
  Vlasov *param  = (Vlasov*) solver->physics;

  int ierr = VlasovEField(u, solver, waqt);
  if (ierr) return ierr;

  return 0;
}

int VlasovPostStage	(	double *	u,
		void *	s,
		void *	m,
		double	waqt
	)

Vlasov-specific function called at the end of each stage in a multistage time integrator: Calls the function to set the electric field

Parameters

u	Solution (conserved variables)
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current solution time

Definition at line 33 of file VlasovFunctions.c.

{
  HyPar  *solver = (HyPar*) s;
  Vlasov *param  = (Vlasov*) solver->physics;

  int ierr = VlasovEField(u, solver, waqt);
  if (ierr) return ierr;

  return 0;
}

int VlasovEField	(	double *	u,
		void *	s,
		double	t
	)

Compute electric field.

Parameters

u	Conserved solution
s	Solver object of type HyPar
t	Current time

Definition at line 258 of file VlasovEField.c.

{
  HyPar  *solver = (HyPar*)  s;
  Vlasov *param  = (Vlasov*) solver->physics;

  int ierr;

  if (param->self_consistent_electric_field) {
    ierr = SetEFieldSelfConsistent(u, s, t);
  } else {
    ierr = SetEFieldPrescribed(u, s, t);
  }

  return ierr;
}

int VlasovInitialize	(	void *	s,
		void *	m
	)

Initialize the Vlasov physics module - allocate and set physics-related parameters, read physics-related inputs from file, and set the physics-related function pointers in HyPar

Parameters

s	Solver object of type HyPar
m	Object of type MPIVariables containing MPI-related info

Definition at line 43 of file VlasovInitialize.c.

{
  HyPar        *solver    = (HyPar*)        s;
  MPIVariables *mpi       = (MPIVariables*) m; 
  Vlasov       *physics   = (Vlasov*)       solver->physics;

  int *dim_global = solver->dim_global;
  int *dim_local = solver->dim_local;
  int ghosts = solver->ghosts;

  if (solver->nvars != _MODEL_NVARS_) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in VlasovInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
    }
    return(1);
  }
  if (solver->ndims != _MODEL_NDIMS_) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in VlasovInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
    }
    return(1);
  }

  /* default is prescribed electric field */
  physics->self_consistent_electric_field = 0;
  physics->ndims_x = 1;
  physics->ndims_v = 1;

  /* reading physical model specific inputs - all processes */
  if (!mpi->rank) {
    FILE *in;
    in = fopen("physics.inp","r");
    if (in) {
      printf("Reading physical model inputs from file \"physics.inp\".\n");
      char word[_MAX_STRING_SIZE_];
      int ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
      if (!strcmp(word, "begin")){
        while (strcmp(word, "end")){
          int ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
          if (!strcmp(word, "self_consistent_electric_field")) {
            /* read whether electric field is self-consistent or prescribed */
            int ferr = fscanf(in,"%d", &physics->self_consistent_electric_field);
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "x_ndims")) {
            /* read number of spatial dimensions */
            int ferr = fscanf(in,"%d", &physics->ndims_x);
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "v_ndims")) {
            /* read number of velocity dimensions */
            int ferr = fscanf(in,"%d", &physics->ndims_v);
            if (ferr != 1) return(1);
          } else if (strcmp(word,"end")) {
            char useless[_MAX_STRING_SIZE_];
            int ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
            printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",
                   word, useless);
            printf("recognized or extraneous. Ignoring.\n");
          }
        }
      } else {
        fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
        return(1);
      }
    }
    fclose(in);
  }

  if ((physics->ndims_x+physics->ndims_v) != solver->ndims) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in VlasovInitialize:\n");
      fprintf(stderr, "  space + vel dims not equal to ndims!\n");
    }
    return(1);
  }

  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in VlasovInitialize: This physical model does not have a splitting ");
      fprintf(stderr,"of the hyperbolic term defined.\n");
    }
    return(1);
  }

#ifndef serial
  /* Broadcast parsed problem data */
  MPIBroadcast_integer(&physics->ndims_x,1,0,&mpi->world);
  MPIBroadcast_integer(&physics->ndims_v,1,0,&mpi->world);
  MPIBroadcast_integer((int *) &physics->self_consistent_electric_field,
                       1,0,&mpi->world);
#endif

  /* compute local number of x-space points with ghosts */
  physics->npts_local_x_wghosts = 1;
  physics->npts_local_x = 1;
  physics->npts_global_x_wghosts = 1;
  physics->npts_global_x = 1;
  for (int d=0; d<physics->ndims_x; d++) {
    physics->npts_local_x_wghosts *= (dim_local[d]+2*ghosts);
    physics->npts_local_x *= dim_local[d];
    physics->npts_global_x_wghosts *= (dim_global[d]+2*ghosts);
    physics->npts_global_x *= dim_global[d];
  }

  /* allocate array to hold the electric field (needs to have ghosts);
     note that number of electric field components is the number of
     spatial dimensions */
  physics->e_field = (double*) calloc(  physics->npts_local_x_wghosts
                                      * physics->ndims_x,
                                        sizeof(double)  );
  physics->potential = (double*) calloc(  physics->npts_local_x_wghosts
                                      * physics->ndims_x, 
                                        sizeof(double)  );
  
  /* Put the mpi object in the params for access in other functions */
  physics->m = m;
  
  if (physics->self_consistent_electric_field) {
#ifdef fftw
    /* If using FFTW, make sure MPI is enabled
       Currently we only support the distrubuted memory version
    */
#ifdef serial
    fprintf(stderr,"Error in VlasovInitialize(): Using FFTW requires MPI to be enabled.\n");
    return(1);
#else

    if (physics->ndims_x > 1) {
      fprintf(stderr,"Error in VlasovInitialize():\n");
      fprintf(stderr,"  Self-consistent electric field is implemented for only 1 space dimension.\n");
      return 1;
    }

    /* Create a scratch buffer for moving between real and complex values */
    physics->sum_buffer = (double*) calloc(dim_local[0], sizeof(double));
  
    /* Initialize FFTW and set up data buffers used for the transforms */
    fftw_mpi_init();
    physics->alloc_local = fftw_mpi_local_size_1d(dim_global[0], mpi->comm[0],
                                                  FFTW_FORWARD, 0,
                                                  &physics->local_ni,
                                                  &physics->local_i_start,
                                                  &physics->local_no,
                                                  &physics->local_o_start);
    if (dim_local[0] != physics->local_ni) {
      fprintf(stderr,"Error in VlasovInitialize(): The FFTW data distribution is incompatible with the HyPar one.\n");
      fprintf(stderr,"Decompose the spatial dimension so that the degrees of freedom are evenly divided.\n");
      return(1);
    }
  
    physics->phys_buffer_e = fftw_alloc_complex(physics->alloc_local);
    physics->fourier_buffer_e = fftw_alloc_complex(physics->alloc_local);

    physics->plan_forward_e = fftw_mpi_plan_dft_1d(dim_global[0],
                                                 physics->phys_buffer_e,
                                                 physics->fourier_buffer_e,
                                                 mpi->comm[0],
                                                 FFTW_FORWARD,
                                                 FFTW_ESTIMATE);

    physics->plan_backward_e = fftw_mpi_plan_dft_1d(dim_global[0],
                                                  physics->fourier_buffer_e,
                                                  physics->phys_buffer_e,
                                                  mpi->comm[0],
                                                  FFTW_BACKWARD,
                                                  FFTW_ESTIMATE);

    physics->phys_buffer_phi = fftw_alloc_complex(physics->alloc_local);
    physics->fourier_buffer_phi = fftw_alloc_complex(physics->alloc_local);

    physics->plan_forward_phi = fftw_mpi_plan_dft_1d(dim_global[0],
                                                 physics->phys_buffer_phi,
                                                 physics->fourier_buffer_phi,
                                                 mpi->comm[0],
                                                 FFTW_FORWARD,
                                                 FFTW_ESTIMATE);

    physics->plan_backward_phi = fftw_mpi_plan_dft_1d(dim_global[0],
                                                  physics->fourier_buffer_phi,
                                                  physics->phys_buffer_phi,
                                                  mpi->comm[0],
                                                  FFTW_BACKWARD,
                                                  FFTW_ESTIMATE);
#endif
#else
  fprintf(stderr,"Error in VlasovInitialize():\n");
  fprintf(stderr,"  Self-consistent electric field requires FFTW library.\n");
  return(1);
#endif
  }

  /* initializing physical model-specific functions */
  solver->PreStep       = VlasovPreStep;
  solver->ComputeCFL    = VlasovComputeCFL;
  solver->FFunction     = VlasovAdvection;
  solver->Upwind        = VlasovUpwind;
  solver->PhysicsOutput = VlasovWriteEFieldAndPotential;
  solver->PostStage     = VlasovPostStage;

  int ierr = VlasovEField(solver->u, solver, 0.0);
  if (ierr) return ierr;

  return 0;
}