FPPowerSystem1BusInitialize.c File Reference

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <physicalmodels/fppowersystem1bus.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
double	FPPowerSystem1BusComputeCFL (void *, void *, double, double)
double	FPPowerSystem1BusComputeDiffNumber (void *, void *, double, double)
int	FPPowerSystem1BusAdvection (double *, double *, int, void *, double)
int	FPPowerSystem1BusDiffusionLaplacian (double *, double *, int, void *, double)
int	FPPowerSystem1BusDiffusionGeneral (double *, double *, int, int, void *, double)
int	FPPowerSystem1BusUpwind (double *, double *, double *, double *, double *, double *, int, void *, double)
int	FPPowerSystem1BusPostStep (double *, void *, void *, double, int)
int	FPPowerSystem1BusPrintStep (void *, void *, double)
int	FPPowerSystem1BusInitialize (void *s, void *m)

Function Documentation

double FPPowerSystem1BusComputeCFL	(	void *	,
		void *	,
		double	,
		double
	)

Definition at line 11 of file FPPowerSystem1BusComputeCFL.c.

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

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

  double  max_cfl = 0;
  int     index[ndims];
  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    double x;     _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->x,x);
    double y;     _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->x,y);
    double dxinv; _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
    double dyinv; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
    double drift_x= FPPowerSystem1BusDriftFunction(_XDIR_,params,x,y,t);
    double drift_y= FPPowerSystem1BusDriftFunction(_YDIR_,params,x,y,t);

    double local_cfl_x = absolute(drift_x) * dt * dxinv;
    double local_cfl_y = absolute(drift_y) * dt * dyinv;

    if (local_cfl_x > max_cfl) max_cfl = local_cfl_x;
    if (local_cfl_y > max_cfl) max_cfl = local_cfl_y;

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

  return(max_cfl);
}

double FPPowerSystem1BusComputeDiffNumber	(	void *	,
		void *	,
		double	,
		double
	)

Definition at line 11 of file FPPowerSystem1BusComputeDiffNumber.c.

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

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

  double  max_diff = 0;
  int     index[ndims];
  int     done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    double dxinv; _GetCoordinate_(0,index[0],dim,ghosts,solver->dxinv,dxinv);
    double dyinv; _GetCoordinate_(1,index[1],dim,ghosts,solver->dxinv,dyinv);
    double dissp_yx= FPPowerSystem1BusDissipationFunction(_YDIR_,_XDIR_,params,t);
    double dissp_yy= FPPowerSystem1BusDissipationFunction(_YDIR_,_YDIR_,params,t);

    double local_diff_yx = absolute(dissp_yx) * dt * dyinv * dxinv;
    double local_diff_yy = absolute(dissp_yy) * dt * dyinv * dyinv;

    if (local_diff_yx > max_diff) max_diff = local_diff_yx;
    if (local_diff_yy > max_diff) max_diff = local_diff_yy;

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

  return(max_diff);
}

int FPPowerSystem1BusAdvection	(	double *	,
		double *	,
		int	,
		void *	,
		double
	)

Definition at line 7 of file FPPowerSystem1BusAdvection.c.

{
  HyPar *solver = (HyPar*) s;
  int   *dim    = solver->dim_local;
  int   ghosts  = solver->ghosts;
  int   ndims   = solver->ndims;
  int   nvars   = solver->nvars;

  /* calculate total size of arrays */
  int bounds[ndims]; _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  int size;          _ArrayProduct1D_(bounds,ndims,size); size *= nvars;

  /* f = u */
  _ArrayCopy1D_(u,f,size);

  return(0);
}

int FPPowerSystem1BusDiffusionLaplacian	(	double *	,
		double *	,
		int	,
		void *	,
		double
	)

Definition at line 9 of file FPPowerSystem1BusDiffusion.c.

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

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

  /* calculate total size of arrays */
  int bounds[ndims]; _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  int size;          _ArrayProduct1D_(bounds,ndims,size); size *= nvars;

  /* calculate dissipation coefficient  -- constant in x and y */
  double dissipation = FPPowerSystem1BusDissipationFunction(dir,dir,params,t);

  /* f = dissipation * u */
  _ArrayScaleCopy1D_(u,dissipation,f,size);

  return(0);
}

