FPPowerSystem3BusUpwind.c File Reference

Upwinding function for the 3-bus power system model. More...

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <physicalmodels/fppowersystem3bus.h>
#include <hypar.h>

Go to the source code of this file.

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

Detailed Description

Upwinding function for the 3-bus power system model.

Author

Debojyoti Ghosh

Definition in file FPPowerSystem3BusUpwind.c.

Function Documentation

int FPPowerSystem3BusDriftFunction	(	int	dir,
		void *	p,
		double *	x,
		double	t,
		double *	drift
	)

Compute the drift (advection) coefficients for the 3-bus power system

Parameters

dir	Spatial dimension (not used)
p	Object of type FPPowerSystem3Bus
x	Spatial coordinates
t	Current simulation time
drift	Array to hold the drift velocities

Definition at line 76 of file FPPowerSystem3BusFunctions.c.

{
  FPPowerSystem3Bus *params = (FPPowerSystem3Bus*) p;

  double theta1 = x[0];
  double theta2 = x[1];
  double Omega1 = x[2];
  double Omega2 = x[3];
  
  double omegaB     = params->omegaB;
  double Pm1_avg    = params->Pm1_avg;
  double Pm2_avg    = params->Pm2_avg;
  double Pmref_avg  = params->Pmref_avg;
  double H1         = params->H1;
  double H2         = params->H2;
  double Href       = params->Href;
  double gamma      = params->gamma;
  
  double Pe1, Pe2, Peref;
  ComputeElectricalPower(theta1,theta2,params,&Pe1,&Pe2,&Peref);

  double F1 = Pm1_avg / (2*H1) - Pmref_avg / (2*Href);
  double F2 = Pm2_avg / (2*H2) - Pmref_avg / (2*Href);
  double S1 = Pe1 / (2*H1) - Peref / (2*Href);
  double S2 = Pe2 / (2*H2) - Peref / (2*Href);

  drift[0] = omegaB * Omega1;
  drift[1] = omegaB * Omega2;
  drift[2] = F1 - gamma*Omega1 - S1;
  drift[3] = F2 - gamma*Omega2 - S2;

  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);
}