Srfaces const Position

Examples/testFEMSeiQuadKiteGal.jl

@@ -166,7 +166,7 @@ Examples.InitialProfile!(Model,Problem,Param,Phys)
 
 #Tri
 GridType = "CubedSphere"
-GridType = "DelaunaySphere"
+#GridType = "DelaunaySphere"
 ns=50
 ChangeOrient=2
 if GridType == "TriangularSphere"

Examples/testFEMSeiQuadKiteLinSh.jl

@@ -141,17 +141,24 @@ Phys = DyCore.PhysParameters{FTB}()
 Model = DyCore.ModelStruct{FTB}()
 
 RefineLevel = 5
-RadEarth = 1.0
 nz = 1
-nPanel = 40
+nPanel = 80
+nLat = 0
+nLon = 0
+LatB = 0.0
+ns = 40
 nQuadM = 4 #2
 nQuadS = 4 #3
+nQuad = 4 #3
 Decomp = "EqualArea"
 Problem = "LinearBlob"
 RadEarth = 1.0
-dtau = 0.00025
-nAdveVel = 100
-@show nAdveVel
+uMax = 1.0
+dtau = 2*pi*RadEarth/4/nPanel/uMax*0.1
+EndTime = 2 * pi
+nAdveVel = ceil(EndTime / dtau)
+nAdveVel = 1000
+@show dtau,nAdveVel
 Param = Examples.Parameters(FTB,Problem)
 Examples.InitialProfile!(Model,Problem,Param,Phys)
 
@@ -160,15 +167,15 @@ GridType = "CubedSphere"
 #GridType = "TriangularSphere"
 ChangeOrient=2
 if GridType == "TriangularSphere"
-  GridC, Exchange = Grids.InitGridSphere(backend,FTB,OrdPoly,nz,nPanel,RefineLevel,GridType,Decomp,
-    RadEarth,Model,ParallelCom;order=false,ChangeOrient=2)
+  GridC, Exchange = Grids.InitGridSphere(backend,FTB,OrdPoly,nz,nPanel,ns,nLat,nLon,LatB,
+    RefineLevel,GridType,Decomp,RadEarth,Model,ParallelCom;order=false,ChangeOrient=2)
 else  
-  GridC, Exchange = Grids.InitGridSphere(backend,FTB,OrdPoly,nz,nPanel,RefineLevel,GridType,Decomp,
-    RadEarth,Model,ParallelCom;order=false)
+  GridC, Exchange = Grids.InitGridSphere(backend,FTB,OrdPoly,nz,nPanel,ns,nLat,nLon,LatB,
+    RefineLevel,GridType,Decomp,RadEarth,Model,ParallelCom;order=false)
 end  
 
 Grid = Grids.Grid2KiteGrid(backend,FTB,GridC,Grids.OrientFaceSphere)
-vtkSkeletonMesh = Outputs.vtkStruct{Float64}(backend,Grid)
+vtkSkeletonMesh = Outputs.vtkStruct{Float64}(backend,Grid,Grid.NumFaces,Flat)
 
 CG1KiteP = FEMSei.CG1KitePrimalStruct{FTB}(Grids.Quad(),backend,Grid)
 CG1KiteDHDiv = FEMSei.CG1KiteDualHDiv{FTB}(Grids.Quad(),backend,Grid)

src/Examples/initial.jl

@@ -631,12 +631,7 @@ function (::HeldSuarezMoistExample)(Param,Phys)
     DeltaT =  kT * (T - Teq)
     F[5]  += - U[5] / T * DeltaT 
   end
-  @inline function TSurf(x)
-    FT = eltype(x)
-    (Lon,Lat,R)= Grids.cart2sphere(x[1],x[2],x[3])
-    TS = Param.DeltaTS * exp(-FT(0.5) * Lat^2 / Param.DeltaLat^2) + Param.TSMin
-  end
-  return profile, Force, TSurf
+  return profile, Force
 end
 
 Base.@kwdef struct SchaerSphereExample <: Example end

src/FEMSei/TimestepperFEM.jl

