Tutorial

This page will explain the set-up, run and visualization of a typical case. For a quick overview of running the code, see Quickrun. At a glance, domain dimension and control parameters are handled in file inputs (or similar name) while problem description is in * prob.h

It is recommended to cp the tutorial into a wrk or exm/tmp directory. Is not neccesary, but if the code is cloned, changes in mian files will register as to commit. The previous directories will always be ignored by git.

Problem Set-up

The problem involves a heavy fluid descending into a lighter fluid, creating conditions for a Rayleigh-Taylor instability. The upper half of the domain contains a fluid with a density of 2, while the lower half is filled with a fluid of density 1. The initial pressure distribution is hydrostatic, and a velocity perturbation is introduced to trigger the instability The initial conditions can be summarised as

where Inline equation: $y_0$ is the middle of the domain, $P_0=2$ is the pressure and gravity is taken as 1 pointing downwards and the perturbation is $\epsilon=0.025$ . The initial conditions followed the paper by Shi et al . Both fluids are miscible and follow the ideal gas law with adiabatic coefficient of 5/3.

prob.h

Set-up of Problem Parameters

The parameters of the problem (not the domain) are wrapped into a ProbParm structure

// problem parameters
struct ProbParm {
  Real p_int = 2.0;
  Real rho_1 = 1.0;
  Real rho_2 = 2.0;
  Real grav = -1.0; 
  Real eps =  0.025;
};

Set-up of initial conditions

