Pressure shifted from the input file

[AnyGa]

Hello! I am running a 3D hydrodynamical simulation. This is based on the jet.cpp file, but I converted it to a spherical coordinate system and added a point mass at the centre to see how the stream evolves as it approaches the BH. However, the stream pressure was kept several orders of magnitude higher than I set in the input file once they were injected into the simulation box.

Here is part of my cpp file:

void MeshBlock::ProblemGenerator(ParameterInput *pin) {
  gm1 = peos->GetGamma() - 1.0;
  #pragma omp simd
    for (int k=ks; k<=ke; ++k) {
      for (int j=js; j<=je; ++j) {
       for (int i=is; i<=ie; ++i) {
        phydro->u(IDN,k,j,i) = d_amb;
        phydro->u(IM1,k,j,i) = d_amb*vx_amb;
        phydro->u(IM2,k,j,i) = d_amb*vy_amb;
        phydro->u(IM3,k,j,i) = d_amb*vz_amb;
        if (NON_BAROTROPIC_EOS) {
          phydro->u(IEN,k,j,i) = p_amb/gm1
                                 + 0.5*d_amb*(SQR(vx_amb)+SQR(vy_amb)+SQR(vz_amb));
        }
      }
    }
  } 
}
void JetInnerX3(MeshBlock *pmb, Coordinates *pco, AthenaArray<Real> &prim, FaceField &b,
                Real time, Real dt,
                int il, int iu, int jl, int ju, int kl, int ku, int ngh) {
  Real t_start = 0.0;  
  Real t_end = 1e10;   
  bool jet_active = (time >= t_start && time <= t_end);
    
  // set primitive variables in inlet ghost zones
  #pragma omp simd
    for (int i=il; i<=iu; ++i) {
      for (int j=jl; j<=ju; ++j) {
          for (int k=1; k<=ngh; ++k) {
        Real zk = (pco->x1v(i))*std::cos(pco->x2v(j));
        Real ri = x1_0-(pco->x1v(i));	
        Real rad = std::sqrt(SQR(ri) + SQR(zk));
        
        if ( jet_active && r_jet < rad && rad <= 3*r_jet) {
            Real f = (rad-r_jet) / (2*r_jet);
            Real log_den = (1.0-f) * std::log(d_jet) + f * std::log(d_amb);
            prim(IDN,kl-k,j,i) = std::exp(log_den);
            Real log_pres = (1.0-f) * std::log(p_jet) + f * std::log(p_amb);
            prim(IPR,kl-k,j,i) = std::exp(log_pres);
            prim(IVX,kl-k,j,i) = vx_jet;
            prim(IVY,kl-k,j,i) = vy_jet;
            prim(IVZ,kl-k,j,i) = vz_jet;
        }
        if (jet_active && rad < r_jet) {
              prim(IPR,kl-k,j,i) = p_jet;
              prim(IDN,kl-k,j,i) = d_jet;
              prim(IVX,kl-k,j,i) = vx_jet;
              prim(IVY,kl-k,j,i) = vy_jet;
              prim(IVZ,kl-k,j,i) = vz_jet;
        } else {
            prim(IDN,kl-k,j,i) = d_amb;//d_amb;
            prim(IVX,kl-k,j,i) = prim(IVX,kl,j,i);//vx_amb;//
            prim(IVY,kl-k,j,i) = prim(IVY,kl,j,i);//vy_amb;//
            prim(IVZ,kl-k,j,i) = prim(IVZ,kl,j,i);//vz_amb;//
            prim(IPR,kl-k,j,i) = p_amb;//p_amb;//
        }
      }
    }
  }    

And some of the parameters in my input file are as follows:

<hydro>
gamma   = 1.6666666666667 # gamma = C_p/C_v
dfloor = 1e-6
pfloor = 1e-12

<problem>
four_pi_G  = 8.109059746158522e-09
GM = 4.999999999
r0 = 21.90268868177918
d = 2e-3
p = 1e-8
vx = 0.0
vy = 0.0
vz = 0.0

djet = 2.528072758537737
pjet = 2.006938980190448e-08
vxjet = -0.4761556
vyjet = 0
vzjet = 0.477792

rjet = 0.518834

And here is one of the outputs:
The pressure of the stream is about 1e-6 when they are injected, but in the input file, the 'pjet' is 2e-8.

I think the high ratio between the kinetic energy and internal energy may contribute to this, but no matter how I change the pressure of the inject stream, this phenomenon keeps appearing.

[tomidakn]

From your input file, I cannot tell your whole grid configuration. Does your plot include the ghost cells or only active cells? If not, use ghost_zones=true in your <output?> block. If the ghost cell values are correctly set, your boundary condition is working but the solver cannot handle it well - this may well be possible because of the extremely cold gas, as you mentioned by yourself.

Eulerian hydrodynamic codes are not good at handling very high Mach number. I guess you do not really need to set the pressure so low, and it is OK as long as the pressure is low enough and does not affect the jet dynamics.

[AnyGa]

Thank you for your reply.

I just enabled the ghost_zone and verified that the values in there match those in my input file.

I'm currently studying the evolution of the stream geometry and the possible shocks caused by vertical compression as the stream approaches the black hole, which is related to the balance between tidal force and gas pressure. After being shifted, the pressure can easily become comparable to the tidal force (~1e-5) in the simulation, especially in the region near the pericenter (pressure ~1e-4), which make the balance between them appears much earlier than the situation where the pressure is not shifted (evolves from 2e-8), resulting in significant difference (e.g. finding shock or not).

I was still wondering if there might be a way to reduce such shifts. Any suggestions would be greatly appreciated. Thank you again for the help!

[tomidakn]

Sorry for my late response, I am currently busy (well, as usual).

This is not an easy calculation. From your figure, I guess the injected flow near the boundary is in a quasi steady state. The injected low density flow immediately hits the upstream gas and compressed. It might help if you put the jet in the active zone in the initial condition. Also, you could try lower ambient density, although it may make the simulation unstable.

Is there any literature using similar setup (other than SPH)?

[AnyGa]

Thanks for the reply. It is true that the interaction between the stream elements and the resulting compress along their orbit direction could possibly result in the pressure change. However, before the pericenter, the stream is accelerated by the gravity, and injected elements along the orbit should be farther and farther before pericenter, thus would not be able to hit together (If I understand correctly).

I also tried different inputs through changing the units, and put the stream in the active zone using an user defined SourceFunction like this (here x1_0, x2_0, x3_0 is the inject position in a sphere_polar coordinate):

void FixJetRegionSource(
    MeshBlock *pmb, const Real time, const Real dt,
    const AthenaArray<Real> &prim, const AthenaArray<Real> &prim_scalar,
    const AthenaArray<Real> &bcc, AthenaArray<Real> &cons,
    AthenaArray<Real> &cons_scalar) {
    Real x0 = x1_0*std::sin(x2_0)*std::cos(x3_0);
    Real y0 = x1_0*std::sin(x2_0)*std::sin(x3_0);
    Real z0 = x1_0*std::cos(x2_0);
    for (int k = pmb->ks; k <= pmb->ke; ++k) {
        for (int j = pmb->js; j <= pmb->je; ++j) {
          for (int i = pmb->is; i <= pmb->ie; ++i) {
           Real x = pmb->pcoord->x1v(i) * std::sin(pmb->pcoord->x2v(j)) * std::cos(pmb->pcoord->x3v(k));
           Real y = pmb->pcoord->x1v(i) * std::sin(pmb->pcoord->x2v(j)) * std::sin(pmb->pcoord->x3v(k));
           Real z = pmb->pcoord->x1v(i) * std::cos(pmb->pcoord->x2v(j));
           Real rad = std::sqrt(SQR(x - x0) + SQR(y - y0) + SQR(z - z0));
            if (rad < r_jet) {
              cons(IPR, k, j, i) = p_jet;
              cons(IDN, k, j, i) = d_jet;
              cons(IM1, k, j, i) = d_jet * vx_jet;                   
              cons(IM2, k, j, i) = d_jet * vy_jet; 
              cons(IM3, k, j, i) = d_jet * vz_jet; 
              cons(IEN, k, j, i) = p_jet/ gm1+0.5 * d_jet * ( vx_jet*vx_jet + vy_jet*vy_jet + vz_jet*vz_jet );
            }
           } 
          }
        }
      }
    return;
    }

