eschnett/cactusamrex - Change GSGI6HIWST43XFEORCMKH6VPEGXHEZG5EK4JUUNWP3XFYUKGNA3AC

HydroToyAMReX: Make fluxes work

Created by Erik Schnetter on October 10, 2019

GSGI6HIWST43XFEORCMKH6VPEGXHEZG5EK4JUUNWP3XFYUKGNA3AC

Dependencies

In channels

master

main

Change contents

File addition: hydro-shocktube.par (-rw-r--r--)

[2.614]

ActiveThorns = "
    AMReX
    Formaline
    IOUtil
    HydroToyAMReX
    SystemTopology
    TimerReport
"
 
$nlevels = 1
$ncells = 32
Cactus::cctk_show_schedule = no
# Cactus::terminate = "iteration"
# Cactus::cctk_itlast = 0
Cactus::terminate = "time"
Cactus::cctk_final_time = 2.0
AMReX::verbose = yes
AMReX::xmin = -1.0
AMReX::ymin = -1.0
AMReX::zmin = -1.0
AMReX::xmax = +1.0
AMReX::ymax = +1.0
AMReX::zmax = +1.0
AMReX::ncells_x = $ncells
AMReX::ncells_y = $ncells
AMReX::ncells_z = $ncells
AMReX::max_num_levels = $nlevels
AMReX::regrid_every = 16
AMReX::regrid_error_threshold = 0.1
AMReX::dtfac = 0.5
HydroToyAMReX::setup = "shock tube"
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 32
# AMReX::out_tsv = yes
TimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32
TimerReport::out_filename = "TimerReport"
TimerReport::output_all_timers_together = yes
TimerReport::output_all_timers_readable = yes
TimerReport::n_top_timers = 50

File addition: hydro-sphericalshock.par (-rw-r--r--)

[2.614]

ActiveThorns = "
    AMReX
    Formaline
    IOUtil
    HydroToyAMReX
    SystemTopology
    TimerReport
"
 
$nlevels = 1
$ncells = 32
Cactus::cctk_show_schedule = no
# Cactus::terminate = "iteration"
# Cactus::cctk_itlast = 0
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
AMReX::verbose = no
AMReX::xmin = -1.0
AMReX::ymin = -1.0
AMReX::zmin = -1.0
AMReX::xmax = +1.0
AMReX::ymax = +1.0
AMReX::zmax = +1.0
AMReX::ncells_x = $ncells
AMReX::ncells_y = $ncells
AMReX::ncells_z = $ncells
AMReX::max_num_levels = $nlevels
AMReX::regrid_every = 16
AMReX::regrid_error_threshold = 0.1
AMReX::dtfac = 0.25
HydroToyAMReX::setup = "spherical shock"
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 32
# AMReX::out_tsv = yes
TimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32
TimerReport::out_filename = "TimerReport"
TimerReport::output_all_timers_together = yes
TimerReport::output_all_timers_readable = yes
TimerReport::n_top_timers = 50

Insertion in HydroToyAMReX/param.ccl at line 7 [3.5941]
[2.2666]
[2.2666]
```
  "shock tube" :: ""
  "spherical shock" :: ""
```

Insertion in HydroToyAMReX/param.ccl at line 16 [3.5941]

[2.2762]

CCTK_REAL shock_radius "Shock radius" STEERABLE=never
{
  0.0:* :: ""
} 0.1

Insertion in HydroToyAMReX/param.ccl at line 22 [3.5941]
[2.2764]
[2.2764]
Replacement in HydroToyAMReX/schedule.ccl at line 12 [3.6006]
B:BD[3.6112] → [2.2933:2965]
```
  WRITES: conserved(everywhere)
```
[3.6112]
[3.6112]
```
  WRITES: conserved(interior)
  SYNC: conserved
```

Insertion in HydroToyAMReX/schedule.ccl at line 40 [3.6006]

[2.3553]

SCHEDULE HydroToyAMReX_Output AT poststep
{
  LANG: C
  READS: conserved(everywhere)
} "Output grid data"

Insertion in HydroToyAMReX/schedule.ccl at line 61 [3.6006]

[3.6230]


