Skip to content

Commit

Permalink
some more changes to the reservoir code
Browse files Browse the repository at this point in the history
  • Loading branch information
@isriva
isriva committed Aug 9, 2023
1 parent d2a6dae commit af236b7
Show file tree
Hide file tree
Showing 2 changed files with 120 additions and 96 deletions.
146 changes: 73 additions & 73 deletions src_compressible_stag/main_driver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,12 +248,12 @@ void main_driver(const char* argv)
std::string filename = "crossMeans";
std::ofstream outfile;

std::string filenamereslo = "consResLo";
std::ofstream outfilereslo;
if (ParallelDescriptor::IOProcessor()) outfilereslo.open(filenamereslo);
std::string filenamereshi = "consResHi";
std::ofstream outfilereshi;
if (ParallelDescriptor::IOProcessor()) outfilereshi.open(filenamereshi);
// std::string filenamereslo = "consResLo";
// std::ofstream outfilereslo;
// if (ParallelDescriptor::IOProcessor()) outfilereslo.open(filenamereslo);
// std::string filenamereshi = "consResHi";
// std::ofstream outfilereshi;
// if (ParallelDescriptor::IOProcessor()) outfilereshi.open(filenamereshi);
/////////////////////////////////////////////
// Setup Structure factor variables & scaling
/////////////////////////////////////////////
Expand Down Expand Up @@ -920,7 +920,7 @@ void main_driver(const char* argv)
// timer
Real ts2 = ParallelDescriptor::second() - ts1;
ParallelDescriptor::ReduceRealMax(ts2, ParallelDescriptor::IOProcessorNumber());
if (step%1 == 0) {
if (step%100 == 0) {
amrex::Print() << "Advanced step " << step << " in " << ts2 << " seconds\n";
}

Expand Down Expand Up @@ -988,7 +988,7 @@ void main_driver(const char* argv)
statsCount, geom);
}
statsCount++;
if (step%1 == 0) {
if (step%100 == 0) {
amrex::Print() << "Mean Density: " << ComputeSpatialMean(cu, 0) << " Mean Momentum (x):" << ComputeSpatialMean(cumom[0], 0) << " Mean Energy:" << ComputeSpatialMean(cu, 4) << "\n";
}

Expand Down Expand Up @@ -1249,72 +1249,72 @@ void main_driver(const char* argv)
}
}

// reservoir debug output
Box dom(geom.Domain());

