diff --git a/exec/thinFilm/main_driver.cpp b/exec/thinFilm/main_driver.cpp index df115084..997c1783 100644 --- a/exec/thinFilm/main_driver.cpp +++ b/exec/thinFilm/main_driver.cpp @@ -137,7 +137,9 @@ void main_driver(const char* argv) MultiFab height(ba, dmap, 1, 1); MultiFab Laph (ba, dmap, 1, 1); + MultiFab disjoining(ba, dmap, 1, 1); Laph.setVal(0.); // prevent intermediate NaN calculations behind physical boundaries + disjoining.setVal(0.); // for statsitics @@ -169,6 +171,11 @@ void main_driver(const char* argv) gradLaph[1].define(convert(ba,nodal_flag_y), dmap, 1, 0);, gradLaph[2].define(convert(ba,nodal_flag_z), dmap, 1, 0);); + std::array< MultiFab, AMREX_SPACEDIM > gradDisjoining; + AMREX_D_TERM(gradDisjoining[0].define(convert(ba,nodal_flag_x), dmap, 1, 0);, + gradDisjoining[1].define(convert(ba,nodal_flag_y), dmap, 1, 0);, + gradDisjoining[2].define(convert(ba,nodal_flag_z), dmap, 1, 0);); + std::array< MultiFab, AMREX_SPACEDIM > flux; AMREX_D_TERM(flux[0] .define(convert(ba,nodal_flag_x), dmap, 1, 0);, flux[1] .define(convert(ba,nodal_flag_y), dmap, 1, 0);, @@ -205,6 +212,7 @@ 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 Const3dx_nogamma = 1. / (3.*visc_coef); Real time = 0.; @@ -365,12 +373,16 @@ void main_driver(const char* argv) } - // compute Laph + // compute Laph and disjoining for ( MFIter mfi(Laph,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { const Box& bx = mfi.tilebox(); const Array4 & L = Laph.array(mfi); + + const Array4 & h = height.array(mfi); + + const Array4 & Disjoining = disjoining.array(mfi); AMREX_D_TERM(const Array4 & gradhx = gradh[0].array(mfi);, const Array4 & gradhy = gradh[1].array(mfi);, @@ -380,11 +392,14 @@ void main_driver(const char* argv) { L(i,j,k) = x_flux_fac * (gradhx(i+1,j,k) - gradhx(i,j,k)) / dx[0] + y_flux_fac * (gradhy(i,j+1,k) - gradhy(i,j,k)) / dx[1]; + + Disjoining(i,j,k) = thinfilm_hamaker / (6.*M_PI*std::pow(h(i,j,k),3.)); }); } Laph.FillBoundary(geom.periodicity()); + disjoining.FillBoundary(geom.periodicity()); - // compute gradLaph + // compute gradLaph and gradDisjoining for ( MFIter mfi(height,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { AMREX_D_TERM(const Box & bx_x = mfi.nodaltilebox(0);, @@ -394,17 +409,25 @@ void main_driver(const char* argv) AMREX_D_TERM(const Array4 & gradLaphx = gradLaph[0].array(mfi);, const Array4 & gradLaphy = gradLaph[1].array(mfi);, const Array4 & gradLaphz = gradLaph[2].array(mfi);); + + AMREX_D_TERM(const Array4 & gradDisjoiningx = gradDisjoining[0].array(mfi);, + const Array4 & gradDisjoiningy = gradDisjoining[1].array(mfi);, + const Array4 & gradDisjoiningz = gradDisjoining[2].array(mfi);); const Array4 & L = Laph.array(mfi); + const Array4 & Disjoining = disjoining.array(mfi); + amrex::ParallelFor(bx_x, bx_y, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { gradLaphx(i,j,k) = ( L(i,j,k) - L(i-1,j,k) ) / dx[0]; + gradDisjoiningx(i,j,k) = ( Disjoining(i,j,k) -Disjoining(i-1,j,k) ) / dx[0]; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { gradLaphy(i,j,k) = ( L(i,j,k) - L(i,j-1,k) ) / dx[1]; + gradDisjoiningy(i,j,k) = ( Disjoining(i,j,k) -Disjoining(i,j-1,k) ) / dx[1]; }); } @@ -431,18 +454,24 @@ void main_driver(const char* argv) const Array4 & randfacey = randface[1].array(mfi);, const Array4 & randfacez = randface[2].array(mfi);); + AMREX_D_TERM(const Array4 & gradDisjoiningx = gradDisjoining[0].array(mfi);, + const Array4 & gradDisjoiningy = gradDisjoining[1].array(mfi);, + const Array4 & gradDisjoiningz = gradDisjoining[2].array(mfi);); + amrex::ParallelFor(bx_x, bx_y, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { fluxx(i,j,k) = x_flux_fac * ( std::sqrt(ConstNoise*std::pow(hfacex(i,j,k),3.) / (dt*dVol)) * randfacex(i,j,k) - + Const3dx * std::pow(hfacex(i,j,k),3.)*gradLaphx(i,j,k) ); + + Const3dx * std::pow(hfacex(i,j,k),3.)*gradLaphx(i,j,k) + + Const3dx_nogamma * std::pow(hfacex(i,j,k),3.)*gradDisjoiningx(i,j,k)); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { fluxy(i,j,k) = y_flux_fac * ( std::sqrt(ConstNoise*std::pow(hfacey(i,j,k),3.) / (dt*dVol)) * randfacey(i,j,k) - + Const3dx * std::pow(hfacey(i,j,k),3.)*gradLaphy(i,j,k) ); + + Const3dx * std::pow(hfacey(i,j,k),3.)*gradLaphy(i,j,k) + + Const3dx_nogamma * std::pow(hfacex(i,j,k),3.)*gradDisjoiningy(i,j,k) ); }); // lo x-faces @@ -489,8 +518,11 @@ void main_driver(const char* argv) const Array4 & h = height.array(mfi); + // amrex::Print{} << "HEIGHT " << time << " " << h(0,0,0) << " " << h(31,0,0) << std::endl; + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + h(i,j,k) -= dt * ( (fluxx(i+1,j,k) - fluxx(i,j,k)) / dx[0] +(fluxy(i,j+1,k) - fluxy(i,j,k)) / dx[1]); }); @@ -511,6 +543,8 @@ void main_driver(const char* argv) // copy distributed data into 1D data height_onegrid.ParallelCopy(height, 0, 0, 1); + amrex::Real sumh = 0.; + for ( MFIter mfi(height_onegrid,false); mfi.isValid(); ++mfi ) { std::ofstream hstream; @@ -530,9 +564,14 @@ void main_driver(const char* argv) 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) << " "; + sumh += mfdata(i,j,0); + if(j==0 && i==0){ + amrex::Print{} << "HEIGHT " << time << " " << mfdata(0,0,0) << " " << mfdata(31,0,0) << std::endl; + } } hstream << "\n"; } + amrex::Print{} << "SUM " << sumh << std::endl; } // end MFIter diff --git a/exec/thinFilm/thinfilm_functions.cpp b/exec/thinFilm/thinfilm_functions.cpp index 3a4917d2..dde70c0b 100644 --- a/exec/thinFilm/thinfilm_functions.cpp +++ b/exec/thinFilm/thinfilm_functions.cpp @@ -5,6 +5,7 @@ AMREX_GPU_MANAGED amrex::Real thinfilm::thinfilm_h0; AMREX_GPU_MANAGED amrex::Real thinfilm::thinfilm_gamma; AMREX_GPU_MANAGED amrex::Real thinfilm::thinfilm_pertamp; +AMREX_GPU_MANAGED amrex::Real thinfilm::thinfilm_hamaker; AMREX_GPU_MANAGED int thinfilm::thinfilm_icorr; AMREX_GPU_MANAGED int thinfilm::thinfilm_jcorr; @@ -19,6 +20,7 @@ void InitializeThinfilmNamespace() { thinfilm_icorr = 0; thinfilm_jcorr = 0; thinfilm_pertamp = 0.; + thinfilm_hamaker = 0.; do_fft_diag = 1; @@ -28,6 +30,7 @@ void InitializeThinfilmNamespace() { pp.get("thinfilm_gamma",thinfilm_gamma); pp.query("thinfilm_pertamp",thinfilm_pertamp); + pp.query("thinfilm_hamaker",thinfilm_hamaker); pp.query("thinfilm_icorr",thinfilm_icorr); pp.query("thinfilm_jcorr",thinfilm_jcorr); diff --git a/exec/thinFilm/thinfilm_namespace.H b/exec/thinFilm/thinfilm_namespace.H index 59091542..16f3afef 100644 --- a/exec/thinFilm/thinfilm_namespace.H +++ b/exec/thinFilm/thinfilm_namespace.H @@ -3,6 +3,7 @@ namespace thinfilm { extern AMREX_GPU_MANAGED amrex::Real thinfilm_h0; extern AMREX_GPU_MANAGED amrex::Real thinfilm_gamma; extern AMREX_GPU_MANAGED amrex::Real thinfilm_pertamp; + extern AMREX_GPU_MANAGED amrex::Real thinfilm_hamaker; extern AMREX_GPU_MANAGED int thinfilm_icorr; extern AMREX_GPU_MANAGED int thinfilm_jcorr;