From 33988e957a5d5b1d91b7182956231d1f864adf0b Mon Sep 17 00:00:00 2001 From: Hannah Klion Date: Tue, 5 Sep 2023 12:46:34 -0700 Subject: [PATCH] bugfixes and debugging --- Source/TimeIntegration/ROMSX_Advance.cpp | 4 +++- Source/TimeIntegration/ROMSX_prsgrd.cpp | 8 ++++---- Source/TimeIntegration/ROMSX_rho_eos.cpp | 21 +++++++++++++++++++-- Source/TimeIntegration/ROMSX_rhs_uv_3d.cpp | 9 +++++---- 4 files changed, 31 insertions(+), 11 deletions(-) diff --git a/Source/TimeIntegration/ROMSX_Advance.cpp b/Source/TimeIntegration/ROMSX_Advance.cpp index dae8fdf9..c9371df1 100644 --- a/Source/TimeIntegration/ROMSX_Advance.cpp +++ b/Source/TimeIntegration/ROMSX_Advance.cpp @@ -161,6 +161,7 @@ ROMSX::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle { Array4 const& DC = mf_DC.array(mfi); Array4 const& Akv = (vec_Akv[lev])->array(mfi); + Array4 const& h = (vec_hOfTheConfusingName[lev])->array(mfi); Array4 const& Hz = (vec_Hz[lev])->array(mfi); Array4 const& Huon = (vec_Huon[lev])->array(mfi); Array4 const& Hvom = (vec_Hvom[lev])->array(mfi); @@ -325,7 +326,7 @@ ROMSX::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle set_massflux_3d(gbx2,1,0,uold,Huon,Hz,on_u,nnew); set_massflux_3d(gbx2,0,1,vold,Hvom,Hz,om_v,nnew); - rho_eos(gbx2,temp,salt,rho,rhoA,rhoS,pden,Hz,z_w,nrhs,N); + rho_eos(gbx2,temp,salt,rho,rhoA,rhoS,pden,Hz,z_w,h,nrhs,N); } if(solverChoice.use_prestep) { @@ -552,6 +553,7 @@ ROMSX::Advance (int lev, Real time, Real dt_lev, int /*iteration*/, int /*ncycle if(solverChoice.use_uv3dmix) { //u=u+(contributions from S-surfaces viscosity not scaled by dt)*dt*dx*dy //rufrc=rufrc + (contributions from S-surfaces viscosity not scaled by dt*dx*dy) + //For debugging version: uv3dmix(bx, u, v, rufrc, rvfrc, visc2_p, visc2_r, Hz, om_r, on_r, om_p, on_p, pm, pn, nrhs, nnew, dt_lev,false); uv3dmix(bx, u, v, rufrc, rvfrc, visc2_p, visc2_r, Hz, om_r, on_r, om_p, on_p, pm, pn, nrhs, nnew, dt_lev); if (verbose > 2) { amrex::PrintToFile("ru_afteruvmix").SetPrecision(18)<=0;k--) { Real cff= phi_bxD.contains(i+1,j,0) ? 2.0*aux(i,j,k)*aux(i+1,j,k) : 2.0*aux(i,j,k)*aux(i,j,k); if (cff>eps) { - Real cff1= 1.0_rt/(aux(i+1,j,k)+aux(i,j,k)); + Real cff1= phi_bxD.contains(i+1,j,0) ? 1.0_rt/(aux(i+1,j,k)+aux(i,j,k)) : 1.0_rt/(2.0*aux(i,j,k)); dZx(i,j,k)=cff*cff1; } else { dZx(i,j,k)=0.0; } Real cff1= phi_bxD.contains(i+1,j,0) ? 2.0*FC(i,j,k)*FC(i+1,j,k) : 2.0*FC(i,j,k)*FC(i,j,k); if (cff1>eps) { - Real cff2= 1.0_rt/(FC(i,j,k)+FC(i+1,j,k)); + Real cff2= phi_bxD.contains(i+1,j,0) ? 1.0_rt/(FC(i,j,k)+FC(i+1,j,k)) : 1.0_rt/(2.0*FC(i,j,k)); dRx(i,j,k)=cff1*cff2; } else { dRx(i,j,k)=0.0; @@ -188,14 +188,14 @@ ROMSX::prsgrd (const Box& phi_bx, for(int k=N;k>=0;k--) { Real cff= phi_bxD.contains(i,j+1,0) ? 2.0*aux(i,j,k)*aux(i,j+1,k) : 2.0*aux(i,j,k)*aux(i,j,k); if (cff>eps) { - Real cff1= 1.0_rt/(aux(i,j+1,k)+aux(i,j,k)); + Real cff1= phi_bxD.contains(i,j+1,0) ? 1.0_rt/(aux(i,j+1,k)+aux(i,j,k)) : 1.0_rt/(2.0*aux(i,j,k)); dZx(i,j,k)=cff*cff1; } else { dZx(i,j,k)=0.0; } Real cff1= phi_bxD.contains(i,j+1,0) ? 2.0*FC(i,j,k)*FC(i,j+1,k) : 2.0*FC(i,j,k)*FC(i,j,k); if (cff1>eps) { - Real cff2= 1.0_rt/(FC(i,j,k)+FC(i,j+1,k)); + Real cff2= phi_bxD.contains(i,j+1,0) ? 1.0_rt/(FC(i,j,k)+FC(i,j+1,k)) : 1.0_rt/(2.0*FC(i,j,k)); dRx(i,j,k)=cff1*cff2; } else { dRx(i,j,k)=0.0; diff --git a/Source/TimeIntegration/ROMSX_rho_eos.cpp b/Source/TimeIntegration/ROMSX_rho_eos.cpp index 5e0ab656..778c3f2b 100644 --- a/Source/TimeIntegration/ROMSX_rho_eos.cpp +++ b/Source/TimeIntegration/ROMSX_rho_eos.cpp @@ -12,6 +12,7 @@ ROMSX::rho_eos (const Box& phi_bx, Array4 pden, Array4 Hz, Array4 z_w, + Array4 h, const int nrhs, const int N) { @@ -40,12 +41,16 @@ ROMSX::rho_eos (const Box& phi_bx, amrex::ParallelFor(phi_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + //printf("%d %d %d %15.15g temp\n", i,j,k, temp(i,j,k,nrhs)); + //printf("%d %d %d %15.15g salt\n", i,j,k, salt(i,j,k,nrhs)); rho(i,j,k)=R0- R0*Tcoef*(temp(i,j,k,nrhs)-T0); rho(i,j,k)=rho(i,j,k)+ R0*Scoef*(salt(i,j,k,nrhs)-S0); rho(i,j,k)=rho(i,j,k)-1000.0_rt; pden(i,j,k)=rho(i,j,k); + //printf("%d %d %d %15.15g rhorhoeos\n", i,j,k, rho(i,j,k)); + //printf("%d %d %d %15.15g Hzrhoeos\n", i,j,k, Hz(i,j,k)); }); // //----------------------------------------------------------------------- @@ -56,9 +61,18 @@ ROMSX::rho_eos (const Box& phi_bx, amrex::ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real cff1=rho(i,j,N-1)*Hz(i,j,N); + //printf("%d %d %15.15g Hzstart\n", i,j, Hz(i,j,N)); + Real cff1=rho(i,j,N)*Hz(i,j,N); rhoS(i,j,0)=0.5_rt*cff1*Hz(i,j,N); rhoA(i,j,0)=cff1; + //printf("%d %d %d %15.15g cff1\n", i,j,N, cff1); + //printf("%d %d %d %15.15g rhoN\n", i,j,N, rho(i,j,N)); + //printf("%d %d %d %15.15g rhoA\n", i,j,N, rhoA(i,j,0)); + //printf("%d %d %d %15.15g rhoS\n", i,j,N, rhoS(i,j,0)); + //printf("%d %d %d %15.15g Hzend \n", i,j,N, Hz(i,j,N)); + //printf("%d %d %15.15g Hzend\n", i,j, Hz(i,j,N)); + printf("%d %d %15.15g %15.15g %15.15g %15.15g %15.15g cff1 rhoN rhoA rhoS Hz rhoeos\n", + i,j, cff1, rho(i,j,N), rhoA(i,j,0), rhoS(i,j,0), Hz(i,j,N)); }); amrex::ParallelFor(phi_bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -67,15 +81,18 @@ ROMSX::rho_eos (const Box& phi_bx, Real cff1=rho(i,j,N-k)*Hz(i,j,N-k); rhoS(i,j,0)=rhoS(i,j,0)+Hz(i,j,N-k)*(rhoA(i,j,0)+0.5_rt*cff1); rhoA(i,j,0)=rhoA(i,j,0)+cff1; + printf("%d %d %d %15.15g %15.15g %15.15g %15.15g cff1 rhoA rhoS Hz rhoeos\n", i,j,k, cff1, rhoA(i,j,0), rhoS(i,j,0), Hz(i,j,N-k)); } }); amrex::ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real cff2=1.0_rt/rho0; - Real cff1=1.0_rt/(z_w(i,j,N)-z_w(i,j,0)); + Real cff1=1.0_rt/(z_w(i,j,N)+h(i,j,0,0)); rhoA(i,j,0)=cff2*cff1*rhoA(i,j,0); rhoS(i,j,0)=2.0_rt*cff1*cff1*cff2*rhoS(i,j,0); + printf("%d %d %15.15g %15.15g %15.15g %15.15g z_wN z_w0\n", i,j, z_w(i,j,N), z_w(i,j,0), h(i,j,0,0),h(i,j,0,1)); + printf("%d %d %15.15g %15.15g %15.15g %15.15g cff2 cff1 rhoA rhoS rhoeos\n", i,j,cff2, cff1, rhoA(i,j,0), rhoS(i,j,0)); }); } diff --git a/Source/TimeIntegration/ROMSX_rhs_uv_3d.cpp b/Source/TimeIntegration/ROMSX_rhs_uv_3d.cpp index 47096a07..83bb3a58 100644 --- a/Source/TimeIntegration/ROMSX_rhs_uv_3d.cpp +++ b/Source/TimeIntegration/ROMSX_rhs_uv_3d.cpp @@ -234,10 +234,10 @@ ROMSX::rhs_3d (const Box& bx, const Box& gbx, ( cff1*( W(i ,j,N-1)+ W(i-1,j,N-1)) -cff2*( W(i+1,j,N-1)+ W(i-2,j,N-1)) ); - FC(i,j,0)=( cff1*(uold(i ,j,1,nrhs)+ uold(i,j,2,nrhs)) - -cff2*(uold(i ,j,1,nrhs)+ uold(i,j,3,nrhs)) )* - ( cff1*( W(i ,j,1)+ W(i-1,j,1)) - -cff2*( W(i+1,j,1)+ W(i-2,j,1)) ); + FC(i,j,0)=( cff1*(uold(i ,j,0,nrhs)+ uold(i,j,1,nrhs)) + -cff2*(uold(i ,j,0,nrhs)+ uold(i,j,2,nrhs)) )* + ( cff1*( W(i ,j,0)+ W(i-1,j,0)) + -cff2*( W(i+1,j,0)+ W(i-2,j,0)) ); // FC(i,0,-1)=0.0; } @@ -254,6 +254,7 @@ ROMSX::rhs_3d (const Box& bx, const Box& gbx, } ru(i,j,k,nrhs) -= cff; + //printf("%d %d %d %d %15.15g %15.15g %15.15g ru cff\n", i,j,k,nrhs,ru(i,j,k,nrhs), cff, FC(i,j,k)); }); Gpu::synchronize(); //This uses W being an extra grow cell sized amrex::ParallelFor(gbx1,