It turns out using higher jet density (with specific internal energy fixed) input could help reduce the pressure change, but the problem still exist and, of course, the simulation gets more unstable, as dt becomes smaller or changes more drastically.

Huang et al. 2024 did a similar simulation as mine, injeting gas in several Rperi along a eccentric orbit around BH but considering the radiation transfer also. They don't seem to have the same problems, which may possibly attribute to a higher temperature choice (their work ~1e5 K; mine ~1e3 - 1e4 K), or the radiation transfer module. I tried radiation module, but the simulation becomes very slow (obviously), and the pressure seems to getting much higher.

[tomidakn]

Sorry for the confusion, but what I meant by "upstream" is the ambient gas. The jet injected from the boundary hits the envelope gas and forms a shock at its tip. For this shock, the ambient is the upstream.

The ram pressure from the ambient is O(1e-3) in your configuration. This is larger than the jet pressure, and therefore the jet gets easily compressed by the interaction with the envelope. I guess you need to lower the ambient density and/or increase the jet pressure. I do not think the radiation transfer is directly relevant.

@swdavis, can you provide some advice?

[AnyGa]

Thank you for the further explanation, now I understand better. There may be some misleading information in my configuration documents. In fact, I had set the ambient pressure to 1e-8, and the injected stream's pressure to 2e-8. Of course, the difference between these two values is indeed very small.

Therefore, I tried another run with Pamb ~ 1e-9, damb ~ 1e-5; Pjet ~ 1e-7, djet ~ 1, and set the ambient velocity to zero. Attached is the corresponding snapshot taken when the stream has just been injected into the active zone.

The ring-like structure appearing in the center of the diagram is also puzzling. This might possibly be related to the GM flag being turned on in the section of the input file. However, aside from the inner boundary in the X3 direction (phi) being user-defined, where gas outside rjet is set to ambient conditions with zero velocity, all other boundary conditions are set to 'outflow'. This makes me wonder why the gas still accumulates in such a ring-like pattern, leading to an increase in Pamb.

I also attempted to lower pamb (and accordingly pfloor) to around 1e-10 or even 1e-12, but the overall behavior did not change significantly, while the simulation became increasingly unstable.

[tomidakn]

The structure near the center is likely due to the gravity and boundary condition. Athena++'s outflow boundary conditions are not truly outflowing. They are just "zero gradient". So some structures may form near the boundary. Please check the velocity field. There might be a shock propagating outward from the inner boundary.

[AnyGa]

Thank you! I just checked the velocity field and found that this ring-like structure is formed by compressed ambient gas. The strong tidal forces near the pericenter compress the ambient gas in the vertical direction, triggering a shock and causing an increase in pressure. As the pressure rises, the gas near the pericenter tries to expand in all directions, leading to collisions with the ambient gas falling in from larger radii. This interaction causes the structure to extend also in the radial direction.

Interestingly, this behaviour, the shock induced by compression near the pericenter, is exactly what I had hoped to observe in the injected stream. Unfortunately, as shown above, all attempts so far have failed to reproduce such dynamics.