C++ input::PRef方法代码示例

/*

In the diagram below Ui represents the interior state. This should be represnted by (*this). Ug1 and Ug2 are the first and second layers of ghost cells respectively. They
are the output of the function. The input "layer" determines whether Ug1 or Ug2 is output for BCs that use extrapolation. Reflected boundaries (slipWall, viscousWall) do 
not use the layer input. Instead Ui+1 must be specified as (*this) to output Ug2.

____________________|___________
|         |         |          |
|         |         |          |
|   Ug2   |   Ug1   |   Ui     |
|         |         |          |
|         |         |          |
|_________|_________|__________|
                    |
             boundary face

Currently the following boundary conditions are supported: slipWall, viscousWall, characteristic, stagnationInlet, pressureOutlet, subsonicInflow, subsonicOutflow,
supersonicInflow, supersonicOutflow

*/
primVars primVars::GetGhostState( const string &bcType, const vector3d<double> &areaVec, const string &surf, const input &inputVars, const idealGas &eqnState, const int layer )const{
  //bcType -- type of boundary condition to supply ghost cell for
  //areaVec -- area vector of boundary face
  //surf -- i, j, k surface of boundary
  //inputVar -- all input variables
  //eqnState -- equation of state
  //layer -- layer of ghost cell to return (first (closest) or second (farthest))
  //the instance of primVars being acted upon should be the interior cell bordering the boundary

  primVars ghostState = (*this);  //set ghost state equal to boundary state to start

  //check to see that ghost layer corresponds to allowable number
  if ( !(layer == 1 || layer == 2) ){
    cerr << "ERROR: Error in primVars::GetGhostState. Requesting ghost state at a ghost layer " << layer << ". Please choose either 1 or 2" << endl;
    exit(0);
  }

  //normalize area vector (should always point out of domain)
  vector3d<double> normArea;
  if (surf == "il" || surf == "jl" || surf == "kl"){
    normArea = -1.0 * areaVec / areaVec.Mag(); //at lower surface normal should point out of domain for ghost cell calculation
  }
  else if (surf == "iu" || surf == "ju" || surf == "ku"){
    normArea = areaVec / areaVec.Mag(); 
  }

  double normVelCellCenter = 0;

  //slip wall boundary condition ----------------------------------------------------------------------------------------------------------------
  //slipWall is implemented as a reflection of the interior states so that there is no mass flow across the boundary.
  if (bcType == "slipWall"){             //for slip wall state should be reflected across boundary face, density and pressure stay equal to the boundary cell
    vector3d<double> stateVel = (*this).Velocity();
    normVelCellCenter = stateVel.DotProd(normArea);

    //for a slip wall the velocity of the boundary cell center is reflected across the boundary face to get the velocity at the ghost cell center
    vector3d<double> ghostVel (stateVel.X() - 2.0 * normArea.X() * normVelCellCenter,
			       stateVel.Y() - 2.0 * normArea.Y() * normVelCellCenter,
			       stateVel.Z() - 2.0 * normArea.Z() * normVelCellCenter);

    ghostState.data[1] = ghostVel.X();
    ghostState.data[2] = ghostVel.Y();
    ghostState.data[3] = ghostVel.Z();

    //numerical BCs for rho and pressure, same as boundary state

  }
  //viscous wall boundary condition ----------------------------------------------------------------------------------------------------------------
  //viscous wall uses the interior density and pressure, but flips the sign on the velocity so that the velocity at the boundary is 0.
  else if (bcType == "viscousWall"){             //for viscous wall velocity at face should be 0.0, density and pressure stay equal to the boundary cell
    vector3d<double> stateVel = (*this).Velocity();

    //ghost cell velocity at cell center is set to opposite of velocity at boundary cell center so that velocity at face will be zero
    vector3d<double> ghostVel (-1.0 * stateVel.X(),
			       -1.0 * stateVel.Y(),
			       -1.0 * stateVel.Z() );

    ghostState.data[1] = ghostVel.X();
    ghostState.data[2] = ghostVel.Y();
    ghostState.data[3] = ghostVel.Z();

    //numerical BCs for rho and pressure, same as boundary state

  }
  //subsonic inflow boundary condition --------------------------------------------------------------------------------------------------------------
  //this boundary condition enforces density and velocity as freestream inputs (constant) and extrapolates from the interior state to get pressure
  //this is a primative implementation, stagnationInlet or characteristic are better options
  else if (bcType == "subsonicInflow"){     //set velocity and density to freestream values
    double sos = eqnState.GetSoS( inputVars.PRef(), inputVars.RRef() );
    vector3d<double> ghostVel = inputVars.VelRef() / sos;   //nondimensionalize velocity

    ghostState.data[0] = 1.0;
    ghostState.data[1] = ghostVel.X();
    ghostState.data[2] = ghostVel.Y();
    ghostState.data[3] = ghostVel.Z();

    //numerical bc for pressure, same as boundary state

    if (layer == 2){ //extrapolate to get ghost state at 2nd layer
      ghostState = 2.0 * ghostState - (*this);
    }
//.........这里部分代码省略.........