int FPPowerSystem1BusDiffusionGeneral	(	double *	,
		double *	,
		int	,
		int	,
		void *	,
		double
	)

Definition at line 32 of file FPPowerSystem1BusDiffusion.c.

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

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

  /* calculate total size of arrays */
  int bounds[ndims]; _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  int size;          _ArrayProduct1D_(bounds,ndims,size); size *= nvars;

  /* calculate dissipation coefficient  -- constant in x and y */
  double dissipation = FPPowerSystem1BusDissipationFunction(dir1,dir2,params,t);

  /* f = dissipation * u */
  _ArrayScaleCopy1D_(u,dissipation,f,size);

  return(0);
}

int FPPowerSystem1BusUpwind	(	double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		double *	,
		int	,
		void *	,
		double
	)

Definition at line 9 of file FPPowerSystem1BusUpwind.c.

{
  HyPar             *solver = (HyPar*)              s;
  FPPowerSystem1Bus *params = (FPPowerSystem1Bus*)  solver->physics;
  int               done,v;

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

  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);
      double x = 0,y = 0; /* x,y coordinates of the interface */
      if (dir == 0) {
        /* x-interface */
        double x1, x2;
        _GetCoordinate_(_XDIR_,index_inter[_XDIR_]-1,dim,ghosts,solver->x,x1);
        _GetCoordinate_(_XDIR_,index_inter[_XDIR_]  ,dim,ghosts,solver->x,x2);
        x = 0.5 * ( x1 + x2 );
        _GetCoordinate_(_YDIR_,index_inter[_YDIR_],dim,ghosts,solver->x,y);
      } else if (dir == 1) {
        /* y-interface */
        _GetCoordinate_(_XDIR_,index_inter[_XDIR_],dim,ghosts,solver->x,x);
        double y1, y2;
        _GetCoordinate_(_YDIR_,index_inter[_YDIR_]-1,dim,ghosts,solver->x,y1);
        _GetCoordinate_(_YDIR_,index_inter[_YDIR_]  ,dim,ghosts,solver->x,y2);
        y = 0.5 * ( y1 + y2 );
      }
      double drift = FPPowerSystem1BusDriftFunction(dir,params,x,y,t);
      for (v = 0; v < nvars; v++)  
        fI[nvars*p+v] = drift * (drift > 0 ? fL[nvars*p+v] : fR[nvars*p+v] );
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

int FPPowerSystem1BusPostStep	(	double *	,
		void *	,
		void *	,
		double	,
		int
	)

Definition at line 8 of file FPPowerSystem1BusPostStep.c.

{
  HyPar             *solver = (HyPar*)              s;
  MPIVariables      *mpi    = (MPIVariables*)       m;
  FPPowerSystem1Bus *params = (FPPowerSystem1Bus*)  solver->physics;
  _DECLARE_IERR_;

  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int ndims   = solver->ndims;
  int index[ndims];

  double local_sum = 0;
  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
    double dxinv; _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
    double dyinv; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
    local_sum     += (u[p] / (dxinv * dyinv));
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }
  double local_integral = local_sum;
  double global_integral = 0; 
  IERR MPISum_double(&global_integral,&local_integral,1,&mpi->world); CHECKERR(ierr);
  params->pdf_integral = global_integral;

  return(0);
}

int FPPowerSystem1BusPrintStep	(	void *	,
		void *	,
		double
	)

Definition at line 5 of file FPPowerSystem1BusPrintStep.c.

{
  HyPar             *solver = (HyPar*)              s;
  FPPowerSystem1Bus *params = (FPPowerSystem1Bus*)  solver->physics;
  printf("Domain integral of the probability density function: %1.16E\n",params->pdf_integral);
  return(0);
}

int FPPowerSystem1BusInitialize	(	void *	s,
		void *	m
	)

Definition at line 21 of file FPPowerSystem1BusInitialize.c.

