FPPowerSystem3BusInitialize.c File Reference

Function to initialize the 3-bus power system model. More...

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

Go to the source code of this file.

Functions
double	FPPowerSystem3BusComputeCFL (void *, void *, double, double)
double	FPPowerSystem3BusComputeDiffNumber (void *, void *, double, double)
int	FPPowerSystem3BusAdvection (double *, double *, int, void *, double)
int	FPPowerSystem3BusDiffusion (double *, double *, int, int, void *, double)
int	FPPowerSystem3BusUpwind (double *, double *, double *, double *, double *, double *, int, void *, double)
int	FPPowerSystem3BusInitialize (void *s, void *m)

Detailed Description

Function to initialize the 3-bus power system model.

Author

Debojyoti Ghosh

Definition in file FPPowerSystem3BusInitialize.c.

Function Documentation

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

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 FPPowerSystem3BusComputeCFL.c.

{
  HyPar             *solver = (HyPar*)              s;
  FPPowerSystem3Bus *params = (FPPowerSystem3Bus*)  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[ndims], dxinv[ndims], drift[ndims];
    _GetCoordinate_(0,index[0],dim,ghosts,solver->x,x[0]);
    _GetCoordinate_(1,index[1],dim,ghosts,solver->x,x[1]);
    _GetCoordinate_(2,index[2],dim,ghosts,solver->x,x[2]);
    _GetCoordinate_(3,index[3],dim,ghosts,solver->x,x[3]);
    _GetCoordinate_(0,index[0],dim,ghosts,solver->dxinv,dxinv[0]);
    _GetCoordinate_(1,index[1],dim,ghosts,solver->dxinv,dxinv[1]);
    _GetCoordinate_(2,index[2],dim,ghosts,solver->dxinv,dxinv[2]);
    _GetCoordinate_(3,index[3],dim,ghosts,solver->dxinv,dxinv[3]);

    FPPowerSystem3BusDriftFunction(0,params,x,t,drift);

    double local_cfl[ndims];
    local_cfl[0] = absolute(drift[0]) * dt * dxinv[0];
    local_cfl[1] = absolute(drift[1]) * dt * dxinv[1];
    local_cfl[2] = absolute(drift[2]) * dt * dxinv[2];
    local_cfl[3] = absolute(drift[3]) * dt * dxinv[3];

    if (local_cfl[0] > max_cfl) max_cfl = local_cfl[0];
    if (local_cfl[1] > max_cfl) max_cfl = local_cfl[1];
    if (local_cfl[2] > max_cfl) max_cfl = local_cfl[2];
    if (local_cfl[3] > max_cfl) max_cfl = local_cfl[3];

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

  return(max_cfl);
}

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

Computes the maximum diffusion number over the domain. Note that the diffusion 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 FPPowerSystem3BusComputeDiffNumber.c.

{
  HyPar             *solver = (HyPar*)              s;
  FPPowerSystem3Bus *params = (FPPowerSystem3Bus*)  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[ndims],dissp[ndims*ndims];
    _GetCoordinate_(0,index[0],dim,ghosts,solver->dxinv,dxinv[0]);
    _GetCoordinate_(1,index[1],dim,ghosts,solver->dxinv,dxinv[1]);
    _GetCoordinate_(2,index[2],dim,ghosts,solver->dxinv,dxinv[2]);
    _GetCoordinate_(3,index[3],dim,ghosts,solver->dxinv,dxinv[3]);
    FPPowerSystem3BusDissipationFunction(0,0,params,t,dissp);

    double local_diff[ndims];
    local_diff[0] = absolute(dissp[0*ndims+0]) * dt * dxinv[0] * dxinv[0];
    local_diff[1] = absolute(dissp[1*ndims+1]) * dt * dxinv[1] * dxinv[1];
    local_diff[2] = absolute(dissp[2*ndims+2]) * dt * dxinv[2] * dxinv[2];
    local_diff[3] = absolute(dissp[3*ndims+3]) * dt * dxinv[3] * dxinv[3];

    if (local_diff[0] > max_diff) max_diff = local_diff[0];
    if (local_diff[1] > max_diff) max_diff = local_diff[1];
    if (local_diff[2] > max_diff) max_diff = local_diff[2];
    if (local_diff[3] > max_diff) max_diff = local_diff[3];

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

  return(max_diff);
}

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

