From 2a39ddb9f125bcddf71bc9c9d6f1dcba039fa298 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 6 Nov 2023 17:10:39 -0800 Subject: [PATCH 1/5] disjoining pressure --- exec/thinFilm/main_driver.cpp | 42 +++++++++++++++++++++++----- exec/thinFilm/thinfilm_functions.cpp | 4 ++- exec/thinFilm/thinfilm_namespace.H | 1 + 3 files changed, 39 insertions(+), 8 deletions(-) diff --git a/exec/thinFilm/main_driver.cpp b/exec/thinFilm/main_driver.cpp index df115084..614cc4ce 100644 --- a/exec/thinFilm/main_driver.cpp +++ b/exec/thinFilm/main_driver.cpp @@ -135,9 +135,11 @@ void main_driver(const char* argv) // define DistributionMapping dmap.define(ba); - MultiFab height(ba, dmap, 1, 1); - MultiFab Laph (ba, dmap, 1, 1); + 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.; @@ -301,7 +309,7 @@ void main_driver(const char* argv) AMREX_D_TERM(const Array4 & gradhx = gradh[0].array(mfi);, const Array4 & gradhy = gradh[1].array(mfi);, const Array4 & gradhz = gradh[2].array(mfi);); - + const Array4 & h = height.array(mfi); amrex::ParallelFor(bx_x, bx_y, @@ -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,13 @@ 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()); - // 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 +408,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 +453,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 diff --git a/exec/thinFilm/thinfilm_functions.cpp b/exec/thinFilm/thinfilm_functions.cpp index 3a4917d2..2c466f32 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; @@ -26,8 +28,8 @@ void InitializeThinfilmNamespace() { pp.get("thinfilm_h0",thinfilm_h0); 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; From 9dff1494157e7dade625d33cf5b5e4fa7a0d05de Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 6 Nov 2023 17:13:44 -0800 Subject: [PATCH 2/5] readability --- exec/thinFilm/main_driver.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/exec/thinFilm/main_driver.cpp b/exec/thinFilm/main_driver.cpp index 614cc4ce..f6b78432 100644 --- a/exec/thinFilm/main_driver.cpp +++ b/exec/thinFilm/main_driver.cpp @@ -135,8 +135,8 @@ void main_driver(const char* argv) // define DistributionMapping dmap.define(ba); - MultiFab height (ba, dmap, 1, 1); - MultiFab Laph (ba, dmap, 1, 1); + 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.); @@ -462,14 +462,14 @@ void main_driver(const char* argv) { 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) ); }); From f7bc077ac88102a790cc602d9606dd6d0a79c9e2 Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 6 Nov 2023 17:15:17 -0800 Subject: [PATCH 3/5] more readability --- exec/thinFilm/main_driver.cpp | 2 +- exec/thinFilm/thinfilm_functions.cpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/exec/thinFilm/main_driver.cpp b/exec/thinFilm/main_driver.cpp index f6b78432..d2b2e492 100644 --- a/exec/thinFilm/main_driver.cpp +++ b/exec/thinFilm/main_driver.cpp @@ -309,7 +309,7 @@ void main_driver(const char* argv) AMREX_D_TERM(const Array4 & gradhx = gradh[0].array(mfi);, const Array4 & gradhy = gradh[1].array(mfi);, const Array4 & gradhz = gradh[2].array(mfi);); - + const Array4 & h = height.array(mfi); amrex::ParallelFor(bx_x, bx_y, diff --git a/exec/thinFilm/thinfilm_functions.cpp b/exec/thinFilm/thinfilm_functions.cpp index 2c466f32..dde70c0b 100644 --- a/exec/thinFilm/thinfilm_functions.cpp +++ b/exec/thinFilm/thinfilm_functions.cpp @@ -28,6 +28,7 @@ void InitializeThinfilmNamespace() { pp.get("thinfilm_h0",thinfilm_h0); pp.get("thinfilm_gamma",thinfilm_gamma); + pp.query("thinfilm_pertamp",thinfilm_pertamp); pp.query("thinfilm_hamaker",thinfilm_hamaker); From 9cfa60320dc0d3c73de667dc63ac259cf1ae1a1a Mon Sep 17 00:00:00 2001 From: Andy Nonaka Date: Mon, 6 Nov 2023 17:16:25 -0800 Subject: [PATCH 4/5] final cleanup - some testing now required --- exec/thinFilm/main_driver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exec/thinFilm/main_driver.cpp b/exec/thinFilm/main_driver.cpp index d2b2e492..56cf91b5 100644 --- a/exec/thinFilm/main_driver.cpp +++ b/exec/thinFilm/main_driver.cpp @@ -136,7 +136,7 @@ void main_driver(const char* argv) dmap.define(ba); MultiFab height(ba, dmap, 1, 1); - MultiFab Laph(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.); From 8d4d230f2e83e2443b2276bc75c3a67c31cf41c4 Mon Sep 17 00:00:00 2001 From: jbb Date: Thu, 9 Nov 2023 13:56:08 -0800 Subject: [PATCH 5/5] fixed bug in disjoining pressure --- exec/thinFilm/main_driver.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/exec/thinFilm/main_driver.cpp b/exec/thinFilm/main_driver.cpp index 56cf91b5..997c1783 100644 --- a/exec/thinFilm/main_driver.cpp +++ b/exec/thinFilm/main_driver.cpp @@ -397,6 +397,7 @@ void main_driver(const char* argv) }); } Laph.FillBoundary(geom.periodicity()); + disjoining.FillBoundary(geom.periodicity()); // compute gradLaph and gradDisjoining for ( MFIter mfi(height,TilingIfNotGPU()); mfi.isValid(); ++mfi ) { @@ -517,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]); }); @@ -539,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; @@ -558,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