Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/main' into bulk_visc
Browse files Browse the repository at this point in the history
  • Loading branch information
@isriva
isriva committed Sep 22, 2023
2 parents 1fd3499 + 46f74a3 commit 9271cbf
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 47 deletions.
124 changes: 79 additions & 45 deletions exec/thinFilm/main_driver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,26 +182,20 @@ void main_driver(const char* argv)
// initialize height
height.setVal(thinfilm_h0);

Real time = 0.;

// update height using forward Euler
if(thinfilm_pertamp > 0.){
for ( MFIter mfi(height,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
if (thinfilm_pertamp > 0.) {
for ( MFIter mfi(height,TilingIfNotGPU()); mfi.isValid(); ++mfi ) {

const Box& bx = mfi.tilebox();

const Array4<Real> & h = height.array(mfi);
const Box& bx = mfi.tilebox();
const Array4<Real> & h = height.array(mfi);

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
h(i,j,k) += thinfilm_pertamp * thinfilm_h0*std::cos(2.*3.14159265358979323846264338327950288*(i+.5)*dx[0]/(prob_hi[0]-prob_lo[0]));
if(i==0 && j ==0){
amrex::Print() << " HEIGHT " << time << " " << h(i,j,k)-thinfilm_h0 << std::endl;
}
});
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
h(i,j,k) += thinfilm_pertamp * thinfilm_h0*std::cos(2.*M_PI*(i+.5)*dx[0]/(prob_hi[0]-prob_lo[0]));
});

}
}
}
}

height.FillBoundary(geom.periodicity());

Expand All @@ -211,7 +205,9 @@ void main_driver(const char* argv)
// constant factor in noise term
Real ConstNoise = 2.*k_B*T_init[0] / (3.*visc_coef);
Real Const3dx = thinfilm_gamma / (3.*visc_coef);


Real time = 0.;

Real dx_min;
if (algorithm_type == 0) {
dx_min = dx[0];
Expand All @@ -222,36 +218,50 @@ void main_driver(const char* argv)
}

//////////////////////////////////////
// data structures to take 1D ffts for 1D mode
// for ASCII writeout and 1D ffts in 1D mode
BoxArray ba_onegrid(domain);
DistributionMapping dmap_onegrid(ba_onegrid);

MultiFab height_onegrid(ba_onegrid, dmap_onegrid, 1, 0);

// make the "1D" MultiFab x-based, regardless of whether problem is x or y-based
IntVect dom_lo_flat(0,0);
IntVect dom_hi_flat(n_cells[algorithm_type]-1,0);
// data structures to take 1D ffts in 1D mode
MultiFab height_flat_onegrid;
MultiFab fft_sum;
MultiFab fft_avg;

Geometry geom_fft;

long npts_flat;

Box domain_flat(dom_lo_flat,dom_hi_flat);
BoxArray ba_flat_onegrid(domain_flat);
if (algorithm_type == 0 || algorithm_type == 1) {
// make the "1D" MultiFab x-based, regardless of whether problem is x or y-based
IntVect dom_lo_flat(0,0);
IntVect dom_hi_flat(n_cells[algorithm_type]-1,0);

// use same dmap as height_onegrid so we can easily copy a strip into here
MultiFab height_flat_onegrid(ba_flat_onegrid,dmap_onegrid,1,0);
Box domain_flat(dom_lo_flat,dom_hi_flat);
BoxArray ba_flat_onegrid(domain_flat);

// fft contains n_cell/2+1 points (index space 0 to n_cell/2
IntVect dom_lo_fft(0,0);
IntVect dom_hi_fft(n_cells[algorithm_type]/2,0);
Box domain_fft(dom_lo_fft,dom_hi_fft);
BoxArray ba_fft(domain_fft);
npts_flat = domain_flat.numPts();

// use same dmap as height_onegrid so we can easily copy a strip into here
height_flat_onegrid.define(ba_flat_onegrid,dmap_onegrid,1,0);

// fft contains n_cell/2+1 points (index space 0 to n_cell/2
IntVect dom_lo_fft(0,0);
IntVect dom_hi_fft(n_cells[algorithm_type]/2,0);
Box domain_fft(dom_lo_fft,dom_hi_fft);
BoxArray ba_fft(domain_fft);

// magnitude of fft, sum and average
// use same dmap as height_onegrid
MultiFab fft_sum(ba_fft,dmap_onegrid,1,0);
MultiFab fft_avg(ba_fft,dmap_onegrid,1,0);
// magnitude of fft, sum and average
// use same dmap as height_onegrid
fft_sum.define(ba_fft,dmap_onegrid,1,0);
fft_avg.define(ba_fft,dmap_onegrid,1,0);

fft_sum.setVal(0.);
fft_sum.setVal(0.);

Geometry geom_fft(domain_fft,&real_box,CoordSys::cartesian,is_periodic.data());
geom_fft.define(domain_fft,&real_box,CoordSys::cartesian,is_periodic.data());
}