Compute the advection term for the FPPowerSystem3Bus system: Since the advection coefficient is a function of x and not the solution, here the flux is set to the solution. The advection velocity is multiplied during upwinding FPPowerSystem3BusUpwind().

Parameters

f	Array to hold the computed flux vector (same layout as u)
u	Array with the solution vector
dir	Spatial dimension for which to compute the flux
s	Solver object of type HyPar
t	Current simulation time

Definition at line 16 of file FPPowerSystem3BusAdvection.c.

{
  HyPar *solver = (HyPar*) s;
  _ArrayCopy1D_(u,f,solver->npoints_local_wghosts);
  return(0);
}

int FPPowerSystem3BusDiffusion	(	double *	f,
		double *	u,
		int	dir1,
		int	dir2,
		void *	s,
		double	t
	)

Compute the dissipation term for the FPPowerSystem3Bus system

Parameters

f	Array to hold the computed dissipation term vector (same layout as u)
u	Array with the solution vector
dir1	First spatial dimension for the dissipation term being computed
dir2	Second spatial dimension for the dissipation term being computed
s	Solver object of type HyPar
t	Current simulation time

Definition at line 15 of file FPPowerSystem3BusDiffusion.c.

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

  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  static int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
  static double dissipation[_MODEL_NDIMS_*_MODEL_NDIMS_]; 

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

  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);

  int done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
    FPPowerSystem3BusDissipationFunction(dir1,dir2,params,t,dissipation);
    for (v = 0; v < _MODEL_NVARS_; v++) f[_MODEL_NVARS_*p+v] = dissipation[dir1*_MODEL_NDIMS_+dir2] * u[_MODEL_NVARS_*p+v];
    _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
  }

  return(0);
}

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

Compute the upwind flux at the interface, based on the drift velocity at that interface, from the left and right biased approximations to the interface flux. The drift (advection) velocity is multiplied to the solution in this function to get the advective flux.

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 (x or y)
s	Solver object of type HyPar
t	Current solution time

Definition at line 18 of file FPPowerSystem3BusUpwind.c.

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

  return(0);
}

int FPPowerSystem3BusInitialize	(	void *	s,
		void *	m
	)