for ( MFIter mfi(cu); mfi.isValid(); ++mfi) {

const Box& bx = mfi.validbox();
const auto lo = amrex::lbound(bx);
const auto hi = amrex::ubound(bx);

AMREX_D_TERM(Array4<Real const> const& xmom = cumom[0].array(mfi);,
Array4<Real const> const& ymom = cumom[1].array(mfi);,
Array4<Real const> const& zmom = cumom[2].array(mfi););

const Array4<const Real>& cufab = cu.array(mfi);
const Array4<const Real>& primfab = prim.array(mfi);

// Reservoir in LO X
if ((bc_mass_lo[0] == 4) and (bx.smallEnd(0) <= dom.smallEnd(0))) {
for (auto k = lo.z; k <= hi.z; ++k) {
for (auto j = lo.y; j <= hi.y; ++j) {
for (auto i = lo.x; i <= hi.x; ++i) {
if (i == dom.smallEnd(0) and (j==0) and (k==0)) {
Real dens = cufab(i,j,k,0);
Real en = cufab(i,j,k,4);
Real rho1 = cufab(i,j,k,5);
Real rho2 = cufab(i,j,k,6);
Real temp = primfab(i,j,k,4);
Real press = primfab(i,j,k,5);
Real moment = xmom(i,j,k);
if (ParallelDescriptor::IOProcessor()) {
outfilereslo << dens << " " << en << " " << rho1 << " " << rho2 << " " << temp << " " << press << " " << moment << std::endl;
}
}
}
}
}
}

// Reservoir in HI X
if ((bc_mass_hi[0] == 4) and (bx.bigEnd(0) >= dom.bigEnd(0))) {
for (auto k = lo.z; k <= hi.z; ++k) {
for (auto j = lo.y; j <= hi.y; ++j) {
for (auto i = lo.x; i <= hi.x; ++i) {
if (i == dom.bigEnd(0) and (j==0) and (k==0)) {
Real dens = cufab(i,j,k,0);
Real en = cufab(i,j,k,4);
Real rho1 = cufab(i,j,k,5);
Real rho2 = cufab(i,j,k,6);
Real temp = primfab(i,j,k,4);
Real press = primfab(i,j,k,5);
Real moment = xmom(i+1,j,k);
if (ParallelDescriptor::IOProcessor()) {
outfilereshi << dens << " " << en << " " << rho1 << " " << rho2 << " " << temp << " " << press << " " << moment << std::endl;
}
}
}
}
}
}
}
// // reservoir debug output
// Box dom(geom.Domain());
//
// for ( MFIter mfi(cu); mfi.isValid(); ++mfi) {
//
// const Box& bx = mfi.validbox();
// const auto lo = amrex::lbound(bx);
// const auto hi = amrex::ubound(bx);
//
// AMREX_D_TERM(Array4<Real const> const& xmom = cumom[0].array(mfi);,
// Array4<Real const> const& ymom = cumom[1].array(mfi);,
// Array4<Real const> const& zmom = cumom[2].array(mfi););
//
// const Array4<const Real>& cufab = cu.array(mfi);
// const Array4<const Real>& primfab = prim.array(mfi);
//
// // Reservoir in LO X
// if ((bc_mass_lo[0] == 4) and (bx.smallEnd(0) <= dom.smallEnd(0))) {
// for (auto k = lo.z; k <= hi.z; ++k) {
// for (auto j = lo.y; j <= hi.y; ++j) {
// for (auto i = lo.x; i <= hi.x; ++i) {
// if (i == dom.smallEnd(0) and (j==0) and (k==0)) {
// Real dens = cufab(i,j,k,0);
// Real en = cufab(i,j,k,4);
// Real rho1 = cufab(i,j,k,5);
// Real rho2 = cufab(i,j,k,6);
// Real temp = primfab(i,j,k,4);
// Real press = primfab(i,j,k,5);
// Real moment = xmom(i,j,k);
// if (ParallelDescriptor::IOProcessor()) {
// outfilereslo << dens << " " << en << " " << rho1 << " " << rho2 << " " << temp << " " << press << " " << moment << std::endl;
// }
// }
// }
// }
// }
// }
//
// // Reservoir in HI X
// if ((bc_mass_hi[0] == 4) and (bx.bigEnd(0) >= dom.bigEnd(0))) {
// for (auto k = lo.z; k <= hi.z; ++k) {
// for (auto j = lo.y; j <= hi.y; ++j) {
// for (auto i = lo.x; i <= hi.x; ++i) {
// if (i == dom.bigEnd(0) and (j==0) and (k==0)) {
// Real dens = cufab(i,j,k,0);
// Real en = cufab(i,j,k,4);
// Real rho1 = cufab(i,j,k,5);
// Real rho2 = cufab(i,j,k,6);
// Real temp = primfab(i,j,k,4);
// Real press = primfab(i,j,k,5);
// Real moment = xmom(i+1,j,k);
// if (ParallelDescriptor::IOProcessor()) {
// outfilereshi << dens << " " << en << " " << rho1 << " " << rho2 << " " << temp << " " << press << " " << moment << std::endl;
// }
// }
// }
// }
// }
// }
// }


