Controlling Initial Pressure

On this page of the wiki in the problem generators section it says that primitive variables (such as pressure) are automatically generated by the code. What I'm wondering is: how exactly is the code generating these primitive variables and is there any way they (specifically pressure) can be set as initial conditions in the problem generator.

So far I have tried using phydro->u(IEN,k,j,i) and phydro->w(IPR,k,j,i) to control the initial pressure but when I look at the t=0 output the pressure profile doesn't look like what I think it should.

I think the problem may be that I am trying to evolve an isothermal equilibrium using an adiabatic equation of state (ie there are adiabatic perturbations on an isothermal equilibrium) to see how instabilities might emerge. I think I may have to use a general equation of state to achieve this but I thought that since the isothermal part is only in the initial conditions I would be able to have a simple adiabatic equation of state.

I have attached a section of my problem generator and an image of the initial pressure (and density) profile of the t=0 output.

//========================================================================================
//! \fn void MeshBlock::ProblemGenerator(ParameterInput *pin)
//  \brief Problem Generator for the pressure bound cylinder test
//========================================================================================

void MeshBlock::ProblemGenerator(ParameterInput *pin) {
  Real gmma  = peos->GetGamma();
  Real gmma1 = gmma - 1.0;
  //Real delta_p = 0;
  // Read input parameters

  // Set paramters in ambient medium 
  Real Gamma_z = pin->GetReal("problem","Gamma_z");
  Real Gamma_phi =pin->GetReal("problem","Gamma_phi");

  Real sigma_sq1 = pin->GetReal("problem","sigma_sq1");
  Real sigma_sq2 = pin->GetReal("problem","sigma_sq2");
  Real P_inf = pin->GetReal("problem","P_inf");
  Real P_s = pin->GetReal("problem","P_s");
  Real r_0 = pin->GetReal("problem","r_0");

  Real size = pin->GetReal("problem","size");
  Real tol = pin->GetReal("problem","tol");
  Real rad = 0;
  Real press1 =0;

  Real rho;
  Real num;
  Real den;
  Real press;

  int step =0;

  //We want to find r_s which is the radius at which press = P_s.  Instead of finding the inverse of 
  //P(r) and solving for r_s we instead step through r and check P(r) to see if it's near P_s.
  //Because P(r) is monotonic decreasing (see our MATLAB program) we check if P((step-1)*size) or
  //P(step*size) is closer to P_s.  Then rad is an approximation of r_s within the tolerance set by
  //size (the step size).

  while (rad < 0.5) {
    //num = (-sigma_sq1 + sqrt(SQR(sigma_sq1) + (1/(2*M_PI))*P_inf*(SQR(Gamma_z)+SQR(Gamma_phi*step*size))));
    num = -sigma_sq1 + sqrt( SQR(sigma_sq1)  + 2*P_inf*(SQR(Gamma_z)+SQR(Gamma_phi*step*size)));
    //den = (1/(4*M_PI))*( SQR(Gamma_z) + SQR(Gamma_phi*step*size) );
    den = (SQR(Gamma_z)+ SQR(Gamma_phi*step*size));
    rho = num/den;
    press = sigma_sq1*rho;

    if (press - P_s < 0) {
      //num = (-sigma_sq1 + sqrt(SQR(sigma_sq1) + (1/(2*M_PI))*P_inf*(SQR(Gamma_z)+SQR(Gamma_phi*step*size))));
      num = -sigma_sq1 + sqrt( SQR(sigma_sq1)  + 2*P_inf*(SQR(Gamma_z)+SQR(Gamma_phi*(step-1)*size)));
      //den = (1/(4*M_PI))*( SQR(Gamma_z) + SQR(Gamma_phi*step*size) );
      den = (SQR(Gamma_z)+ SQR(Gamma_phi*(step-1)*size));
      rho = num/den;
      press1 = sigma_sq1*rho;

      if (abs(press - P_s) < abs(press1-P_s) ) {
        rad = 0.6;
      }
      else if (abs(press - P_s) >= abs(press1 - P_s) ) {
        rad = 0.7;
      }
    }

    step++;
  }

  if (rad = 0.6){
    rad = step*size;
  } else if (rad = 0.7) {
    rad = (step-1)*size; 
  }

  std::cout << size <<"\n";
  std::cout << press <<"\n";
  std::cout << rad << "\n";
  std::cout << step;

  Real B_z;
  Real B_phi;
  Real v_phi;

  Real B_x;
  Real B_y;

  Real v_x;
  Real v_y;

  // Initialize the grid
  for (int k=ks; k<=ke; k++) {
    for (int j=js; j<=je; j++) {
      for (int i=is; i<=ie; i++) {

        Real x = pcoord->x1v(i);
        Real y = pcoord->x2v(j);
        Real z = pcoord->x3v(k);
        Real r = sqrt(SQR(x)+SQR(y));

        // cloud interior
        //num = (-sigma_sq1 + sqrt(SQR(sigma_sq1) + (1/(2*M_PI))*P_inf*(SQR(Gamma_z)+SQR(Gamma_phi*r))));
        num = -sigma_sq1 + sqrt( SQR(sigma_sq1)  + 2*P_inf*(SQR(Gamma_z)+SQR(Gamma_phi*r)));
        //den = (1/(4*M_PI))*( SQR(Gamma_z) + SQR(Gamma_phi*r) );
        den = (SQR(Gamma_z)+ SQR(Gamma_phi*r));
        rho = num/den;
        press = sigma_sq1*rho;

        //v_phi = (Gamma_phi/std::sqrt(4*M_PI))*std::sqrt(SQR(r)*rho);
        v_phi = Gamma_phi*std::sqrt(SQR(r)*rho);

        v_x = -v_phi*std::sin(std::atan2(y,x));
        v_y = v_phi*std::cos(std::atan2(y,x));

        if (press<P_s) {
          press = P_s;
          rho = P_s/sigma_sq1;
          v_x = v_x*exp(-SQR(r-rad)/r_0);
          v_y = v_y*exp(-SQR(r-rad)/r_0);
        }

        phydro->u(IDN,k,j,i) = rho;
        phydro->u(IEN,k,j,i) = press/gmma1;
        phydro->w(IPR,k,j,i) = press;
        phydro->u(IM1,k,j,i) = v_x*rho;
        phydro->u(IM2,k,j,i) = v_y*rho;

      }
    }
  }

  return;
}

profile

My intention is that there should be a linear relationship between rho and press at t=0 but clearly there isn't. I'm not sure how to manipulate the initial conditions to allow that.

Any help would be greatly appreciated, thank you very much in advance.

Will

PrincetonUniversity / athena-public-version

Controlling Initial Pressure #75