@@ -15,17 +15,18 @@ function TimeStepper(backend,FTB,U,dtau,Fcn,Model,Grid,nQuadM,nQuadS,Jacobi,nAdv
   VelSp = zeros(Grid.NumFaces,2)
   pC = zeros(Grid.NumFaces)
   uCurl = zeros(Model.DG.NumG)
-  @views Curl!(backend,FTB,uCurl,Model.DG,Uu,Model.RT,Model.ND,Grid,Jacobi!,nQuadM,
-    Model.Curl,UCache[uPosS:uPosE])
-  ConvertScalar!(backend,FTB,pC,uCurl,Model.DG,Grid,Jacobi!)
+  @show typeof(Model.DG)
+  ConvertScalar!(backend,FTB,pC,Up,Model.DG,Grid,Jacobi!)
   ConvertVelocityCart!(backend,FTB,VelCa,Uu,Model.RT,Grid,Jacobi!)
   ConvertVelocitySp!(backend,FTB,VelSp,Uu,Model.RT,Grid,Jacobi!)
   Outputs.vtkSkeleton!(vtkSkeletonMesh, GridType, Proc, ProcNumber, [pC VelCa VelSp], FileNumber)
+  stop
 
   time = 0.0
   UNew = similar(U)
   F = similar(U)
   time = 0.0
+  nPrint = ceil(nAdveVel/10)
   for i = 1 : nAdveVel
     @show i,time  
     Fcn(backend,FTB,F,U,Model,Grid,nQuadM,nQuadS,Jacobi!;UCache)
@@ -34,10 +35,8 @@ function TimeStepper(backend,FTB,U,dtau,Fcn,Model,Grid,nQuadM,nQuadS,Jacobi,nAdv
     @. UNew = U + 1/2 * dtau * F
     Fcn(backend,FTB,F,UNew,Model,Grid,nQuadM,nQuadS,Jacobi!;UCache)
     @. U = U + dtau * F
-    if mod(i,1) == 0
-      @views Curl!(backend,FTB,uCurl,Model.DG,Uu,Model.RT,Model.ND,Grid,Jacobi!,nQuadM,
-        Model.Curl,UCache[uPosS:uPosE])
-      ConvertScalar!(backend,FTB,pC,uCurl,Model.DG,Grid,Jacobi)
+    if mod(i,nPrint) == 0
+      ConvertScalar!(backend,FTB,pC,Up,Model.DG,Grid,Jacobi!)
       ConvertVelocityCart!(backend,FTB,VelCa,Uu,Model.RT,Grid,Jacobi)
       ConvertVelocitySp!(backend,FTB,VelSp,Uu,Model.RT,Grid,Jacobi)
       FileNumber += 1

src/Surfaces/Surfaces.jl

@@ -12,15 +12,15 @@ include("functions.jl")
 include("surface.jl")
 include("BoundaryLayer.jl")
 
-TSurfPos = 1
-RhoVSurfPos = 2
-uStarPos = 3
-CMPos = 4
-CTPos = 5
-CHPos = 6
-RiBSurfPos = 7
-hBLPos = 8
-LenSurfaceData = 8
+global const TSurfPos = 1
+global const RhoVSurfPos = 2
+global const uStarPos = 3
+global const CMPos = 4
+global const CTPos = 5
+global const CHPos = 6
+global const RiBSurfPos = 7
+global const hBLPos = 8
+global const LenSurfaceData = 8
 
 
 mutable struct LandUseData{FT<:AbstractFloat,

src/Surfaces/surface.jl

@@ -163,8 +163,8 @@ function SurfaceFlux(Phys,Param,ThPos,RhoPos,RhoVPos)
     Rm = Phys.Rd * RhoD + Phys.Rv * RhoV
     Cpml = Phys.Cpd * RhoD + Phys.Cpv * RhoV
     T = p / Rm
-    LatFlux = -SD[CTPos] * SD[uStarPos] * (RhoV - SD[RhoVSurfPos])
-    SensFlux = -SD[CHPos] * SD[uStarPos] * (T - SD[TSurfPos])
+    LatFlux = -SD[CHPos] * SD[uStarPos] * (RhoV - SD[RhoVSurfPos])
+    SensFlux = -SD[CTPos] * SD[uStarPos] * (T - SD[TSurfPos])
     FRho = LatFlux
     FRhoV = LatFlux
     PrePi=(p / Phys.p0)^(Rm / Cpml)
@@ -185,7 +185,7 @@ function SurfaceFlux(Phys,Param,ThPos,RhoPos)
     Rm = Phys.Rd * RhoD 
     Cpml = Phys.Cpd * RhoD 
     T = p / Rm
-    SensFlux = -CH * SD[uStarPos] * (T - TSurf)
+    SensFlux = -SD[CTPos] * SD[uStarPos] * (T - SD[TSurfPos])
     PrePi=(p / Phys.p0)^(Rm / Cpml)
     FRhoTh = RhoTh * SensFlux / T 
     FU[ThPos] += FRhoTh  / dz