// timer
Real aux2 = ParallelDescriptor::second() - aux1;
ParallelDescriptor::ReduceRealMax(aux2, ParallelDescriptor::IOProcessorNumber());
if (step%1 == 0) {
if (step%100 == 0) {
amrex::Print() << "Aux time (stats, struct fac, plotfiles) " << aux2 << " seconds\n";
}

Expand All @@ -1329,7 +1329,7 @@ void main_driver(const char* argv)
ParallelDescriptor::ReduceLongMin(min_fab_megabytes, IOProc);
ParallelDescriptor::ReduceLongMax(max_fab_megabytes, IOProc);

if (step%1 == 0) {
if (step%100 == 0) {
amrex::Print() << "High-water FAB megabyte spread across MPI nodes: ["
<< min_fab_megabytes << " ... " << max_fab_megabytes << "]\n";
}
Expand All @@ -1340,15 +1340,15 @@ void main_driver(const char* argv)
ParallelDescriptor::ReduceLongMin(min_fab_megabytes, IOProc);
ParallelDescriptor::ReduceLongMax(max_fab_megabytes, IOProc);

if (step%1 == 0) {
if (step%100 == 0) {
amrex::Print() << "Curent FAB megabyte spread across MPI nodes: ["
<< min_fab_megabytes << " ... " << max_fab_megabytes << "]\n";
}
}

if (ParallelDescriptor::IOProcessor()) outfile.close();
if (ParallelDescriptor::IOProcessor()) outfilereslo.close();
if (ParallelDescriptor::IOProcessor()) outfilereshi.close();
// if (ParallelDescriptor::IOProcessor()) outfilereslo.close();
// if (ParallelDescriptor::IOProcessor()) outfilereshi.close();

// timer
Real stop_time = ParallelDescriptor::second() - strt_time;
Expand Down
70 changes: 47 additions & 23 deletions src_compressible_stag/timeStepStag.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,11 +311,13 @@ void RK3stepStag(MultiFab& cu,
geom, stoch_weights,dt);

// add to the total continuum fluxes based on RK3 weight
Real aux1 = ParallelDescriptor::second();
if (do_reservoir) {
for (int d=0;d<AMREX_SPACEDIM;++d) {
MultiFab::Saxpy(faceflux_cont[d],Real(1.0/6.0),faceflux[d],0,0,nvars,0);
}
ComputeFluxMomReservoir(cu,prim,vel,cumom_res,faceflux_res,geom,dt); // compute fluxes and momentum from reservoir particle update
ResetReservoirFluxes(faceflux, faceflux_res, geom); // reset fluxes in the FHD-reservoir interface from particle update
}
Real aux2 = ParallelDescriptor::second() - aux1;
ParallelDescriptor::ReduceRealMax(aux2, ParallelDescriptor::IOProcessorNumber());

if (nreaction>0) {
MultiFab::LinComb(ranchem,
Expand Down Expand Up @@ -409,6 +411,10 @@ void RK3stepStag(MultiFab& cu,
BCMomTrans(cupmom[i], vel[i], geom, i);
}

if (do_reservoir) {
ResetReservoirMom(cupmom, cumom_res, geom); // set momentum at the reservoir interface to its value from particle update
}

// Fill boundaries for conserved variables
for (int d=0; d<AMREX_SPACEDIM; d++) {
cupmom[d].FillBoundary(geom.periodicity());
Expand Down Expand Up @@ -485,11 +491,13 @@ void RK3stepStag(MultiFab& cu,
geom, stoch_weights,dt);

// add to the total continuum fluxes based on RK3 weight
Real aux3 = ParallelDescriptor::second();
if (do_reservoir) {
for (int d=0;d<AMREX_SPACEDIM;++d) {
MultiFab::Saxpy(faceflux_cont[d],Real(1.0/6.0),faceflux[d],0,0,nvars,0);
}
ComputeFluxMomReservoir(cup,prim,vel,cumom_res,faceflux_res,geom,0.25*dt); // compute fluxes and momentum from reservoir particle update
ResetReservoirFluxes(faceflux, faceflux_res, geom); // reset fluxes in the FHD-reservoir interface from particle update
}
Real aux4 = ParallelDescriptor::second() - aux3;
ParallelDescriptor::ReduceRealMax(aux4, ParallelDescriptor::IOProcessorNumber());

if (nreaction>0) {
MultiFab::LinComb(ranchem,
Expand Down Expand Up @@ -588,6 +596,10 @@ void RK3stepStag(MultiFab& cu,
BCMomTrans(cup2mom[i], vel[i], geom, i);
}

if (do_reservoir) {
ResetReservoirMom(cup2mom, cumom_res, geom); // set momentum at the reservoir interface to its value from particle update
}

// Fill boundaries for conserved variables
for (int d=0; d<AMREX_SPACEDIM; d++) {
cup2mom[d].FillBoundary(geom.periodicity());
Expand Down Expand Up @@ -664,11 +676,13 @@ void RK3stepStag(MultiFab& cu,
geom, stoch_weights,dt);

// add to the total continuum fluxes based on RK3 weight
Real aux5 = ParallelDescriptor::second();
if (do_reservoir) {
for (int d=0;d<AMREX_SPACEDIM;++d) {
MultiFab::Saxpy(faceflux_cont[d],Real(2.0/3.0),faceflux[d],0,0,nvars,0);
}
ComputeFluxMomReservoir(cup2,prim,vel,cumom_res,faceflux_res,geom,(2.0/3.0)*dt); // compute fluxes and momentum from reservoir particle update
ResetReservoirFluxes(faceflux, faceflux_res, geom); // reset fluxes in the FHD-reservoir interface from particle update
}
Real aux6 = ParallelDescriptor::second() - aux5;
ParallelDescriptor::ReduceRealMax(aux6, ParallelDescriptor::IOProcessorNumber());

if (nreaction>0) {
MultiFab::LinComb(ranchem,
Expand Down Expand Up @@ -762,23 +776,27 @@ void RK3stepStag(MultiFab& cu,
BCMomNormal(cumom[i], vel[i], cu, geom, i);
BCMomTrans(cumom[i], vel[i], geom, i);
}

/////////////////////////////////////////////////////
// compute fluxes and momentum at the reservoir /////
// over an intergration step dt /////////////////////
// timer
Real aux1 = ParallelDescriptor::second();

if (do_reservoir) {
ComputeFluxMomReservoir(cu0,prim0,vel0,cumom_res,faceflux_res,geom,dt); // compute fluxes and momentum from reservoir particle update
ResetReservoirMom(cumom, cumom_res, geom); // set momentum at the reservoir interface to its value from particle update
ReFluxCons(cu, cu0, faceflux_res, faceflux_cont, geom, dt); // reflux the conservative qtys in the adjacent cell from reservoir flux
}
Real aux2 = ParallelDescriptor::second() - aux1;
ParallelDescriptor::ReduceRealMax(aux2, ParallelDescriptor::IOProcessorNumber());
if (step%1 == 0) {
amrex::Print() << "Step: " << step << " Reservoir generator time: " << aux2 << " seconds\n";
}
/////////////////////////////////////////////////////

// /////////////////////////////////////////////////////
// // compute fluxes and momentum at the reservoir /////
// // over an intergration step dt /////////////////////
// // timer
// Real aux1 = ParallelDescriptor::second();
// if (do_reservoir) {
// ComputeFluxMomReservoir(cu0,prim0,vel0,cumom_res,faceflux_res,geom,dt); // compute fluxes and momentum from reservoir particle update
// ResetReservoirMom(cumom, cumom_res, geom); // set momentum at the reservoir interface to its value from particle update
// ReFluxCons(cu, cu0, faceflux_res, faceflux_cont, geom, dt); // reflux the conservative qtys in the adjacent cell from reservoir flux
// }
// Real aux2 = ParallelDescriptor::second() - aux1;
// ParallelDescriptor::ReduceRealMax(aux2, ParallelDescriptor::IOProcessorNumber());
// if (step%1 == 0) {
// amrex::Print() << "Step: " << step << " Reservoir generator time: " << aux2 << " seconds\n";
// }
// /////////////////////////////////////////////////////

// Fill boundaries for conserved variables
for (int d=0; d<AMREX_SPACEDIM; d++) {
Expand All @@ -803,4 +821,10 @@ void RK3stepStag(MultiFab& cu,

// Correctly set momentum and velocity at the walls & temperature, pressure, density & mass/mole fractions in ghost cells
setBCStag(prim, cu, cumom, vel, geom);

if (do_reservoir) {
if (step%100 == 0) {
amrex::Print() << "Step: " << step << " Reservoir generator time: " << aux2 + aux4 + aux6 << " seconds\n";
}
}
}

0 comments on commit af236b7

Please sign in to comment.