//////////////////////////////////////

// time step
Expand Down Expand Up @@ -298,13 +308,11 @@ void main_driver(const char* argv)
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
hfacex(i,j,k) = 0.5*( h(i-1,j,k) + h(i,j,k) );
hfacex(i,j,k) = thinfilm_h0;
gradhx(i,j,k) = ( h(i,j,k) - h(i-1,j,k) ) / dx[0];
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
hfacey(i,j,k) = 0.5*( h(i,j-1,k) + h(i,j,k) );
hfacey(i,j,k) = thinfilm_h0;
gradhy(i,j,k) = ( h(i,j,k) - h(i,j-1,k) ) / dx[1];
});

Expand Down Expand Up @@ -497,6 +505,37 @@ void main_driver(const char* argv)
{
const std::string& pltfile = amrex::Concatenate("plt",istep,10);
WriteSingleLevelPlotfile(pltfile, height, {"height"}, geom, time, 0);

// ASCII

// copy distributed data into 1D data
height_onegrid.ParallelCopy(height, 0, 0, 1);

for ( MFIter mfi(height_onegrid,false); mfi.isValid(); ++mfi ) {

std::ofstream hstream;
std::string heightBaseName = "height";

std::string heightName = Concatenate(heightBaseName,istep,10);
heightName += ".txt";

hstream.open(heightName);

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

const Array4<Real>& mfdata = height_onegrid.array(mfi);

for (auto j = lo.y; j <= hi.y; ++j) {
for (auto i = lo.x; i <= hi.x; ++i) {
hstream << std::setprecision(15) << mfdata(i,j,0) << " ";
}
hstream << "\n";
}

} // end MFIter

}

// statistics
Expand Down Expand Up @@ -547,10 +586,6 @@ void main_driver(const char* argv)
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
dhstar(i,j,k) += ( h(thinfilm_icorr,j,k)-thinfilm_h0 ) * ( h(i,j,k)-thinfilm_h0 );
// HACK
if(i==0 && j ==0 && thinfilm_pertamp > 0.){
amrex::Print() << " HEIGHT " << time << " " << h(i,j,k)-thinfilm_h0 << std::endl;
}
});
} else {
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
Expand Down Expand Up @@ -648,8 +683,7 @@ void main_driver(const char* argv)
Vector<std::unique_ptr<BaseFab<GpuComplex<Real> > > > spectral_field;
Vector<FFTplan> forward_plan;

long npts = domain_flat.numPts();
Real sqrtnpts = std::sqrt(npts);
Real sqrtnpts_flat = std::sqrt(npts_flat);

// take fft of strip and add magnitude of result to fft_sum
for (MFIter mfi(height_flat_onegrid); mfi.isValid(); ++mfi) {
Expand Down Expand Up @@ -725,7 +759,7 @@ void main_driver(const char* argv)

ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
sum(i,j,k) += ( spectral(i,j,k).real()*spectral(i,j,k).real() + spectral(i,j,k).imag()*spectral(i,j,k).imag() ) / sqrtnpts;
sum(i,j,k) += ( spectral(i,j,k).real()*spectral(i,j,k).real() + spectral(i,j,k).imag()*spectral(i,j,k).imag() ) / sqrtnpts_flat;
if (i == 0) sum(i,j,k) = 0.;
});
}
Expand Down
5 changes: 3 additions & 2 deletions exec/thinFilm/thinfilm_functions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,16 @@ void InitializeThinfilmNamespace() {
// defaults
thinfilm_icorr = 0;
thinfilm_jcorr = 0;
thinfilm_pertamp = 0;
thinfilm_pertamp = 0.;

do_fft_diag = 1;

ParmParse pp;

pp.get("thinfilm_h0",thinfilm_h0);
pp.get("thinfilm_gamma",thinfilm_gamma);
pp.get("thinfilm_pertamp",thinfilm_pertamp);

pp.query("thinfilm_pertamp",thinfilm_pertamp);

pp.query("thinfilm_icorr",thinfilm_icorr);
pp.query("thinfilm_jcorr",thinfilm_jcorr);
Expand Down

0 comments on commit 9271cbf

Please sign in to comment.