Initialize the 3-bus power system model: Sets the default parameters, read in and set physics-related parameters, and set the physics-related function pointers in HyPar.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 27 of file FPPowerSystem3BusInitialize.c.

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

  static int count = 0;

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

  double pi = 4.0*atan(1.0);
  /* default values of model parameters */
  physics->N            = 2;
  physics->Pm1_avg      = 0.8;
  physics->Pm2_avg      = 1.6;
  physics->Pmref_avg    = 0.79330781761651;
  physics->H1           = 3.20;
  physics->H2           = 6.40;
  physics->Href         = 13.60;
  physics->E1           = 1.01556070860155;
  physics->E2           = 1.0491099265981;
  physics->Eref         = 1.05623172878954;
  physics->omegaB       = 2*pi*60.0;
  physics->sigma[0][0]  = 0.0125;
  physics->sigma[0][1]  = 0.0;
  physics->sigma[1][0]  = 0.0;
  physics->sigma[1][1]  = 0.0125;
  physics->lambda[0][0] = 10.0/physics->omegaB;
  physics->lambda[0][1] = 0.0;
  physics->lambda[1][0] = 0.0;
  physics->lambda[1][1] = 10.0/physics->omegaB;
  physics->gamma        = 0.25;

  physics->G   = (double*) calloc (3*3,sizeof(double));
  physics->B   = (double*) calloc (3*3,sizeof(double));
  
  physics->G[0*3+0] = 0.276805493111691;
  physics->G[0*3+1] = 0.213024867595501;
  physics->G[0*3+2] = 0.209205876527443;
  physics->G[1*3+0] = 0.213024867595501;
  physics->G[1*3+1] = 0.419642083051144;
  physics->G[1*3+2] = 0.286592141665043;
  physics->G[2*3+0] = 0.209205876527443;
  physics->G[2*3+1] = 0.286592141665044;
  physics->G[2*3+2] = 0.844559256324453;
  
  physics->B[0*3+0] = -2.36794416971567;
  physics->B[0*3+1] =  1.08817493992579;
  physics->B[0*3+2] =  1.22601259339234;
  physics->B[1*3+0] =  1.08817493992579;
  physics->B[1*3+1] = -2.72352378723346;
  physics->B[1*3+2] =  1.51348094527252;
  physics->B[2*3+0] =  1.22601259339234;
  physics->B[2*3+1] =  1.51348094527252;
  physics->B[2*3+2] = -2.98729895217208;

  /* reading physical model specific inputs - all processes */
  FILE *in;
  if ((!mpi->rank) && (!count)) printf("Reading physical model inputs from file \"physics.inp\".\n");
  in = fopen("physics.inp","r");
  if (!in) {
    if (!mpi->rank) fprintf(stderr,"Error: File \"physics.inp\" not found.\n");
    return(1);
  } 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,"Pm1_avg"   ))  {ferr=fscanf(in,"%lf",&physics->Pm1_avg   ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"Pm2_avg"   ))  {ferr=fscanf(in,"%lf",&physics->Pm2_avg   ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"Pmref_avg" ))  {ferr=fscanf(in,"%lf",&physics->Pmref_avg ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"H1"        ))  {ferr=fscanf(in,"%lf",&physics->H1        ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"H2"        ))  {ferr=fscanf(in,"%lf",&physics->H2        ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"Href"      ))  {ferr=fscanf(in,"%lf",&physics->Href      ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"E1"        ))  {ferr=fscanf(in,"%lf",&physics->E1        ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"E2"        ))  {ferr=fscanf(in,"%lf",&physics->E2        ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"Eref"      ))  {ferr=fscanf(in,"%lf",&physics->Eref      ) ;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,"gamma"     ))  {ferr=fscanf(in,"%lf",&physics->gamma     ) ;if(ferr!=1)return(1);}
        else if (!strcmp(word,"sigma"  ))  {
          ferr=fscanf(in,"%lf",&physics->sigma[0][0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->sigma[0][1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->sigma[1][0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->sigma[1][1]) ;if(ferr!=1)return(1);
        } else if (!strcmp(word,"lambda"))  {
          ferr=fscanf(in,"%lf",&physics->lambda[0][0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->lambda[0][1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->lambda[1][0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->lambda[1][1]) ;if(ferr!=1)return(1);
        } else if (!strcmp(word,"G"))  {
          ferr=fscanf(in,"%lf",&physics->G[0*3+0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[0*3+1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[0*3+2]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[1*3+0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[1*3+1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[1*3+2]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[2*3+0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[2*3+1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->G[2*3+2]) ;if(ferr!=1)return(1);
        } else if (!strcmp(word,"B"))  {
          ferr=fscanf(in,"%lf",&physics->B[0*3+0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[0*3+1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[0*3+2]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[1*3+0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[1*3+1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[1*3+2]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[2*3+0]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[2*3+1]) ;if(ferr!=1)return(1);
          ferr=fscanf(in,"%lf",&physics->B[2*3+2]) ;if(ferr!=1)return(1);
        }
      }
      } else {
        if (!mpi->rank) fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
      return(1);
      }
  }
  fclose(in);

  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
    if (!mpi->rank) {
      fprintf(stderr,"Error in FPPowerSystem3BusInitialize: 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         = FPPowerSystem3BusComputeCFL;
  solver->ComputeDiffNumber  = FPPowerSystem3BusComputeDiffNumber;
  solver->FFunction          = FPPowerSystem3BusAdvection;
  solver->HFunction          = FPPowerSystem3BusDiffusion;
  solver->Upwind             = FPPowerSystem3BusUpwind;

  count++;
  return(0);
}