{
  HyPar               *solver  = (HyPar*)             s;
  MPIVariables        *mpi     = (MPIVariables*)      m; 
  FPPowerSystem1Bus   *physics = (FPPowerSystem1Bus*) solver->physics;
  int                 ferr;
  _DECLARE_IERR_;

  static int count = 0;

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

  /* default values of model parameters */
  physics->omegaS = 1.0;
  physics->omegaB = 120 * (4.0 * atan(1.0));
  physics->H      = 5.0;
  physics->D      = 5.0;
  physics->Pm_avg = 0.9;
  physics->Pmax   = 2.1;
  physics->sigma  = 0.2;
  physics->lambda = 100.0 / physics->omegaB;

  /* reading physical model specific inputs */
  if (!mpi->rank) {
    FILE *in;
    if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
    in = fopen("physics.inp","r");
    if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
    else {
      char word[_MAX_STRING_SIZE_];
      ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
      if (!strcmp(word, "begin")){
          while (strcmp(word, "end")){
              ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
          if      (!strcmp(word, "omegaS")) {ferr=fscanf(in,"%lf",&physics->omegaS);if(ferr!=1)return(1);}
          else if (!strcmp(word, "omegaB")) {ferr=fscanf(in,"%lf",&physics->omegaB);if(ferr!=1)return(1);}
          else if (!strcmp(word, "H"     )) {ferr=fscanf(in,"%lf",&physics->H     );if(ferr!=1)return(1);}
          else if (!strcmp(word, "D"     )) {ferr=fscanf(in,"%lf",&physics->D     );if(ferr!=1)return(1);}
          else if (!strcmp(word, "Pm_avg")) {ferr=fscanf(in,"%lf",&physics->Pm_avg);if(ferr!=1)return(1);}
          else if (!strcmp(word, "Pmax"  )) {ferr=fscanf(in,"%lf",&physics->Pmax  );if(ferr!=1)return(1);}
          else if (!strcmp(word, "sigma" )) {ferr=fscanf(in,"%lf",&physics->sigma );if(ferr!=1)return(1);}
          else if (!strcmp(word, "lambda")) {ferr=fscanf(in,"%lf",&physics->lambda);if(ferr!=1)return(1);}
        }
        } else {
          fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
        return(1);
        }
    }
    fclose(in);
  }

#ifndef serial
  IERR MPIBroadcast_double(&physics->omegaS,1,0,&mpi->world);                        CHECKERR(ierr);
  IERR MPIBroadcast_double(&physics->omegaB,1,0,&mpi->world);                        CHECKERR(ierr);
  IERR MPIBroadcast_double(&physics->H     ,1,0,&mpi->world);                        CHECKERR(ierr);
  IERR MPIBroadcast_double(&physics->D     ,1,0,&mpi->world);                        CHECKERR(ierr);
  IERR MPIBroadcast_double(&physics->Pm_avg,1,0,&mpi->world);                        CHECKERR(ierr);
  IERR MPIBroadcast_double(&physics->Pmax  ,1,0,&mpi->world);                        CHECKERR(ierr);
  IERR MPIBroadcast_double(&physics->sigma ,1,0,&mpi->world);                        CHECKERR(ierr);
  IERR MPIBroadcast_double(&physics->lambda,1,0,&mpi->world);                        CHECKERR(ierr);
#endif

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

  /* initializing physical model-specific functions */
  solver->ComputeCFL         = FPPowerSystem1BusComputeCFL;
  solver->ComputeDiffNumber  = FPPowerSystem1BusComputeDiffNumber;
  solver->FFunction          = FPPowerSystem1BusAdvection;
  solver->GFunction          = FPPowerSystem1BusDiffusionLaplacian;
  solver->HFunction          = FPPowerSystem1BusDiffusionGeneral;
  solver->Upwind             = FPPowerSystem1BusUpwind;
  solver->PostStep           = FPPowerSystem1BusPostStep;
  solver->PrintStep          = FPPowerSystem1BusPrintStep;

  /* check that solver is using the correct diffusion formulation */
/*
  if ((strcmp(solver->spatial_type_par,_NC_2STAGE_)) && (strcmp(solver->spatial_type_par,_NC_1_5STAGE_))) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in FPPowerSystem1BusInitialize(): Parabolic term spatial discretization must be ");
      fprintf(stderr,"\"%s\" or \"%s\".\n",_NC_2STAGE_,_NC_1_5STAGE_);
    }
    return(1);
  }
*/

  count++;
  return(0);
}