Skip to content

Commit

Permalink
output local velocity gradients to a plotfile
Browse files Browse the repository at this point in the history
  • Loading branch information
@isriva
isriva committed Jul 26, 2023
1 parent 8196b93 commit 9e054f5
Show file tree
Hide file tree
Showing 3 changed files with 135 additions and 1 deletion.
5 changes: 5 additions & 0 deletions src_compressible_stag/compressible_functions_stag.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,11 @@ void WriteSpatialCross1D(const MultiFab& spatialCross , int step, const Geome
void WritePlotFilesSF_2D(const amrex::MultiFab& mag, const amrex::MultiFab& realimag, const amrex::Geometry& geom,
const int step, const Real time, const amrex::Vector< std::string >& names, std::string plotfile_base);

void EvaluateWritePlotFileVelGrad(int step,
const amrex::Real time,
const amrex::Geometry& geom,
const std::array<MultiFab, AMREX_SPACEDIM>& vel);

void conservedToPrimitiveStag(MultiFab& prim_in, std::array<MultiFab, AMREX_SPACEDIM>& velStag_in,
MultiFab& cons_in, const std::array<MultiFab, AMREX_SPACEDIM>& momStag_in);

Expand Down
9 changes: 8 additions & 1 deletion src_compressible_stag/main_driver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@ void main_driver(const char* argv)
std::array< MultiFab, AMREX_SPACEDIM > velVars;
std::array< MultiFab, AMREX_SPACEDIM > cumomMeans;
std::array< MultiFab, AMREX_SPACEDIM > cumomVars;

if ((plot_cross) and ((cross_cell < 0) or (cross_cell > n_cells[0]-1))) {
Abort("Cross cell needs to be within the domain: 0 <= cross_cell <= n_cells[0] - 1");
}
Expand Down Expand Up @@ -772,6 +772,9 @@ void main_driver(const char* argv)
if (plot_int > 0) {
WritePlotFileStag(0, 0.0, geom, cu, cuMeans, cuVars, cumom, cumomMeans, cumomVars,
prim, primMeans, primVars, vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, eta, kappa);
if (turbForcing > 0) {
EvaluateWritePlotFileVelGrad(0, 0.0, geom, vel);
}

if (plot_cross) {
if (do_1D) {
Expand Down Expand Up @@ -1297,6 +1300,10 @@ void main_driver(const char* argv)
// cuMeansAv, cuVarsAv, primMeansAv, primVarsAv, spatialCrossAv);
WritePlotFileStag(step, time, geom, cu, cuMeans, cuVars, cumom, cumomMeans, cumomVars,
prim, primMeans, primVars, vel, velMeans, velVars, coVars, surfcov, surfcovMeans, surfcovVars, eta, kappa);

if (turbForcing > 0) {
EvaluateWritePlotFileVelGrad(step, time, geom, vel);
}

if (plot_cross) {
if (do_1D) {
Expand Down
122 changes: 122 additions & 0 deletions src_compressible_stag/writePlotFileStag.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -515,3 +515,125 @@ void WritePlotFilesSF_2D(const amrex::MultiFab& mag, const amrex::MultiFab& real

WriteSingleLevelPlotfile(plotfilename2,realimag,varNames,geom2,time,step);
}

void EvaluateWritePlotFileVelGrad(int step,
const amrex::Real time,
const amrex::Geometry& geom,
const std::array<MultiFab, AMREX_SPACEDIM>& vel)
{
BL_PROFILE_VAR("EvaluateWritePlotFileVelGrad()",EvaluateWritePlotFileVelGrad);

// Evaluate velocity gradient components and divergence and vorticity
MultiFab vel_grad;

// Cell-Centered Velocity Gradient Stats (1,2,3 are directions)
// 0: u_1,1
// 1: u_2,2
// 2: u_3,3
// 3: u_1,2
// 4: u_1,3
// 5: u_2,3
// 6: divergence = u_1,1 + u_2,2 + u_3,3
// 7: vorticity = sqrt(wx + wy + wz)
vel_grad.define(convert(vel[0].boxArray(),IntVect(AMREX_D_DECL(0,0,0))), vel[0].DistributionMap(), 8, 0);
vel_grad.setVal(0.0);

const GpuArray<Real, AMREX_SPACEDIM> dx = geom.CellSizeArray();

for ( MFIter mfi(vel_grad,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {

const Box& bx = mfi.tilebox();

const Array4<Real> & vgrad = vel_grad.array(mfi);

AMREX_D_TERM(Array4<Real const> const& velx = vel[0].array(mfi);,
Array4<Real const> const& vely = vel[1].array(mfi);,
Array4<Real const> const& velz = vel[2].array(mfi););

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
// u_1,1
vgrad(i,j,k,0) = (velx(i+1,j,k) - velx(i,j,k))/dx[0];

// u_2,2
vgrad(i,j,k,1) = (vely(i,j+1,k) - vely(i,j,k))/dx[1];

// u_3,3
vgrad(i,j,k,2) = (velz(i,j,k+1) - velz(i,j,k))/dx[2];

// divergence
vgrad(i,j,k,6) = vgrad(i,j,k,0) + vgrad(i,j,k,1) + vgrad(i,j,k,2);

// on edges: u_1,2 and u_2,1 and curl w1 = u_2,1 - u_1,2
Real u12_mm = (velx(i,j,k) - velx(i,j-1,k))/dx[1];
Real u21_mm = (vely(i,j,k) - vely(i-1,j,k))/dx[0];
Real w1_mm = u21_mm - u12_mm;
Real u12_mp = (velx(i,j+1,k) - velx(i,j,k))/dx[1];
Real u21_mp = (vely(i,j+1,k) - vely(i-1,j+1,k))/dx[0];
Real w1_mp = u21_mp - u12_mp;
Real u12_pm = (velx(i+1,j,k) - velx(i+1,j-1,k))/dx[1];
Real u21_pm = (vely(i+1,j,k) - vely(i,j,k))/dx[0];
Real w1_pm = u21_pm - u12_pm;
Real u12_pp = (velx(i+1,j+1,k) - velx(i+1,j,k))/dx[1];
Real u21_pp = (vely(i+1,j+1,k) - vely(i,j+1,k))/dx[0];
Real w1_pp = u21_pp - u12_pp;

// u_1,2
vgrad(i,j,k,3) = 0.25*(u12_mm + u12_mp + u12_pm + u12_pp);

// on edges: u_1,3 and u_3,1 and curl w2 = u_1,3 - u_3,1
Real u13_mm = (velx(i,j,k) - velx(i,j,k-1))/dx[2];
Real u31_mm = (velz(i,j,k) - velz(i-1,j,k))/dx[0];
Real w2_mm = u13_mm - u31_mm;
Real u13_mp = (velx(i,j,k+1) - velx(i,j,k))/dx[2];
Real u31_mp = (velz(i,j,k+1) - velz(i-1,j,k+1))/dx[0];
Real w2_mp = u13_mp - u31_mp;
Real u13_pm = (velx(i+1,j,k) - velx(i+1,j,k-1))/dx[2];
Real u31_pm = (velz(i+1,j,k) - velz(i,j,k))/dx[0];
Real w2_pm = u13_pm - u31_pm;
Real u13_pp = (velx(i+1,j,k+1) - velx(i+1,j,k))/dx[2];
Real u31_pp = (velz(i+1,j,k+1) - velz(i,j,k+1))/dx[0];
Real w2_pp = u13_pp - u31_pp;

// u_1,3
vgrad(i,j,k,4) = 0.25*(u13_mm + u13_mp + u13_pm + u13_pp);

// on edges: u_2,3 and u_3,2 and curl w2 = u_3,2 - u_2,3
Real u23_mm = (vely(i,j,k) - vely(i,j,k-1))/dx[2];
Real u32_mm = (velz(i,j,k) - velz(i,j-1,k))/dx[1];
Real w3_mm = u32_mm - u23_mm;
Real u23_mp = (vely(i,j,k+1) - vely(i,j,k))/dx[2];
Real u32_mp = (velz(i,j,k+1) - velz(i,j-1,k+1))/dx[1];
Real w3_mp = u32_mp - u23_mp;
Real u23_pm = (vely(i,j+1,k) - vely(i,j+1,k-1))/dx[2];
Real u32_pm = (velz(i,j+1,k) - velz(i,j,k))/dx[1];
Real w3_pm = u32_pm - u23_pm;
Real u23_pp = (vely(i,j+1,k+1) - vely(i,j+1,k))/dx[2];
Real u32_pp = (velz(i,j+1,k+1) - velz(i,j,k+1))/dx[1];
Real w3_pp = u32_pp - u23_pp;

// u_2,3
vgrad(i,j,k,5) = 0.25*(u23_mm + u23_mp + u23_pm + u23_pp);

// vorticity magnitude: sqrt(w1*w1 + w2*w2 + w3*w3)
vgrad(i,j,k,7) = sqrt(0.25*(w1_mm*w1_mm + w1_mp*w1_mp + w1_pm*w1_pm + w1_pp*w1_pp +
w2_mm*w2_mm + w2_mp*w2_mp + w2_pm*w2_pm + w2_pp*w2_pp +
w3_mm*w3_mm + w3_mp*w3_mp + w3_pm*w3_pm + w3_pp*w3_pp));
});
}

// Write on a plotfile
std::string plotfilename = amrex::Concatenate("vel_grad",step,9);
amrex::Vector<std::string> varNames(8);
varNames[0] = "ux_x";
varNames[1] = "uy_y";
varNames[2] = "uz_z";
varNames[3] = "ux_y";
varNames[4] = "ux_z";
varNames[5] = "uy_z";
varNames[6] = "div";
varNames[7] = "vort";
WriteSingleLevelPlotfile(plotfilename,vel_grad,varNames,geom,time,step);
}


0 comments on commit 9e054f5

Please sign in to comment.