PrincetonUniversity / athena-public-version

(MOVED) Athena++ GRMHD code and adaptive mesh refinement (AMR) framework. Active repository --->
https://github.com/PrincetonUniversity/athena
BSD 3-Clause "New" or "Revised" License
160 stars 118 forks source link

Disappearing Mass #36

Closed umgraftw closed 4 years ago

umgraftw commented 4 years ago

I have been attempting to model gravitational instability for filaments with an Ostriker density profile.

My simulations have recently encountered an interesting error where the mass goes to zero after the simulation runs for a period of time. I am using the FFT gravity method thus all the boundary conditions are periodic, so I don't think the mass is flowing out at the boundary.

I should be clear that athena does not output an error message, I am just under the impression that under periodic boundary conditions the mass should be conserved and even if it weren't there would be some time for the mass to `escape', it wouldn't happen instantaneously.

My current suspicion is that there is some maximum density which if exceeded will cause some variable to return NaN and in turn removes the mass somehow, or the self gravity returns some error due to an attempt to resolve sub-voxel structure and that removes the mass.

The .hst output I have returns NaN for the momenta, kinetic energy, and total energy after the mass `escapes'. It gives 0.00000e+00 for the gravitational energy.

Has anyone else encountered behaviour like this? I'm very new to simulation and am not entirely sure where to look to resolve this.

Thanks in advance.

c-white commented 4 years ago

Could you provide more details? It would help to have the problem generator, input file, and time by which things have gone wrong.

changgoo commented 4 years ago

Is this a gravitationally unstable case? It won't be surprising if density becomes NaN due to runaway collapse. Yes, we need more information.

umgraftw commented 4 years ago

Here is the problem generator I used. I'm trying to model adiabatic perturbations on top of a isothermal cylinder. The conceit is the cylinder is initially at some isothermal equilibrium (it has the Ostriker density profile) and there is an adiabatic perturbation on top. That's why we set the initial pressure to be the square of the isothermal sound speed times the density. The sinusoidal perturbation is basically an anzatz (we know how a perturbation would develop). The randomized density perturbations are to simulate anisotropy, they're from an older version where I was hoping the small variations would be enough to seed a collapse on their own. This problem generator is a modified version of a problem generator a colleague of mine was using to make a Bonnor Ebert sphere, which is a modified shock cloud.

//Bonnor Ebert problem generator

// C++ headers
#include <cmath>      // sqrt()
#include <iostream>   // endl
#include <sstream>    // stringstream
#include <stdexcept>  // runtime_error
#include <string>     // c_str()
#include <time.h>     // srand(time)

// Athena++ headers
#include "../athena.hpp"
#include "../athena_arrays.hpp"
#include "../bvals/bvals.hpp"
#include "../coordinates/coordinates.hpp"
#include "../eos/eos.hpp"
#include "../field/field.hpp"
#include "../hydro/hydro.hpp"
#include "../mesh/mesh.hpp"
#include "../parameter_input.hpp"
#include "../gravity/gravity.hpp"

namespace {
Real gmma1, dl, pl, ul;
Real bxl, byl, bzl;
Real gconst;
} // namespace

//========================================================================================
//! \fn void Mesh::InitUserMeshData(ParameterInput *pin)
//  \brief Function to initialize problem-specific data in mesh class.  Can also be used
//  to initialize variables which are global to (and therefore can be passed to) other
//  functions in this file.  Called in Mesh constructor.
//========================================================================================

void Mesh::InitUserMeshData(ParameterInput *pin) {

  if (SELF_GRAVITY_ENABLED) {
    Real fourPiG = pin->GetReal("problem","fourPiG");
    SetFourPiG(fourPiG);
    //Real eps = pin->GetOrAddReal("problem","grav_eps", 0.0);
    //SetGravityThreshold(eps);
    SetMeanDensity(0.0);
  }
  return;
}

//========================================================================================
//! \fn void MeshBlock::ProblemGenerator(ParameterInput *pin)
//  \brief Problem Generator for the shock-cloud interaction test
//========================================================================================

void MeshBlock::ProblemGenerator(ParameterInput *pin) {
  Real gmma  = peos->GetGamma();
  gmma1 = gmma - 1.0;

  // Read input parameters
  Real rad    = pin->GetReal("problem","radius");
  Real mach = pin->GetReal("problem","Mach");
  Real drat = pin->GetReal("problem","drat");
  Real iso_sp = pin->GetReal("hydro","iso_sound_speed");
  Real pil = pin->GetReal("problem","pil");
  Real maxn = pin->GetReal("problem","maxn");
  Real amp = pin->GetReal("problem","amp");

  // Set paramters in ambient medium 
  Real dr = std::pow((1+SQR(rad)/8),-2);
  Real pr = 1.0/(peos->GetGamma());
  Real ur = 0.0;

  // Initialize the grid
  for (int k=ks; k<=ke; k++) {
    for (int j=js; j<=je; j++) {
      for (int i=is; i<=ie; i++) {
         // ambient gas
         phydro->u(IDN,k,j,i) = dr;
         phydro->u(IM1,k,j,i) = ur*dr;
         phydro->u(IM2,k,j,i) = 0.0;
         phydro->u(IM3,k,j,i) = 0.0;
         //phydro->u(IEN,k,j,i) = pr/gmma1 + 0.5*dr*(ur*ur);
         phydro->u(IPR,k,j,i) = SQR(iso_sp)*phydro->u(IDN,k,j,i);

        // cloud interior
        Real diag = std::sqrt(SQR(pcoord->x1v(i)) + SQR(pcoord->x2v(j))
                              );
    Real profile =std::pow((1+SQR(diag)/8),-2);
        if (diag < rad) {
          phydro->u(IDN,k,j,i) = profile*(1+amp*std::cos(pil*maxn*pcoord->x3v(k)));
          phydro->u(IM1,k,j,i) = ur*dr*drat;
          phydro->u(IM2,k,j,i) = 0.0;
          phydro->u(IM3,k,j ,i) = 0.0;
          //phydro->u(IEN,k,j,i) = pr/gmma1 + 0.5*dr*drat*(ur*ur);
      phydro->u(IPR,k,j,i) = SQR(iso_sp)*phydro->u(IDN,k,j,i);
        }
    std::srand(time(0));    
    if (diag < rad/3) {
          phydro->u(IDN,k,j,i) = profile*(1+amp*std::cos(pil*maxn*pcoord->x3v(k)))+0.0001*(rand()/RAND_MAX);
          phydro->u(IM1,k,j,i) = ur*dr*drat;
          phydro->u(IM2,k,j,i) = 0.0;
          phydro->u(IM3,k,j ,i) = 0.0;
          //phydro->u(IEN,k,j,i) = pr/gmma1 + 0.5*dr*drat*(ur*ur);
      phydro->u(IPR,k,j,i) = SQR(iso_sp)*phydro->u(IDN,k,j,i);
        }
      }
    }
  }

  return;
}

//----------------------------------------------------------------------------------------

Here is the input file

<comment>
problem   = shock cloud interaction
reference = Shin,M.-S., Snyder, G., & Stone, J.M.
configure = --prob=shk_cloud_NoWave

<job>
problem_id = Cloud_NoWave      # problem ID: basename of output filenames

<output1>
file_type  = hst        # History data dump
dt         = 0.01       # time increment between outputs

<output2>
file_type  = hdf5       # data dump
variable   = prim       # variables to be output
dt         = 0.1        # time increment between outputs

<time>
cfl_number = 0.4        # The Courant, Friedrichs, & Lewy (CFL) Number
nlim       = 100000     # cycle limit
tlim       = 15.00      # time limit
integrator  = vl2       # time integration algorithm
xorder      = 2         # order of spatial reconstruction
ncycle_out  = 1         # interval for stdout summary info

<mesh>
nx1        = 200         # Number of zones in X1-direction
x1min      = -20.0       # minimum value of X1
x1max      = 20.0        # maximum value of X1
ix1_bc     = periodic     # inner-X1 boundary flag
ox1_bc     = periodic     # outer-X1 boundary flag

nx2        = 200         # Number of zones in X2-direction
x2min      = -20.0       # minimum value of X2
x2max      = 20.0        # maximum value of X2
ix2_bc     = periodic     # inner-X2 boundary flag
ox2_bc     = periodic     # outer-X2 boundary flag

nx3        = 400          # Number of zones in X3-direction
x3min      = -20.0       # minimum value of X3
x3max      = 20.0        # maximum value of X3
ix3_bc     = periodic    # inner-X3 boundary flag 
ox3_bc     = periodic    # outer-X3 boundary flag

<meshblock>
nx1 = 100
nx2 = 100
nx3 = 100

<hydro>
iso_sound_speed = 1.0 
gamma      = 1.3333    # gamma = C_p/C_v

<problem>
Mach       = 0.0       # Mach number of shock
drat       = 50000.0   # density ratio of cloud
beta       = 1.0       # plasma beta in upstream gas
fourPiG    = 1.0       #four pi g
radius     = 20.0      #cloud radius
maxn       = 4         #number of maxima
pil        = 0.07853981633 #kfactor
amp        = 0.125     #perturbed   

<meshblock>
nx1=100
nx2=100
nx3=100

The hope is that I'll be able to observe the collapse over a few freefall times, once the behaviour of the collapse in this simpler version starts making more sense to me I'll add more complexity (initial rotation and magnetic fields). Generally we use a polytropic equation of state where gamma is even closer to 1, but due to the way gamma is used to find the internal energy of the cloud I lowered it only a little from the 5/3 to 4/3. The mass tends to disappear around t=10, the instability has definitely had time to develop but the change is so drastic. I was hoping somehow there would be some indication that mass was escaping or the coalescing cores were getting too dense.

Thank you very much for your prompt replies. I'm a novice and every bit of advice is very helpful.

pdmullen commented 4 years ago

I found one issue in your problem generator. In Mesh::InitUserMeshData, you set the mean density (i.e., grav_mean_rho) to zero. Instead, it is imperative that you set this value to the mean density of the mesh.

The (default) FFT self-gravity solver in Athena++ (a) employs Jeans' swindle to solve for the gravitational potential subject to fully periodic boundary conditions and (b) applies a momentum conserving numerical technique to apply gravitational accelerations using "gravity fluxes". These fluxes are computed by taking the divergence of the gravitational stress tensor, thus, when in combination with Jeans' swindle, the (non-zero) mean density actually appears in the equations of motion.

After setting the mean density to the true mean density on the mesh, I was able to run your problem out to tlim = 15.00. The corresponding history file is attached.

Cloud_NoWave.hst.txt

umgraftw commented 4 years ago

Thank you very much, I have just finished a simulation out to tlim =15.00 as well and found that the mass did not disappear.

Thank you all for your suggestions.