# SCHEDULE HydroToyAMReX_Evolve1 AT evol
# {
#   LANG: C
#   READS: conserved_p(everywhere)
#   READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)
#   WRITES: conserved(everywhere)
# } "Evolve the hydro equations, step 1"
# 
# SCHEDULE HydroToyAMReX_Pressure AT evol AFTER HydroToyAMReX_Evolve1
# {
#   LANG: C
#   READS: conserved(everywhere)
#   WRITES: primitive(everywhere)
# } "Calculate pressure"
# 
# SCHEDULE HydroToyAMReX_Fluxes AT evol AFTER HydroToyAMReX_Pressure
# {
#   LANG: C
#   READS: conserved(everywhere)
#   READS: primitive(everywhere)
#   WRITES: flux_x(interior) flux_y(interior) flux_z(interior)
#   SYNC: flux_x flux_y flux_z
# } "Calculate the hydro fluxes"
# 
# SCHEDULE HydroToyAMReX_Evolve AT evol AFTER HydroToyAMReX_Fluxes
# {
#   LANG: C
#   READS: conserved_p(everywhere)
#   READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)
#   WRITES: conserved(everywhere)
# } "Evolve the hydro equations, step 2"

Replacement in HydroToyAMReX/src/hydrotoy.cxx at line 29 [3.6253]

B:BD[2.4140] → [2.4140:4208]

    Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {

[2.4140]

[2.4208]

    Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {

Replacement in HydroToyAMReX/src/hydrotoy.cxx at line 39 [3.6253]

B:BD[2.4391] → [2.4391:4459]

    Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {

[2.4391]

[2.4459]

    Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {

Replacement in HydroToyAMReX/src/hydrotoy.cxx at line 44 [3.6253]

B:BD[2.4588] → [2.4588:4613]

      etot[p.idx] = 1.0;

[2.4588]

[2.4613]

      etot[p.idx] = 1.0; // should add kinetic energy here
    });
  } else if (CCTK_EQUALS(setup, "shock tube")) {
    Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
      if (p.x <= 0.0) {
        rho[p.idx] = 2.0;
        momx[p.idx] = 0.0;
        momy[p.idx] = 0.0;
        momz[p.idx] = 0.0;
        etot[p.idx] = 2.0;
      } else {
        rho[p.idx] = 1.0;
        momx[p.idx] = 0.0;
        momy[p.idx] = 0.0;
        momz[p.idx] = 0.0;
        etot[p.idx] = 1.0;
      }
    });
  } else if (CCTK_EQUALS(setup, "spherical shock")) {
    Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
      CCTK_REAL r2 = pow(p.x, 2) + pow(p.y, 2) + pow(p.z, 2);
      if (r2 <= pow(shock_radius, 2)) {
        rho[p.idx] = 2.0;
        momx[p.idx] = 0.0;
        momy[p.idx] = 0.0;
        momz[p.idx] = 0.0;
        etot[p.idx] = 2.0;
      } else {
        rho[p.idx] = 1.0;
        momx[p.idx] = 0.0;
        momy[p.idx] = 0.0;
        momz[p.idx] = 0.0;
        etot[p.idx] = 1.0;
      }

Replacement in HydroToyAMReX/src/hydrotoy.cxx at line 195 [3.6253]

B:BD[2.8558] → [2.8558:8751]

      auto lambda_m = -dx / dt;
      auto lambda_p = dx / dt;
      auto var_m = u(p.I - p.DI(0));
      auto var_p = u(p.I);
      auto flux_m = f(p.I - p.DI(0));
      auto flux_p = f(p.I);

[2.8558]

[2.8751]

      auto I_m = p.I - p.DI(0);
      auto I_p = p.I;
      auto lambda_m = 1.0;
      auto lambda_p = -1.0;
      auto var_m = u(I_m);
      auto var_p = u(I_p);
      auto flux_m = f(I_m);
      auto flux_p = f(I_p);

Replacement in HydroToyAMReX/src/hydrotoy.cxx at line 219 [3.6253]

B:BD[2.9377] → [2.9377:9570]

      auto lambda_m = -dy / dt;
      auto lambda_p = dy / dt;
      auto var_m = u(p.I - p.DI(1));
      auto var_p = u(p.I);
      auto flux_m = f(p.I - p.DI(1));
      auto flux_p = f(p.I);

[2.9377]

[2.9570]

      auto I_m = p.I - p.DI(1);
      auto I_p = p.I;
      auto lambda_m = 1.0;
      auto lambda_p = -1.0;
      auto var_m = u(I_m);
      auto var_p = u(I_p);
      auto flux_m = f(I_m);
      auto flux_p = f(I_p);

Replacement in HydroToyAMReX/src/hydrotoy.cxx at line 243 [3.6253]

B:BD[2.10196] → [2.10196:10389]

      auto lambda_m = -dz / dt;
      auto lambda_p = dz / dt;
      auto var_m = u(p.I - p.DI(2));
      auto var_p = u(p.I);
      auto flux_m = f(p.I - p.DI(2));
      auto flux_p = f(p.I);

[2.10196]

[2.10389]

      auto I_m = p.I - p.DI(2);
      auto I_p = p.I;
      auto lambda_m = 1.0;
      auto lambda_p = -1.0;
      auto var_m = u(I_m);
      auto var_p = u(I_p);
      auto flux_m = f(I_m);
      auto flux_p = f(I_p);

Insertion in HydroToyAMReX/src/hydrotoy.cxx at line 265 [3.6253]

[3.7057]


extern "C" void HydroToyAMReX_Evolve1(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

Insertion in HydroToyAMReX/src/hydrotoy.cxx at line 270 [3.6253]

[3.7058]

  const CCTK_REAL dt = CCTK_DELTA_TIME;
  const CCTK_REAL dx = CCTK_DELTA_SPACE(0);
  const CCTK_REAL dy = CCTK_DELTA_SPACE(1);
  const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> rho_p_(cctkGH, rho_p);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momx_p_(cctkGH, momx_p);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momy_p_(cctkGH, momy_p);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momz_p_(cctkGH, momz_p);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> etot_p_(cctkGH, etot_p);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxrho_(cctkGH, fxrho);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomx_(cctkGH, fxmomx);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomy_(cctkGH, fxmomy);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomz_(cctkGH, fxmomz);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxetot_(cctkGH, fxetot);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fyrho_(cctkGH, fyrho);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fymomx_(cctkGH, fymomx);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fymomy_(cctkGH, fymomy);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fymomz_(cctkGH, fymomz);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 1> fyetot_(cctkGH, fyetot);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzrho_(cctkGH, fzrho);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzmomx_(cctkGH, fzmomx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzmomy_(cctkGH, fzmomy);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzmomz_(cctkGH, fzmomz);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 0> fzetot_(cctkGH, fzetot);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> rho_(cctkGH, rho);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> momx_(cctkGH, momx);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> momy_(cctkGH, momy);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> momz_(cctkGH, momz);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> etot_(cctkGH, etot);
  // Transport
  // dt rho + d_i (rho vel^i) = 0
  // dt mom_j + d_i (mom_j vel^i) = 0
  // dt etot + d_i (etot vel^i) = 0
  Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    auto mkupdate{[&](auto &fx, auto &fy, auto &fz) {
      return 0.5 * dt *
             ((fx(p.I + p.DI(0)) - fx(p.I)) / dx +
              (fy(p.I + p.DI(1)) - fy(p.I)) / dy +
              (fz(p.I + p.DI(2)) - fz(p.I)) / dz);
    }};
    rho_(p.I) = rho_p_(p.I) - mkupdate(fxrho_, fyrho_, fzrho_);
    momx_(p.I) = momx_p_(p.I) - mkupdate(fxmomx_, fymomx_, fzmomx_);
    momy_(p.I) = momy_p_(p.I) - mkupdate(fxmomy_, fymomy_, fzmomy_);
    momz_(p.I) = momz_p_(p.I) - mkupdate(fxmomz_, fymomz_, fzmomz_);
    etot_(p.I) = etot_p_(p.I) - mkupdate(fxetot_, fyetot_, fzetot_);
  });
}