This is done in the function prob_initdata

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
prob_initdata(int i, int j, int k, Array4<Real> const &state,
              GeometryData const &geomdata, ProbClosures const &cls,
              ProbParm const &prob_parm) {
  const Real *prob_lo = geomdata.ProbLo();
  const Real *prob_hi = geomdata.ProbHi();
  const Real *dx = geomdata.CellSize();

  Real x = prob_lo[0] + (i + Real(0.5)) * dx[0];
  Real y = prob_lo[1] + (j + Real(0.5)) * dx[1];

where the spatial coordinates x and y of the cells are determined from the mesh sizes and problem dimensions (which are defined in inputs).

The initial conditions are defined as

  Real Pt, rhot, uxt,uyt;
  Real Lint = prob_hi[0] / 2;
  Real Pint = prob_parm.p_int; // interface Pressure

  const Real freq = Real(8)*Real(3.14159265359); // wavelength = x-domain

  Real yrel = y - Lint;
  Real delta= 0.2*Lint;   // region size where perturbation is significant
  Real delta2  = dx[1]/5; // transition region between top/bottom
  Real step = Real(0.5) + Real(0.5)*tanh(yrel/delta2);
  rhot = step*prob_parm.rho_2 + (Real(1.0) -step)*prob_parm.rho_1;
  Pt = Pint + rhot*prob_parm.grav*(y - Lint); // hydrostatic pressure

  uxt = Real(0.0);
  uyt = -prob_parm.eps*cos(freq*x)*aux; // perturbation in y-component

  state(i, j, k, cls.URHO) = rhot;
  state(i, j, k, cls.UMX) = rhot * uxt;
  state(i, j, k, cls.UMY) = rhot * uyt;
  Real et = Pt / (cls.gamma - Real(1.0));
  state(i, j, k, cls.UET) = et + Real(0.5) * rhot * (uxt * uxt + uyt * uyt);

The array state holds all the information required in conserved variables

to initialise the problem.

Solvers

The solvers is selected in the line

typedef rhs_dt<riemann_t<false, ProbClosures>, no_diffusive_t, user_source_t<ProbClosures>>
    ProbRHS;

which defines which problem to solve

Source Term

Function

template <typename cls_t>
class user_source_t {
  public:
  void inline src(const amrex::MFIter &mfi,
                  const amrex::Array4<const amrex::Real> &prims,
                  const amrex::Array4<amrex::Real> &rhs, const cls_t *cls_d,
                  amrex::Real dt){

The lines below loop over the cells and gravity is added to the momentum and energy equation. The prob_parm.grav reads the gravity value defined in ProbParm structure.

amrex::ParallelFor(bx,
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
    const auto& cls = *cls_d;
    Real rho = prims(i, j, k, cls.QRHO);
    rhs(i,j,k,cls.UMY) += prob_parm.grav*rho;
    rhs(i,j,k,cls.UET) += prob_parm.grav*rho*prims(i, j, k, cls.QV);
    });

AMR

This function defines which points need to be "tagged" to be refined. AMReX will then create patches around that cell to create the refined regions.

AMREX_GPU_DEVICE AMREX_FORCE_INLINE void
user_tagging(int i, int j, int k, int nt_level, auto &tagfab,
             const auto &sdatafab, const auto &geomdata,
             const ProbParm &prob_parm, int level) {

The actual tagging is produced in the follwoing lines, where the array tagfab(i,j,k) is set ot true or false based on the value of density.

  Real rhot = sdatafab(i,j,k,ProbClosures::URHO);
  if (nt_level > 0)
  {
  // tag cells based on density values
  switch (level)
   {
    case 0:
      tagfab(i,j,k) = (rhot > 1.1 && rhot < 1.9);
      break;
    case 1:
      tagfab(i,j,k) = (rhot > 1.2 && rhot < 1.8);
      break;
    default:
      tagfab(i,j,k) = (rhot > 1.3 && rhot < 1.7);
      break;
    }
  }

Input file

Time steps

max_step  = 20000
stop_time = 2.0
cfl = 0.3

Boundary conditions

The boudantry conditions are defined More complex boudnarty conditions can be defined

cns.lo_bc = -1 5
cns.hi_bc = -1 5

Mesh and refinement

amr.n_cell          = 128 512
amr.max_level       = 1     # maximum level number allowed
amr.regrid_int      = 2 2 # how often to regrid per level

These lines will select size of the mesh and how often to regrid

Compile and Running

To compile use

$ make

and run case using 2 cores (or any number you want)

$ mpirun -np 2 ./main2d.gnu.MPI.ex inputs

Post-processing

The input file is set up to plot every 100 steps

$ ls plot
plt00000	plt01300	plt02600	plt03900	plt05200
plt00100	plt01400	plt02700	plt04000	plt05300
plt00200	plt01500	plt02800	plt04100	plt05400
plt00300	plt01600	plt02900	plt04200	plt05500
plt00400	plt01700	plt03000	plt04300	plt05600
plt00500	plt01800	plt03100	plt04400	plt05700
plt00600	plt01900	plt03200	plt04500	plt05800
plt00700	plt02000	plt03300	plt04600	plt05900
plt00800	plt02100	plt03400	plt04700	plt06000
plt00900	plt02200	plt03500	plt04800	plt06100
plt01000	plt02300	plt03600	plt04900	plt06200
plt01100	plt02400	plt03700	plt05000	plt06300
plt01200	plt02500	plt03800	plt05100	plt06390

In the following sectionsm we will show how to visualize the results using Python, Visit and Paraview.

Python

The following assumes that yt is installed. Check Tips for installation. A Python script for easy of use,by invoking

$ python plot.py

The last density snaphot at t=2 will be saved into a png file, which looks like (with default color scheme in yt).

Visit

To see the results with VisIt, you can use the script cerisse, to open all directories at the same time (to make an animation for example). See Tips to set-up the script.

$ cerisse visit

This will create a file movie.visit in the tutorial directory. To open Visit, type visit (or use the appropiate icon).

Once opened, use File/Open to load the movie.visit file. If it loads correctly, it should look something like:

Adding a plot is easy by just Add/Pseudocolor/Density

The results shows the pseudocolor for density with default visualization options (see Visit manual and webpage for details on costumization). By clicking the arrow in the panel

we obtain a quick animation of the results. The final snapshot should be

By Add/Subset/Levels we can get an idea of the refinement

Double-cliking on the Subset/Levels we can edit the attributes of the plot

And the final plot should look like, where the solid lines are the limits of refinmement

There is no need to open the complete database and tiem steps can be opened individually by open the individual header files (for example opening directly plot/plt02600/Header) VisIt is good for exploring the data interactively, Visit can also be scripted with Python. For a VisIt tutorial check (http://visitusers.org/index.php?title=VisIt_Tutorial)

Paraview

To use Paraview, start Paraview in the usual way and open the plot\plt... folder

Select the AMReX/Boxlib Grid Reader

Tick the density box and click Apply, it would look something like

Select Density and Surface

And then use the final time to show the final plot (with default options)

Similar to Visit, Paraview is good for exploring the data interactively. can also be scripted with Python. Both softwares are similar and can use HPC, remote visualization and support large number of points (billions) A Paraview Tutorial manual can be found in Paraview Manual

PreviousQuickStart NextTips

Last updated 8 months ago