Hyperdiffusion for rho

CliMA · Jun 22, 2023 · 4ff5871 · 4ff5871
1 parent 486693e
commit 4ff5871
Show file tree

Hide file tree

Showing 4 changed files with 42 additions and 15 deletions.
diff --git a/Examples/testNHSphere.jl b/Examples/testNHSphere.jl
@@ -221,6 +221,6 @@ if ModelType == "VectorInvariant" || ModelType == "Advection"
 elseif ModelType == "Conservative"
   DiscType = Val(:Conservative)  
 end  
-nT = max(5,NumV + NumTr + 3)
+nT = max(9 + NumTr, NumV + NumTr)
 CGDycore.InitExchangeData3D(nz,nT,Global.Exchange)
 CGDycore.TimeStepper!(U,CGDycore.Fcn!,CGDycore.TransSphereX,CG,Global,Param,DiscType)
diff --git a/Jobs/JobNHHeldSuarezDrySphere b/Jobs/JobNHHeldSuarezDrySphere
@@ -1,4 +1,4 @@
-mpirun -n 1 julia --project Examples/testNHSphere.jl \
+mpirun -n 6 julia --project Examples/testNHSphere.jl \
   --Problem="HeldSuarezDrySphere" \
   --NumV=5 \
   --NumTr=0 \
@@ -14,9 +14,9 @@ mpirun -n 1 julia --project Examples/testNHSphere.jl \
   --VerticalDiffusionMom=false \
   --HorLimit=false \
   --Decomp="EqualArea" \
-  --SimDays=0 \
-  --SimSeconds=900 \
-  --PrintDays=20 \
+  --SimDays=100 \
+  --SimSeconds=0 \
+  --PrintDays=1 \
   --StartAverageDays=200 \
   --dtau=300 \
   --IntMethod="Rosenbrock" \

diff --git a/src/DyCore/Fcn.jl b/src/DyCore/Fcn.jl
@@ -196,18 +196,20 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
   @views Rot2 = Global.Cache.Temp1[:,:,2]
   @views Grad1 = Global.Cache.Temp1[:,:,3]
   @views Grad2 = Global.Cache.Temp1[:,:,4]
-  @views Div = Global.Cache.Temp1[:,:,NumV]
-  @views DivTr = Global.Cache.Temp1[:,:,NumV+1:NumV+NumTr]
-  @views Divw = Global.Cache.Temp1[:,:,NumV+NumTr+1]
-  @views JRho = Global.Cache.Temp1[:,:,NumV+NumTr+2]
-  @views JRhoF = Global.Cache.Temp1[:,:,NumV+NumTr+3]
+  @views Div = Global.Cache.Temp1[:,:,5]
+  @views DivTh = Global.Cache.Temp1[:,:,6]
+  @views Divw = Global.Cache.Temp1[:,:,7]
+  @views JRho = Global.Cache.Temp1[:,:,8]
+  @views JRhoF = Global.Cache.Temp1[:,:,9]
+  @views DivTr = Global.Cache.Temp1[:,:,9+1:9+NumTr]
   FCG=Global.Cache.FCC
   FwCG=Global.Cache.FwCC
   Rot1CG=Global.Cache.Rot1C
   Rot2CG=Global.Cache.Rot2C
   Grad1CG=Global.Cache.Grad1C
   Grad2CG=Global.Cache.Grad2C
   DivCG=Global.Cache.DivC
+  DivThCG=Global.Cache.DivThC
   DivwCG=Global.Cache.DivwC
   DivTrCG=Global.Cache.DivC
   @views RhoCG = Global.Cache.RhoCG[:,:,:]
@@ -235,6 +237,7 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
   @. Grad1 = 0.0
   @. Grad2 = 0.0
   @. Div = 0.0
+  @. DivTh = 0.0
   @. Divw = 0.0
   @. DivTr = 0.0
   @. F = 0.0
