From 7499f89bb6e28a7dbca02faac55c72877a07b598 Mon Sep 17 00:00:00 2001 From: OsKnoth <50015520+OsKnoth@users.noreply.github.com> Date: Fri, 22 Nov 2024 11:33:25 +0100 Subject: [PATCH] Srfaces const Position --- Examples/testFEMSeiQuadKiteGal.jl | 2 +- Examples/testFEMSeiQuadKiteLinSh.jl | 27 +++++++++++++++++---------- src/Examples/initial.jl | 7 +------ src/FEMSei/TimestepperFEM.jl | 13 ++++++------- src/Surfaces/Surfaces.jl | 18 +++++++++--------- src/Surfaces/surface.jl | 6 +++--- 6 files changed, 37 insertions(+), 36 deletions(-) diff --git a/Examples/testFEMSeiQuadKiteGal.jl b/Examples/testFEMSeiQuadKiteGal.jl index f561800..c602819 100644 --- a/Examples/testFEMSeiQuadKiteGal.jl +++ b/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" diff --git a/Examples/testFEMSeiQuadKiteLinSh.jl b/Examples/testFEMSeiQuadKiteLinSh.jl index 520947c..4120635 100644 --- a/Examples/testFEMSeiQuadKiteLinSh.jl +++ b/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) diff --git a/src/Examples/initial.jl b/src/Examples/initial.jl index c886dd2..2c116dc 100755 --- a/src/Examples/initial.jl +++ b/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 diff --git a/src/FEMSei/TimestepperFEM.jl b/src/FEMSei/TimestepperFEM.jl index a43177d..9f41ea0 100644 --- a/src/FEMSei/TimestepperFEM.jl +++ b/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 diff --git a/src/Surfaces/Surfaces.jl b/src/Surfaces/Surfaces.jl index 91db658..00aded2 100644 --- a/src/Surfaces/Surfaces.jl +++ b/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, diff --git a/src/Surfaces/surface.jl b/src/Surfaces/surface.jl index 84c4614..4c7fbe5 100644 --- a/src/Surfaces/surface.jl +++ b/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