Insertion in HydroToyAMReX/src/hydrotoy.cxx at line 395 [3.6253]

[2.13583]

}
extern "C" void HydroToyAMReX_Output(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> rho_(cctkGH, rho);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momx_(cctkGH, momx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momy_(cctkGH, momy);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> momz_(cctkGH, momz);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 1> etot_(cctkGH, etot);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxrho_(cctkGH, fxrho);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomx_(cctkGH, fxmomx);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomy_(cctkGH, fxmomy);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxmomz_(cctkGH, fxmomz);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> fxetot_(cctkGH, fxetot);
  Loop::loop_all<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    if (p.j == 0 && p.k == 0) {
#pragma omp critical(HydroToyAMReX_Output)
      {
        cout << "[" << p.i << "," << p.j << "," << p.k << "] (" << p.x << ","
             << p.y << "," << p.z << ") rho=" << rho_(p.I)
             << " momx=" << momx_(p.I) << " etot=" << etot_(p.I) << "\n";
      }
    }
  });
  Loop::loop_all<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    if (p.j == 0 && p.k == 0) {
#pragma omp critical(HydroToyAMReX_Output)
      {
        cout << "[" << p.i << "," << p.j << "," << p.k << "] (" << p.x << ","
             << p.y << "," << p.z << ") fxrho=" << fxrho_(p.I)
             << " fxmomx=" << fxmomx_(p.I) << " fxetot=" << fxetot_(p.I)
             << "\n";
      }
    }
  });