@@ -262,7 +265,9 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
     @views GradDiv!(Grad1CG,Grad2CG,v1CG,v2CG,CG,
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
-    @views DivRhoGrad!(DivCG,ThCG,RhoCG,CG,
+    @views DivRhoGrad!(DivCG,RhoCG,RhoCG,CG,
+     Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
+    @views DivRhoGrad!(DivThCG,ThCG,RhoCG,CG,
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
     @views DivGradF!(DivwCG,wCG,CG,
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
@@ -276,6 +281,7 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
           Grad1[iz,ind] += Grad1CG[iP,jP,iz] 
           Grad2[iz,ind] += Grad2CG[iP,jP,iz] 
           Div[iz,ind] += DivCG[iP,jP,iz] 
+          DivTh[iz,ind] += DivThCG[iP,jP,iz] 
           JRho[iz,ind] += (J[iP,jP,1,iz,iF] + J[iP,jP,2,iz,iF]) * RhoCG[iP,jP,iz]
         end
         @inbounds for iz=1:nz-1
@@ -327,7 +333,9 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
     @views GradDiv!(Grad1CG,Grad2CG,v1CG,v2CG,CG,
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
-    @views DivRhoGrad!(DivCG,ThCG,RhoCG,CG,
+    @views DivRhoGrad!(DivCG,RhoCG,RhoCG,CG,
+     Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
+    @views DivRhoGrad!(DivThCG,ThCG,RhoCG,CG,
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
     @views DivGradF!(DivwCG,wCG,CG,
      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache)
@@ -341,6 +349,7 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
           Grad1[iz,ind] += Grad1CG[iP,jP,iz] 
           Grad2[iz,ind] += Grad2CG[iP,jP,iz] 
           Div[iz,ind] += DivCG[iP,jP,iz] 
+          DivTh[iz,ind] += DivThCG[iP,jP,iz] 
           JRho[iz,ind] += (J[iP,jP,1,iz,iF] + J[iP,jP,2,iz,iF]) * RhoCG[iP,jP,iz]
         end
         @inbounds for iz=1:nz-1
@@ -388,6 +397,7 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
           Grad1CG[iP,jP,iz] = Grad1[iz,ind] / CG.M[iz,ind]
           Grad2CG[iP,jP,iz] = Grad2[iz,ind] / CG.M[iz,ind]
           DivCG[iP,jP,iz] = Div[iz,ind] / CG.M[iz,ind]
+          DivThCG[iP,jP,iz] = DivTh[iz,ind] / CG.M[iz,ind]
           zPG[iP,jP,iz] = zP[iz,ind]
           pBGrdCG[iP,jP,iz] = Global.pBGrd[iz,ind]
           RhoBGrdCG[iP,jP,iz] = Global.RhoBGrd[iz,ind]
@@ -412,7 +422,10 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
     @views GradDiv!(FCG[:,:,:,uPos],FCG[:,:,:,vPos],Grad1CG,Grad2CG,RhoCG,CG,
       Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache,
       Global.Model.HyperDGrad)
-    @views DivRhoGrad!(FCG[:,:,:,ThPos],DivCG,RhoCG,CG,
+    @views DivRhoGrad!(FCG[:,:,:,RhoPos],DivCG,RhoCG,CG,
+      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache,
+      Global.Model.HyperDDiv)
+    @views DivRhoGrad!(FCG[:,:,:,ThPos],DivThCG,RhoCG,CG,
       Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache,
       Global.Model.HyperDDiv)
     @views DivGradF!(FwCG,DivwCG,RhoCG,CG,
@@ -582,6 +595,7 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
           Grad1CG[iP,jP,iz] = Grad1[iz,ind] / CG.M[iz,ind]
           Grad2CG[iP,jP,iz] = Grad2[iz,ind] / CG.M[iz,ind]
           DivCG[iP,jP,iz] = Div[iz,ind] / CG.M[iz,ind]
+          DivThCG[iP,jP,iz] = DivTh[iz,ind] / CG.M[iz,ind]
           pBGrdCG[iP,jP,iz] = Global.pBGrd[iz,ind]
           RhoBGrdCG[iP,jP,iz] = Global.RhoBGrd[iz,ind]
           PresCG[iP,jP,iz] = PresG[iz,ind]
@@ -606,7 +620,10 @@ function Fcn!(F,U,CG,Global,Param,DiscType::Val{:VectorInvariant})
     @views GradDiv!(FCG[:,:,:,uPos],FCG[:,:,:,vPos],Grad1CG,Grad2CG,RhoCG,CG,
       Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache,
       Global.Model.HyperDGrad)
-    @views DivRhoGrad!(FCG[:,:,:,ThPos],DivCG,RhoCG,CG,
+    @views DivRhoGrad!(FCG[:,:,:,RhoPos],DivCG,RhoCG,CG,
+      Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache,
+      Global.Model.HyperDDiv)
+    @views DivRhoGrad!(FCG[:,:,:,ThPos],DivThCG,RhoCG,CG,
       Global.Metric.dXdxI[:,:,:,:,:,:,iF],Global.Metric.J[:,:,:,:,iF],Global.ThreadCache,
       Global.Model.HyperDDiv)
     @views DivGradF!(FwCG,DivwCG,RhoCG,CG,

diff --git a/src/Model/GlobalVariables.jl b/src/Model/GlobalVariables.jl
@@ -45,6 +45,7 @@ Rot2CG::Array{FT, 3}
 Grad1CG::Array{FT, 3}
 Grad2CG::Array{FT, 3}
 DivCG::Array{FT, 3}
+DivThCG::Array{FT, 3}
 DivwCG::Array{FT, 3}
 zPG::Array{FT, 3}
 pBGrdCG::Array{FT, 3}
@@ -54,6 +55,7 @@ Rot2C::Array{FT, 3}
 Grad1C::Array{FT, 3}
 Grad2C::Array{FT, 3}
 DivC::Array{FT, 3}
+DivThC::Array{FT, 3}
 DivwC::Array{FT, 3}
 KVCG::Array{FT, 3}
 Temp1::Array{FT, 3}
@@ -119,6 +121,7 @@ Rot2CG=zeros(FT,0,0,0)
 Grad1CG=zeros(FT,0,0,0)
 Grad2CG=zeros(FT,0,0,0)
 DivCG=zeros(FT,0,0,0)
+DivThCG=zeros(FT,0,0,0)
 DivwCG=zeros(FT,0,0,0)
 zPG=zeros(FT,0,0,0)
 pBGrdCG=zeros(FT,0,0,0)
@@ -128,6 +131,7 @@ Rot2C=zeros(FT,0,0,0)
 Grad1C=zeros(FT,0,0,0)
 Grad2C=zeros(FT,0,0,0)
 DivC=zeros(FT,0,0,0)
+DivThC=zeros(FT,0,0,0)
 DivwC=zeros(FT,0,0,0)
 KVCG=zeros(FT,0,0,0)
 Temp1=zeros(FT,0,0,0)
@@ -192,6 +196,7 @@ return CacheStruct(
   Grad1CG,
   Grad2CG,
   DivCG,
+  DivThCG,
   DivwCG,
   zPG,
   pBGrdCG,
@@ -201,6 +206,7 @@ return CacheStruct(
   Grad1C,
   Grad2C,
   DivC,
+  DivThC,
   DivwC,
   KVCG,
   Temp1,
@@ -268,6 +274,7 @@ Rot2CG=zeros(OP,OP,nz)
 Grad1CG=zeros(OP,OP,nz)
 Grad2CG=zeros(OP,OP,nz)
 DivCG=zeros(OP,OP,nz)
+DivThCG=zeros(OP,OP,nz)
 DivwCG=zeros(OP,OP,nz+1)
 zPG=zeros(OP,OP,nz)
 pBGrdCG=zeros(OP,OP,nz)
@@ -277,9 +284,10 @@ Rot2C=zeros(OP,OP,nz)
 Grad1C=zeros(OP,OP,nz)
 Grad2C=zeros(OP,OP,nz)
 DivC=zeros(OP,OP,nz)
+DivThC=zeros(OP,OP,nz)
 DivwC=zeros(OP,OP,nz+1)
 KVCG=zeros(OP,OP,nz)
-Temp1=zeros(nz,NumG,NumV+NumTr+3)
+Temp1=zeros(nz,NumG,max(NumV+NumTr,9+NumTr))
 k=zeros(0,0,0,0)
 Ymyn=zeros(0,0,0,0)
 Y=zeros(0,0,0,0)
@@ -341,6 +349,7 @@ return CacheStruct(
   Grad1CG,
   Grad2CG,
   DivCG,
+  DivThCG,
   DivwCG,
   zPG,
   pBGrdCG,
@@ -350,6 +359,7 @@ return CacheStruct(
   Grad1C,
   Grad2C,
   DivC,
+  DivThC,
   DivwC,
   KVCG,
   Temp1,