diff --git a/AA b/AA deleted file mode 100644 index c42a2c3..0000000 --- a/AA +++ /dev/null @@ -1,50 +0,0 @@ -@kernel function HyperViscTracerKernel!(FTr,@Const(Tr),@Const(Rho),@Const(D),@Const(DW),@Const(dXdxI), - @Const(JJ),@Const(M),@Const(Glob)) - - I, J, iz = @index(Local, NTuple) - _,_,Iz,IF = @index(Global, NTuple) - - ColumnTilesDim = @uniform @groupsize()[3] - N = @uniform @groupsize()[1] - Nz = @uniform @ndrange()[3] - NF = @uniform @ndrange()[4] - - ID = I + (J - 1) * N - @inbounds ind = Glob[ID,IF] - - TrCol = @localmem eltype(FTr) (N,N, ColumnTilesDim) - TrCxCol = @localmem eltype(FTr) (N,N, ColumnTilesDim) - TrCyCol = @localmem eltype(FTr) (N,N, ColumnTilesDim) - if Iz <= Nz && IF <= NF - @inbounds TrCol[I,J,iz] = Tr[Iz,ind] / Rho[Iz,ind] - end - @synchronize - - ID = I + (J - 1) * N - @inbounds ind = Glob[ID,IF] - - if Iz <= Nz && IF <= NF - @inbounds Dxc = D[I,1] * TrCol[1,J,iz] - @inbounds Dyc = D[J,1] * TrCol[I,1,iz] - for k = 2 : N - @inbounds Dxc += D[I,k] * TrCol[k,J,iz] - @inbounds Dyc += D[J,k] * TrCol[I,k,iz] - end - @views @inbounds (GradDx, GradDy) = Grad12(Dxc,Dyc,dXdxI[1:2,1:2,:,ID,Iz,IF],JJ[ID,:,Iz,IF]) - @views @inbounds (tempx, tempy) = Contra12(GradDx,GradDy,dXdxI[1:2,1:2,:,ID,Iz,IF]) - @inbounds TrCxCol[I,J,iz] = tempx - @inbounds TrCyCol[I,J,iz] = tempy - end - - @synchronize - - ID = I + (J - 1) * N - @inbounds ind = Glob[ID,IF] - if Iz <= Nz && IF <= NF - @inbounds DivTr = DW[I,1] * TrCxCol[1,J,iz] + DW[J,1] * TrCyCol[I,1,iz] - for k = 2 : N - @inbounds DivTr += DW[I,k] * TrCxCol[k,J,iz] + DW[J,k] * TrCyCol[I,k,iz] - end - @inbounds @atomic FTr[Iz,ind] += DivTr / M[Iz,ind] - end -end diff --git a/Examples/testCartGrid.jl b/Examples/testCartGrid.jl new file mode 100644 index 0000000..8aa1b01 --- /dev/null +++ b/Examples/testCartGrid.jl @@ -0,0 +1,155 @@ +import CGDycore: + Thermodynamics, Examples, Parallels, Models, Grids, Outputs, Integration, GPU, DyCore +using MPI +using Base +using CUDA +using AMDGPU +using Metal +using KernelAbstractions +using StaticArrays +using ArgParse +using MPI + + +# Model +parsed_args = DyCore.parse_commandline() +Problem = parsed_args["Problem"] +ProfRho = parsed_args["ProfRho"] +ProfTheta = parsed_args["ProfTheta"] +PertTh = parsed_args["PertTh"] +ProfVel = parsed_args["ProfVel"] +ProfpBGrd = parsed_args["ProfpBGrd"] +ProfRhoBGrd = parsed_args["ProfRhoBGrd"] +HorLimit = parsed_args["HorLimit"] +Upwind = parsed_args["Upwind"] +Damping = parsed_args["Damping"] +Relax = parsed_args["Relax"] +StrideDamp = parsed_args["StrideDamp"] +NumV = parsed_args["NumV"] +NumTr = parsed_args["NumTr"] +BoundaryWE = parsed_args["BoundaryWE"] +BoundarySN = parsed_args["BoundarySN"] +BoundaryBT = parsed_args["BoundaryBT"] +Thermo = parsed_args["Thermo"] +RefProfile = parsed_args["RefProfile"] +Profile = parsed_args["Profile"] +Curl = parsed_args["Curl"] +ModelType = parsed_args["ModelType"] +Equation = parsed_args["Equation"] +Microphysics = parsed_args["Microphysics"] +TypeMicrophysics = parsed_args["TypeMicrophysics"] +RelCloud = parsed_args["RelCloud"] +Rain = parsed_args["Rain"] +#Orography +TopoS = parsed_args["TopoS"] +P1 = parsed_args["P1"] +P2 = parsed_args["P2"] +P3 = parsed_args["P3"] +P4 = parsed_args["P4"] + +# Parallel +Decomp = parsed_args["Decomp"] +SimDays = parsed_args["SimDays"] +SimHours = parsed_args["SimHours"] +SimMinutes = parsed_args["SimMinutes"] +SimSeconds = parsed_args["SimSeconds"] +SimTime = parsed_args["SimTime"] +dtau = parsed_args["dtau"] +IntMethod = parsed_args["IntMethod"] +Table = parsed_args["Table"] +GridType = parsed_args["GridType"] +Coriolis = parsed_args["Coriolis"] +CoriolisType = parsed_args["CoriolisType"] +Source = parsed_args["Source"] +VerticalDiffusion = parsed_args["VerticalDiffusion"] +JacVerticalDiffusion = parsed_args["JacVerticalDiffusion"] +JacVerticalAdvection = parsed_args["JacVerticalAdvection"] +VerticalDiffusionMom = parsed_args["VerticalDiffusionMom"] +SurfaceFlux = parsed_args["SurfaceFlux"] +SurfaceFluxMom = parsed_args["SurfaceFluxMom"] +# Grid +nx = parsed_args["nx"] +ny = parsed_args["ny"] +nz = parsed_args["nz"] +H = parsed_args["H"] +Stretch = parsed_args["Stretch"] +StretchType = parsed_args["StretchType"] +OrdPoly = parsed_args["OrdPoly"] +Lx = parsed_args["Lx"] +Ly = parsed_args["Ly"] +x0 = parsed_args["x0"] +y0 = parsed_args["y0"] +# Viscosity +HyperVisc = parsed_args["HyperVisc"] +HyperDCurl = parsed_args["HyperDCurl"] +HyperDGrad = parsed_args["HyperDGrad"] +HyperDDiv = parsed_args["HyperDDiv"] +# Output +PrintDays = parsed_args["PrintDays"] +PrintHours = parsed_args["PrintHours"] +PrintMinutes = parsed_args["PrintMinutes"] +PrintSeconds = parsed_args["PrintSeconds"] +PrintTime = parsed_args["PrintTime"] +PrintStartTime = parsed_args["PrintStartTime"] +# Device +Device = parsed_args["Device"] +GPUType = parsed_args["GPUType"] +FloatTypeBackend = parsed_args["FloatTypeBackend"] +NumberThreadGPU = parsed_args["NumberThreadGPU"] + +MPI.Init() + +if Device == "CPU" || Device == "CPU_P" + backend = CPU() +elseif Device == "GPU" || Device == "GPU_P" + if GPUType == "CUDA" + backend = CUDABackend() + CUDA.allowscalar(false) + CUDA.device!(MPI.Comm_rank(MPI.COMM_WORLD)) + elseif GPUType == "AMD" + backend = ROCBackend() + AMDGPU.allowscalar(false) + elseif GPUType == "Metal" + backend = MetalBackend() + Metal.allowscalar(true) + end +else + backend = CPU() +end + +if FloatTypeBackend == "Float64" + FTB = Float64 +elseif FloatTypeBackend == "Float32" + FTB = Float32 +else + @show "False FloatTypeBackend" + stop +end +Param = Examples.Parameters(FTB,Problem) + +KernelAbstractions.synchronize(backend) + +OrdPolyZ=1 +nx = 10 +ny = 20 +nz = 10 +Lx = 1.0 +Ly = 1.0 +x0 = 0.0 +y0 = 0.0 + +Boundary = Grids.Boundary() +Boundary.WE = "" +Boundary.SN = "" + + +GridCart = Grids.CartGrid(backend,FTB,nx,ny,Lx,Ly,x0,y0,Grids.OrientFaceCart,Boundary,nz) +@show GridCart.NumFaces +vtkSkeletonMeshCart = Outputs.vtkStruct{Float64}(backend,GridCart) +c = ones(FTB,GridCart.NumFaces) +for i = 1 : GridCart.NumFaces + c[i] = i +end +Outputs.vtkSkeleton(vtkSkeletonMeshCart,"Cart", 1, 1,c) + + diff --git a/Examples/testInterpolateMesh.jl b/Examples/testInterpolateMesh.jl index 6253f1c..af4239c 100644 --- a/Examples/testInterpolateMesh.jl +++ b/Examples/testInterpolateMesh.jl @@ -1,5 +1,5 @@ import CGDycore: - Examples, Parallels, Models, Grids, Outputs, Integration, GPU, DyCore + Examples, Parallels, Models, Grids, Outputs, Integration, GPU, DyCore using MPI using Base using CUDA @@ -11,23 +11,7 @@ using ArgParse using MPI -function main() - FTB = Float64 - Problem = "AdvectionSphereSlottedCylinder" - Param = Examples.Parameters(FTB,Problem) - Phys=DyCore.PhysParameters{FTB}() - Profile = Examples.DivergentSphereExample()(Param,Phys) - - nz = 1 - nPanel = 30 - RefineLevel = 5 - RadEarth = 1.0 - FT = Float64 - backend = CPU() - - SrcGrid = Grids.CubedGrid(backend,FT,nPanel,Grids.OrientFaceSphere,RadEarth,nz) -# DestGrid = Grids.TriangularGrid(backend,FT,RefineLevel,RadEarth,nz) - DestGrid = Grids.DelaunayGrid(backend,FT,RefineLevel,RadEarth,nz) +function Interpolate(SrcGrid,DestGrid,Profile,Example) Inter = Grids.interpolate(SrcGrid,DestGrid) #Test @@ -47,35 +31,38 @@ function main() end c2 = Inter * c1 vtkSkeletonMeshSrc = Outputs.vtkStruct{Float64}(backend,SrcGrid) - Outputs.vtkSkeleton(vtkSkeletonMeshSrc,"Source1CuSp", 1, 1, c1) + Outputs.vtkSkeleton(vtkSkeletonMeshSrc,Example*"Source", 1, 1, c1) vtkSkeletonMeshDest = Outputs.vtkStruct{Float64}(backend,DestGrid) - Outputs.vtkSkeleton(vtkSkeletonMeshDest,"Destination1Del", 1, 1, c2) + Outputs.vtkSkeleton(vtkSkeletonMeshDest,Example*"Destin", 1, 1, c2) +end - nPanelSrc = 40 - nPanelDest = 47 - SrcGrid = Grids.CubedGrid(backend,FT,nPanelSrc,Grids.OrientFaceSphere,RadEarth,nz) - DestGrid = Grids.CubedGrid(backend,FT,nPanelDest,Grids.OrientFaceSphere,RadEarth,nz) +FTB = Float64 +backend = CPU() - Inter = Grids.interpolate(SrcGrid,DestGrid) - #Test - c1 = ones(SrcGrid.NumFaces) - c2 = Inter * c1 - @show maximum(c2) - @show minimum(c2) - for iFD = 1 : DestGrid.NumFaces - if c2[iFD] < 0.9 - @show iFD,c2[iFD] - end - end - for iFS = 1 : SrcGrid.NumFaces - Mid = SVector{3}(SrcGrid.Faces[iFS].Mid.x,SrcGrid.Faces[iFS].Mid.y,SrcGrid.Faces[iFS].Mid.z) - _,_,_,_,Tr = Profile(Mid,0.0) - c1[iFS] = Tr - end - c2 = Inter * c1 - vtkSkeletonMeshSrc = Outputs.vtkStruct{Float64}(backend,SrcGrid) - Outputs.vtkSkeleton(vtkSkeletonMeshSrc,"Source2CuSp", 1, 1, c1) - vtkSkeletonMeshDest = Outputs.vtkStruct{Float64}(backend,DestGrid) - Outputs.vtkSkeleton(vtkSkeletonMeshDest,"Destination2CuSp", 1, 1, c2) - return nothing -end +Problem = "AdvectionSphereSlottedCylinder" +Param = Examples.Parameters(FTB,Problem) +Phys=DyCore.PhysParameters{FTB}() +Profile = Examples.DivergentSphereExample()(Param,Phys) + +nz = 1 +nPanel = 30 +RefineLevel = 5 +RadEarth = 1.0 + +SrcGrid = Grids.CubedGrid(backend,FTB,nPanel,Grids.OrientFaceSphere,RadEarth,nz) +DestGrid = Grids.DelaunayGrid(backend,FTB,RefineLevel,RadEarth,nz) +Interpolate(SrcGrid,DestGrid,Profile,"CubeDelau") + +nPanelSrc = 40 +nPanelDest = 47 +SrcGrid = Grids.CubedGrid(backend,FTB,nPanelSrc,Grids.OrientFaceSphere,RadEarth,nz) +DestGrid = Grids.CubedGrid(backend,FTB,nPanelDest,Grids.OrientFaceSphere,RadEarth,nz) +Interpolate(SrcGrid,DestGrid,Profile,"CubeCube") + +nPanelSrc = 40 +nLon = 200 +nLat = 100 +LatB = (1.0 - 0.001) * pi / 2 +SrcGrid = Grids.CubedGrid(backend,FTB,nPanelSrc,Grids.OrientFaceSphere,RadEarth,nz) +DestGrid = Grids.SphericalGrid(backend,FTB,nLon,nLat,LatB,Grids.OrientFaceSphere,RadEarth,nz) +Interpolate(SrcGrid,DestGrid,Profile,"CubeSphere") diff --git a/Examples/testNHSphere.jl b/Examples/testNHSphere.jl index df7ec4a..569d2de 100644 --- a/Examples/testNHSphere.jl +++ b/Examples/testNHSphere.jl @@ -69,6 +69,7 @@ Stretch = parsed_args["Stretch"] StretchType = parsed_args["StretchType"] TopoS = parsed_args["TopoS"] GridType = parsed_args["GridType"] +AdaptGridType = parsed_args["AdaptGridType"] RadEarth = parsed_args["RadEarth"] ScaleFactor = parsed_args["ScaleFactor"] # CG Element @@ -95,6 +96,12 @@ FloatTypeBackend = parsed_args["FloatTypeBackend"] NumberThreadGPU = parsed_args["NumberThreadGPU"] MPI.Init() +comm = MPI.COMM_WORLD +Proc = MPI.Comm_rank(comm) + 1 +ProcNumber = MPI.Comm_size(comm) +ParallelCom = DyCore.ParallelComStruct() +ParallelCom.Proc = Proc +ParallelCom.ProcNumber = ProcNumber if Device == "CPU" backend = CPU() @@ -122,6 +129,7 @@ else @show "False FloatTypeBackend" stop end + Param = Examples.Parameters(FTB,Problem) KernelAbstractions.synchronize(backend) @@ -212,6 +220,9 @@ if RadEarth == 0.0 end end +Grid, Exchange = Grids.InitGrid(backend,FTB,OrdPoly,nz,nPanel,GridType,Decomp,RadEarth,Model,ParallelCom) + + Topography = (TopoS=TopoS,H=H,Rad=RadEarth) #Topography if TopoS == "BaroWaveHill" @@ -222,9 +233,13 @@ else TopoProfile = Examples.Flat()() end +Grid.AdaptGrid = Grids.AdaptGrid(FTB,AdaptGridType,H) + @show "InitSphere" -(CG, Metric, Exchange, Global) = DyCore.InitSphere(backend,FTB,OrdPoly,OrdPolyZ,nz,nPanel,H, - GridType,Topography,Decomp,Model,Phys,RadEarth,TopoProfile) +(CG, Metric, Global) = DyCore.InitSphere(backend,FTB,OrdPoly,OrdPolyZ,H,Topography,Model, + Phys,TopoProfile,Exchange,Grid,ParallelCom) +#(CG, Metric, Exchange, Global) = DyCore.InitSphere(backend,FTB,OrdPoly,OrdPolyZ,nz,nPanel,H, +# GridType,Topography,Decomp,Model,Phys,RadEarth,TopoProfile) # Initial values Examples.InitialProfile!(Model,Problem,Param,Phys) diff --git a/Examples/testSphericalGrid.jl b/Examples/testSphericalGrid.jl new file mode 100644 index 0000000..bd5802e --- /dev/null +++ b/Examples/testSphericalGrid.jl @@ -0,0 +1,149 @@ +import CGDycore: + Thermodynamics, Examples, Parallels, Models, Grids, Outputs, Integration, GPU, DyCore +using MPI +using Base +using CUDA +using AMDGPU +using Metal +using KernelAbstractions +using StaticArrays +using ArgParse +using MPI + + +# Model +parsed_args = DyCore.parse_commandline() +Problem = parsed_args["Problem"] +ProfRho = parsed_args["ProfRho"] +ProfTheta = parsed_args["ProfTheta"] +PertTh = parsed_args["PertTh"] +ProfVel = parsed_args["ProfVel"] +ProfpBGrd = parsed_args["ProfpBGrd"] +ProfRhoBGrd = parsed_args["ProfRhoBGrd"] +HorLimit = parsed_args["HorLimit"] +Upwind = parsed_args["Upwind"] +Damping = parsed_args["Damping"] +Relax = parsed_args["Relax"] +StrideDamp = parsed_args["StrideDamp"] +NumV = parsed_args["NumV"] +NumTr = parsed_args["NumTr"] +BoundaryWE = parsed_args["BoundaryWE"] +BoundarySN = parsed_args["BoundarySN"] +BoundaryBT = parsed_args["BoundaryBT"] +Thermo = parsed_args["Thermo"] +RefProfile = parsed_args["RefProfile"] +Profile = parsed_args["Profile"] +Curl = parsed_args["Curl"] +ModelType = parsed_args["ModelType"] +Equation = parsed_args["Equation"] +Microphysics = parsed_args["Microphysics"] +TypeMicrophysics = parsed_args["TypeMicrophysics"] +RelCloud = parsed_args["RelCloud"] +Rain = parsed_args["Rain"] +#Orography +TopoS = parsed_args["TopoS"] +P1 = parsed_args["P1"] +P2 = parsed_args["P2"] +P3 = parsed_args["P3"] +P4 = parsed_args["P4"] + +# Parallel +Decomp = parsed_args["Decomp"] +SimDays = parsed_args["SimDays"] +SimHours = parsed_args["SimHours"] +SimMinutes = parsed_args["SimMinutes"] +SimSeconds = parsed_args["SimSeconds"] +SimTime = parsed_args["SimTime"] +dtau = parsed_args["dtau"] +IntMethod = parsed_args["IntMethod"] +Table = parsed_args["Table"] +GridType = parsed_args["GridType"] +Coriolis = parsed_args["Coriolis"] +CoriolisType = parsed_args["CoriolisType"] +Source = parsed_args["Source"] +VerticalDiffusion = parsed_args["VerticalDiffusion"] +JacVerticalDiffusion = parsed_args["JacVerticalDiffusion"] +JacVerticalAdvection = parsed_args["JacVerticalAdvection"] +VerticalDiffusionMom = parsed_args["VerticalDiffusionMom"] +SurfaceFlux = parsed_args["SurfaceFlux"] +SurfaceFluxMom = parsed_args["SurfaceFluxMom"] +# Grid +nx = parsed_args["nx"] +ny = parsed_args["ny"] +nz = parsed_args["nz"] +H = parsed_args["H"] +Stretch = parsed_args["Stretch"] +StretchType = parsed_args["StretchType"] +OrdPoly = parsed_args["OrdPoly"] +Lx = parsed_args["Lx"] +Ly = parsed_args["Ly"] +x0 = parsed_args["x0"] +y0 = parsed_args["y0"] +# Viscosity +HyperVisc = parsed_args["HyperVisc"] +HyperDCurl = parsed_args["HyperDCurl"] +HyperDGrad = parsed_args["HyperDGrad"] +HyperDDiv = parsed_args["HyperDDiv"] +# Output +PrintDays = parsed_args["PrintDays"] +PrintHours = parsed_args["PrintHours"] +PrintMinutes = parsed_args["PrintMinutes"] +PrintSeconds = parsed_args["PrintSeconds"] +PrintTime = parsed_args["PrintTime"] +PrintStartTime = parsed_args["PrintStartTime"] +# Device +Device = parsed_args["Device"] +GPUType = parsed_args["GPUType"] +FloatTypeBackend = parsed_args["FloatTypeBackend"] +NumberThreadGPU = parsed_args["NumberThreadGPU"] + +MPI.Init() + +if Device == "CPU" || Device == "CPU_P" + backend = CPU() +elseif Device == "GPU" || Device == "GPU_P" + if GPUType == "CUDA" + backend = CUDABackend() + CUDA.allowscalar(false) + CUDA.device!(MPI.Comm_rank(MPI.COMM_WORLD)) + elseif GPUType == "AMD" + backend = ROCBackend() + AMDGPU.allowscalar(false) + elseif GPUType == "Metal" + backend = MetalBackend() + Metal.allowscalar(true) + end +else + backend = CPU() +end + +if FloatTypeBackend == "Float64" + FTB = Float64 +elseif FloatTypeBackend == "Float32" + FTB = Float32 +else + @show "False FloatTypeBackend" + stop +end +Param = Examples.Parameters(FTB,Problem) + +KernelAbstractions.synchronize(backend) + +OrdPolyZ=1 +nLon = 40 +nLat = 20 +nz = 10 +Rad = 1.0 +LatB = (1.0 - 0.001) * pi / 2 + +GridSpherical = Grids.SphericalGrid(backend,FTB,nLon,nLat,LatB,Grids.OrientFaceSphere,Rad,nz) + +@show GridSpherical.NumFaces +vtkSkeletonMeshSpherical = Outputs.vtkStruct{Float64}(backend,GridSpherical) +c = ones(FTB,GridSpherical.NumFaces) +for i = 1 : GridSpherical.NumFaces + c[i] = i +end +Outputs.vtkSkeleton(vtkSkeletonMeshSpherical,"Spherical", 1, 1,c) + + diff --git a/Examples/testSurfaceFlux.jl b/Examples/testSurfaceFlux.jl new file mode 100644 index 0000000..f7e25ee --- /dev/null +++ b/Examples/testSurfaceFlux.jl @@ -0,0 +1,161 @@ +import CGDycore: + Examples, Parallels, Models, Grids, Outputs, Integration, GPU, DyCore, Surfaces +using MPI +using Base +using CUDA +using AMDGPU +using Metal +using KernelAbstractions +using StaticArrays +using ArgParse + +# Model +parsed_args = DyCore.parse_commandline() +Problem = parsed_args["Problem"] +ProfRho = parsed_args["ProfRho"] +ProfTheta = parsed_args["ProfTheta"] +PertTh = parsed_args["PertTh"] +ProfVel = parsed_args["ProfVel"] +ProfVelGeo = parsed_args["ProfVelGeo"] +RhoVPos = parsed_args["RhoVPos"] +RhoCPos = parsed_args["RhoCPos"] +RhoIPos = parsed_args["RhoIPos"] +RhoRPos = parsed_args["RhoRPos"] +HorLimit = parsed_args["HorLimit"] +Upwind = parsed_args["Upwind"] +Damping = parsed_args["Damping"] +Relax = parsed_args["Relax"] +StrideDamp = parsed_args["StrideDamp"] +Geos = parsed_args["Geos"] +Coriolis = parsed_args["Coriolis"] +CoriolisType = parsed_args["CoriolisType"] +Buoyancy = parsed_args["Buoyancy"] +Equation = parsed_args["Equation"] +RefProfile = parsed_args["RefProfile"] +ProfpBGrd = parsed_args["ProfpBGrd"] +ProfRhoBGrd = parsed_args["ProfRhoBGrd"] +Microphysics = parsed_args["Microphysics"] +TypeMicrophysics = parsed_args["TypeMicrophysics"] +RelCloud = parsed_args["RelCloud"] +Rain = parsed_args["Rain"] +Source = parsed_args["Source"] +Forcing = parsed_args["Forcing"] +VerticalDiffusion = parsed_args["VerticalDiffusion"] +JacVerticalDiffusion = parsed_args["JacVerticalDiffusion"] +JacVerticalAdvection = parsed_args["JacVerticalAdvection"] +SurfaceFlux = parsed_args["SurfaceFlux"] +SurfaceFluxMom = parsed_args["SurfaceFluxMom"] +NumV = parsed_args["NumV"] +NumTr = parsed_args["NumTr"] +Curl = parsed_args["Curl"] +ModelType = parsed_args["ModelType"] +Thermo = parsed_args["Thermo"] +# Parallel +Decomp = parsed_args["Decomp"] +# Time integration +SimDays = parsed_args["SimDays"] +SimHours = parsed_args["SimHours"] +SimMinutes = parsed_args["SimMinutes"] +SimSeconds = parsed_args["SimSeconds"] +StartAverageDays = parsed_args["StartAverageDays"] +dtau = parsed_args["dtau"] +IntMethod = parsed_args["IntMethod"] +Table = parsed_args["Table"] +# Grid +nz = parsed_args["nz"] +nPanel = parsed_args["nPanel"] +H = parsed_args["H"] +Stretch = parsed_args["Stretch"] +StretchType = parsed_args["StretchType"] +TopoS = parsed_args["TopoS"] +GridType = parsed_args["GridType"] +RadEarth = parsed_args["RadEarth"] +# CG Element +OrdPoly = parsed_args["OrdPoly"] +# Viscosity +HyperVisc = parsed_args["HyperVisc"] +HyperDCurl = parsed_args["HyperDCurl"] +HyperDGrad = parsed_args["HyperDGrad"] +HyperDRhoDiv = parsed_args["HyperDRhoDiv"] +HyperDDiv = parsed_args["HyperDDiv"] +HyperDDivW = parsed_args["HyperDDivW"] +# Output +PrintDays = parsed_args["PrintDays"] +PrintHours = parsed_args["PrintHours"] +PrintMinutes = parsed_args["PrintMinutes"] +PrintSeconds = parsed_args["PrintSeconds"] +PrintStartTime = parsed_args["PrintStartTime"] +Flat = parsed_args["Flat"] + +# Device +Device = parsed_args["Device"] +GPUType = parsed_args["GPUType"] +FloatTypeBackend = parsed_args["FloatTypeBackend"] +NumberThreadGPU = parsed_args["NumberThreadGPU"] + +MPI.Init() + +Device = "CPU" +FloatTypeBackend = "Float64" + +if Device == "CPU" + backend = CPU() +elseif Device == "GPU" + if GPUType == "CUDA" + backend = CUDABackend() + CUDA.allowscalar(true) +# CUDA.device!(MPI.Comm_rank(MPI.COMM_WORLD)) + elseif GPUType == "AMD" + backend = ROCBackend() + AMDGPU.allowscalar(false) + elseif GPUType == "Metal" + backend = MetalBackend() + Metal.allowscalar(true) + end +else + backend = CPU() +end + +if FloatTypeBackend == "Float64" + FTB = Float64 +elseif FloatTypeBackend == "Float32" + FTB = Float32 +else + @show "False FloatTypeBackend" + stop +end + +# Physical parameters +Phys = DyCore.PhysParameters{FTB}() + +z0M = 0.05 +z0h = 0.05 +dz = 30.0 +RiB = 0.4 +LandClass = 1 +TSurf = 280.0 +T1 = 280.0 +VT = 2.0 +uf = Surfaces.Businger{FTB}() +theta = FTB(300) +thetaS = FTB(299.0) +#thetaS = FTB(295) +Surfaces.MOSTIteration(uf,z0M,z0h,dz,VT,theta,thetaS,LandClass,Phys) +stop + + +RiB = Phys.Grav/TSurf*((T1-TSurf)*(dz-z0M))/(VT*VT) +@show RiB +Cm, Ch = Surfaces.RichardIteration(z0M,z0h,dz,RiB,LandClass) +@show Cm, Ch + + +LandClass = 10 +TSurf = 265.0 +T1 = 262.0 +VT = 2.0 +RiB = Phys.Grav/TSurf*((T1-TSurf)*(dz-z0M))/(VT*VT) +Cm, Ch = Surfaces.RichardIteration(z0M,z0h,dz,RiB,LandClass) +@show Cm, Ch + + diff --git a/Jobs/NHSphere/JobNHBaroWaveHillDrySphere b/Jobs/NHSphere/JobNHBaroWaveHillDrySphere index 32c5190..208f815 100755 --- a/Jobs/NHSphere/JobNHBaroWaveHillDrySphere +++ b/Jobs/NHSphere/JobNHBaroWaveHillDrySphere @@ -19,24 +19,26 @@ mpirun -n 6 julia --project Examples/testNHSphere.jl \ --Buoyancy=true \ --Damping=true \ --StrideDamp=20000 \ - --Relax=1.0e-1 \ + --Relax=1.0e-3 \ --Decomp="EqualArea" \ --SimDays=10 \ --SimSeconds=0 \ - --PrintSeconds=600 \ + --PrintSeconds=0 \ --PrintMinutes=0 \ - --PrintHours=0 \ + --PrintHours=2 \ --PrintDays=0 \ + --PrintStartTime=0 \ --StartAverageDays=100 \ --Flat=true \ --dtau=150 \ --IntMethod="Rosenbrock" \ --Table="SSP-Knoth" \ --TopoS="BaroWaveHill" \ - --Stretch=true \ + --Stretch=false \ --StretchType="Exp" \ + --AdaptGridType="Sleve" \ --GridType="CubedSphere" \ - --nz=64 \ + --nz=31 \ --nPanel=30 \ --H=31000.0 \ --OrdPoly=3 \ diff --git a/Jobs/NHSphere/JobNHSchaerHillDrySphere b/Jobs/NHSphere/JobNHSchaerHillDrySphere index 484b712..e0b1505 100755 --- a/Jobs/NHSphere/JobNHSchaerHillDrySphere +++ b/Jobs/NHSphere/JobNHSchaerHillDrySphere @@ -39,6 +39,7 @@ mpirun -n 1 julia --project Examples/testNHSphere.jl \ --Stretch=false \ --StretchType="Exp" \ --GridType="CubedSphere" \ + --AdaptGridType="Sleve" \ --ScaleFactor=166.7 \ --nz=64 \ --nPanel=30 \ diff --git a/src/CGDycore.jl b/src/CGDycore.jl index 5677a7b..8efdd2a 100644 --- a/src/CGDycore.jl +++ b/src/CGDycore.jl @@ -29,6 +29,7 @@ include("Integration/Integration.jl") include("DyCore/DyCore.jl") include("GPU/GPU.jl") include("FiniteElements/FiniteElements.jl") +include("Surfaces/Surfaces.jl") OOP = 5 diff --git a/src/DG/QuadTringle.jl b/src/DG/QuadTringle.jl new file mode 100644 index 0000000..54623f4 --- /dev/null +++ b/src/DG/QuadTringle.jl @@ -0,0 +1,41 @@ +function QuadTriangle(Order) + if Order == 2 + w = zeros(3) + ksi = zeros(3,2) + w .= 0.333333333333333 + ksi[1,1] = 0.166666666666667 + ksi[1,2] = 0.666666666666667 + ksi[2,1] = 0.166666666666667 + ksi[2,2] = 0.166666666666667 + ksi[3,1] = 0.666666666666667 + ksi[3,2] = 0.166666666666667 + @. w *= 2 + @. ksi = 2 * ksi - 1 + elseif Order == 4 + w = zeros(6) + ksi = zeros(6,2) + w[1] = 0.223381589678011 + w[2] = 0.223381589678011 + w[3] = 0.223381589678011 + w[4] = 0.109951743655322 + w[5] = 0.109951743655322 + w[6] = 0.109951743655322 + + ksi[1,1] = 0.445948490915965 + ksi[2,1] = 0.445948490915965 + ksi[3,1] = 0.108103018168070 + ksi[4,1] = 0.091576213509771 + ksi[5,1] = 0.091576213509771 + ksi[6,1] = 0.816847572980459 + + ksi[1,2] = .108103018168070 + ksi[2,2] = 0.445948490915965 + ksi[3,2] = 0.445948490915965 + ksi[4,2] = 0.816847572980459 + ksi[5,2] = 0.091576213509771 + ksi[6,2] = 0.091576213509771 + @. w *= 2 + @. ksi = 2 * ksi - 1 + end + return w, ksi +end diff --git a/src/DyCore/DiscretizationCG.jl b/src/DyCore/DiscretizationCG.jl index 7c0ff65..e6862c0 100644 --- a/src/DyCore/DiscretizationCG.jl +++ b/src/DyCore/DiscretizationCG.jl @@ -27,7 +27,7 @@ function DiscretizationCG(backend,FT,Jacobi,CG::CGQuad,Exchange,Global,zs) end copyto!(FGPU,F) if Global.Grid.Form == "Sphere" - Grids.JacobiSphere3GPU!(Metric.X,Metric.dXdxI,Metric.J,CG,FGPU,Grid.z,zs,Grid.Rad) + Grids.JacobiSphere3GPU!(Global.Grid.AdaptGrid,Metric.X,Metric.dXdxI,Metric.J,CG,FGPU,Grid.z,zs,Grid.Rad) else Grids.JacobiDG3GPU!(Metric.X,Metric.dXdxI,Metric.J,CG,FGPU,Grid.z,zs) end diff --git a/src/DyCore/FiniteElementTri.jl b/src/DyCore/FiniteElementTri.jl index d9ea44a..c58bf9b 100644 --- a/src/DyCore/FiniteElementTri.jl +++ b/src/DyCore/FiniteElementTri.jl @@ -71,12 +71,12 @@ function RT0Struct{FT}(backend,::Grids.Tri,Grid) where FT<:AbstractFloat if Type == Grid.Type DoF = 3 Fun = Array{Polynomial,2}(undef,3,2) - Fun[1,1] = Polynomial([0,1]) - Fun[1,2] = Polynomial([-1,1]) - Fun[2,1] = Polynomial([0,0,0]) - Fun[2,2] = Polynomial([0,0,0]) - Fun[3,1] = Polynomial([0,0,0]) - Fun[3,2] = Polynomial([0,0,0]) + Fun[1,1] = Polynomial([sqrt(2)/2,sqrt(2)/2]) + Fun[1,2] = Polynomial([sqrt(2)/2,sqrt(2)/2]) + Fun[2,1] = Polynomial([1/2,-1/2]) + Fun[2,2] = Polynomial([1/2,+1/2]) + Fun[3,1] = Polynomial([1/2,+1/2]) + Fun[3,2] = Polynomial([1/2,-1/2]) Glob = KernelAbstractions.zeros(backend,Int,3,Grid.NumFaces) GlobCPU = zeros(Int,3,Grid.NumFaces) NumG = Grid.NumEdgesI + Grid.NumEdgesB diff --git a/src/DyCore/InitDriver.jl b/src/DyCore/InitDriver.jl index 713c0b6..e91c244 100644 --- a/src/DyCore/InitDriver.jl +++ b/src/DyCore/InitDriver.jl @@ -1,81 +1,42 @@ -function InitSphere(backend,FT,OrdPoly,OrdPolyZ,nz,nPanel,H,GridType,Topography,Decomp,Model,Phys,RadEarth,TopoProfile) - - comm = MPI.COMM_WORLD - Proc = MPI.Comm_rank(comm) + 1 - ProcNumber = MPI.Comm_size(comm) - ParallelCom = ParallelComStruct() - ParallelCom.Proc = Proc - ParallelCom.ProcNumber = ProcNumber +function InitSphere(backend,FT,OrdPoly,OrdPolyZ,H,Topography,Model,Phys,TopoProfile,Exchange,Grid,ParallelCom) + nz = Grid.nz + Proc = ParallelCom.Proc + ProcNumber = ParallelCom.ProcNumber TimeStepper = TimeStepperStruct{FT}(backend) -# Grid = Grids.GridStruct{FT}(backend,nz,Topography) - - if GridType == "HealPix" - # Grid=CGDycore.InputGridH("Grid/mesh_H12_no_pp.nc", - # CGDycore.OrientFaceSphere,Phys.RadEarth,Grid) - Grid=Grids.InputGridH(backend,FT,"Grid/mesh_H24_no_pp.nc", Grids.OrientFaceSphere,RadEarth,nz) - elseif GridType == "SQuadGen" - Grid = Grids.InputGrid("Grid/baroclinic_wave_2deg_x4.g",Grids.OrientFaceSphere,RadEarth,nz) - elseif GridType == "Msh" - Grid = Grids.InputGridMsh(backend,FT,"Grid/Quad.msh",Grids.OrientFaceSphere,RadEarth,nz) - elseif GridType == "CubedSphere" - Grid = Grids.CubedGrid(backend,FT,nPanel,Grids.OrientFaceSphere,RadEarth,nz) - elseif GridType == "TriangularSphere" - IcosahedronGrid = Grids.CreateIcosahedronGrid() - RefineLevel = 0 - for iRef = 1 : RefineLevel - Grids.RefineEdgeTriangularGrid!(IcosahedronGrid) - Grids.RefineFaceTriangularGrid!(IcosahedronGrid) - end - Grids.NumberingTriangularGrid!(IcosahedronGrid) - Grid = Grids.TriangularGridToGrid(IcosahedronGrid,RadEarth,Grid) - else - @show "False GridType" - end - - if Decomp == "Hilbert" - Parallels.HilbertFaceSphere!(Grid) - CellToProc = Grids.Decompose(Grid,ProcNumber) - elseif Decomp == "EqualArea" - CellToProc = Grids.DecomposeEqualArea(Grid,ProcNumber) - else - CellToProc = ones(Int,Grid.NumFaces) - println(" False Decomp method ") - end - SubGrid = Grids.ConstructSubGrid(Grid,CellToProc,Proc) + Output = OutputStruct() + DoF = (OrdPoly + 1) * (OrdPoly + 1) + Global = GlobalStruct{FT}(backend,Grid,Model,TimeStepper,ParallelCom,Output,DoF,nz, + Model.NumV,Model.NumTr,()) + CG = CGQuad{FT}(backend,OrdPoly,OrdPolyZ,Global.Grid) if Model.Stretch - if Model.StretchType == "ICON" + if Model.StretchType == "ICON" sigma = 1.0 lambda = 3.16 - Grids.AddStretchICONVerticalGrid!(SubGrid,nz,H,sigma,lambda) - elseif Model.StretchType == "Exp" - Grids.AddExpStretchVerticalGrid!(SubGrid,nz,H,30.0,1500.0) + Grids.AddStretchICONVerticalGrid!(Global.Grid,nz,H,sigma,lambda) + elseif Model.StretchType == "Exp" + Grids.AddExpStretchVerticalGrid!(Global.Grid,nz,H,30.0,1500.0) else - Grids.AddVerticalGrid!(SubGrid,nz,H) + Grids.AddVerticalGrid!(Global.Grid,nz,H) end - else - Grids.AddVerticalGrid!(SubGrid,nz,H) + else + Grids.AddVerticalGrid!(Global.Grid,nz,H) end - Exchange = Parallels.ExchangeStruct{FT}(backend,SubGrid,OrdPoly,CellToProc,Proc,ProcNumber,Model.HorLimit) - Output = OutputStruct() - DoF = (OrdPoly + 1) * (OrdPoly + 1) - Global = GlobalStruct{FT}(backend,SubGrid,Model,TimeStepper,ParallelCom,Output,DoF,nz, - Model.NumV,Model.NumTr,()) - CG = CGQuad{FT}(backend,OrdPoly,OrdPolyZ,Global.Grid) - if Topography.TopoS == "EarthOrography" zS = Grids.Orography(backend,FT,CG,Exchange,Global) else zS = Grids.Orography(backend,FT,CG,Exchange,Global,TopoProfile) end + (CG,Metric) = DiscretizationCG(backend,FT,Grids.JacobiSphere3,CG,Exchange,Global,zS) # Output Orography - Global.Output.dTol = 2*pi / nPanel + Global.Output.dTol = 2*pi / 30 +# Global.Output.dTol = 1.e-8 Output.Flat=true nzTemp = Global.Grid.nz Global.Grid.nz = 1 @@ -90,7 +51,7 @@ function InitSphere(backend,FT,OrdPoly,OrdPolyZ,nz,nPanel,H,GridType,Topography, Outputs.unstructured_vtkPartition(vtkCachePart,Global.Grid.NumFaces,Proc,ProcNumber) Global.Grid.nz = nzTemp - return CG,Metric,Exchange,Global + return CG, Metric, Global end diff --git a/src/DyCore/parse_commandline.jl b/src/DyCore/parse_commandline.jl index dd77ba3..28e604f 100644 --- a/src/DyCore/parse_commandline.jl +++ b/src/DyCore/parse_commandline.jl @@ -317,6 +317,11 @@ function parse_commandline() arg_type = String default = "" + "--AdaptGridType" + help = "Type of boundary following coordinate" + arg_type = String + default = "Sleve" + "--nx" help = "Number of horizontal grid cell in x-direction" arg_type = Int diff --git a/src/Examples/topography.jl b/src/Examples/topography.jl index 436fb26..2486b4a 100644 --- a/src/Examples/topography.jl +++ b/src/Examples/topography.jl @@ -27,9 +27,9 @@ function (profile::Flat)() end Base.@kwdef struct BaroWaveHill{T} <: Topography - h0::T = 0.25 * 2.e3 - lon1::T = 72 * pi / 180 - lon2::T = 140 * pi / 180 + h0::T = 2.e3 + lon1::T = (72 + 10) * pi / 180 + lon2::T = (140 + 10) * pi / 180 lonBar::T = 7 * pi / 180 c::T = 0.5 * lonBar*(-log(0.1))^(-1/2) lat1::T = pi / 4 diff --git a/src/GPU/FcnGPU.jl b/src/GPU/FcnGPU.jl index 0c1f09a..23179aa 100644 --- a/src/GPU/FcnGPU.jl +++ b/src/GPU/FcnGPU.jl @@ -104,8 +104,6 @@ end function FcnGPU!(F,U,FE,Metric,Phys,Cache,Exchange,Global,Param,DiscType) - - @show sum(abs.(U)) backend = get_backend(F) FT = eltype(F) Glob = FE.Glob diff --git a/src/GPU/surface.jl b/src/GPU/surface.jl index 622fdd3..cf94bf8 100644 --- a/src/GPU/surface.jl +++ b/src/GPU/surface.jl @@ -35,6 +35,7 @@ end end Base.@kwdef struct MOSurface <: SurfaceValues end + function (::MOSurface)(Phys,Param,uPos,vPos,wPos) function SurfaceData(U,p,dXdxI,nS) FT = eltype(U) diff --git a/src/Grids/CartGrid.jl b/src/Grids/CartGrid.jl index dc7f2a9..ec6fa7c 100644 --- a/src/Grids/CartGrid.jl +++ b/src/Grids/CartGrid.jl @@ -211,6 +211,7 @@ function CartGrid(backend,FT,nx::Int,ny::Int,lx::Float64,ly::Float64,x0::Float64 NumBoundaryFaces = 0 NumEdgesI = 0 NumEdgesB = 0 + AdaptGrid = "" return GridStruct{FT, typeof(z)}( nz, @@ -235,6 +236,7 @@ function CartGrid(backend,FT,nx::Int,ny::Int,lx::Float64,ly::Float64,x0::Float64 nBar, colors, NumBoundaryFaces, + AdaptGrid, ) end diff --git a/src/Grids/CubedGrid.jl b/src/Grids/CubedGrid.jl index c971ec2..57b2b22 100644 --- a/src/Grids/CubedGrid.jl +++ b/src/Grids/CubedGrid.jl @@ -339,6 +339,7 @@ function CubedGrid(backend,FT,n,OrientFace,Rad,nz) nBar3 = zeros(0,0) nBar = zeros(0,0) NumBoundaryFaces = 0 + AdaptGrid = "" return GridStruct{FT, typeof(z)}( nz, @@ -363,6 +364,7 @@ function CubedGrid(backend,FT,n,OrientFace,Rad,nz) nBar, colors, NumBoundaryFaces, + AdaptGrid, ) end diff --git a/src/Grids/GridStruct.jl b/src/Grids/GridStruct.jl index e0ac2bc..7191c04 100644 --- a/src/Grids/GridStruct.jl +++ b/src/Grids/GridStruct.jl @@ -38,6 +38,7 @@ mutable struct GridStruct{FT<:AbstractFloat, nBar::Array{FT, 2} colors::Array{Array{Int, 1}, 1} NumBoundaryFaces::Int + AdaptGrid::Any end function GridStruct{FT}(backend,nz) where FT <: AbstractFloat zP=zeros(nz) @@ -61,6 +62,7 @@ function GridStruct{FT}(backend,nz) where FT <: AbstractFloat nBar=zeros(0,0) colors=[[]] NumBoundaryFaces = 0 + AdaptGrid = "" return GridStruct{FT, typeof(z)}( nz, @@ -85,6 +87,7 @@ function GridStruct{FT}(backend,nz) where FT <: AbstractFloat nBar, colors, NumBoundaryFaces, + AdaptGrid, ) end diff --git a/src/Grids/Grids.jl b/src/Grids/Grids.jl index 5c6630d..970873e 100644 --- a/src/Grids/Grids.jl +++ b/src/Grids/Grids.jl @@ -58,5 +58,7 @@ include("polygon.jl") include("intersect.jl") include("interpolate.jl") include("TopographySmoothing.jl") +include("InitGrid.jl") +include("SphericalGrid.jl") end diff --git a/src/Grids/InitGrid.jl b/src/Grids/InitGrid.jl new file mode 100644 index 0000000..8cdf31f --- /dev/null +++ b/src/Grids/InitGrid.jl @@ -0,0 +1,38 @@ +function InitGrid(backend,FT,OrdPoly,nz,nPanel,GridType,Decomp,RadEarth,Model,ParallelCom) + + ProcNumber = ParallelCom.ProcNumber + Proc = ParallelCom.Proc + + if GridType == "HealPix" + # Grid=CGDycore.InputGridH("Grid/mesh_H12_no_pp.nc", + # CGDycore.OrientFaceSphere,Phys.RadEarth,Grid) + Grid=Grids.InputGridH(backend,FT,"Grid/mesh_H24_no_pp.nc", Grids.OrientFaceSphere,RadEarth,nz) + elseif GridType == "SQuadGen" + Grid = Grids.InputGrid("Grid/baroclinic_wave_2deg_x4.g",Grids.OrientFaceSphere,RadEarth,nz) + elseif GridType == "Msh" + Grid = Grids.InputGridMsh(backend,FT,"Grid/Quad.msh",Grids.OrientFaceSphere,RadEarth,nz) + elseif GridType == "CubedSphere" + Grid = Grids.CubedGrid(backend,FT,nPanel,Grids.OrientFaceSphere,RadEarth,nz) + elseif GridType == "TriangularSphere" + Grid = TriangularGrid(backend,FT,RefineLevel,RadEarth,nz) + elseif GridType == "DelaunaySphere" + Grid = DelaunayGrid(backend,FT,RefineLevel,RadEarth,nz) + else + @show "False GridType" + end + + if Decomp == "Hilbert" + Parallels.HilbertFaceSphere!(Grid) + CellToProc = Grids.Decompose(Grid,ProcNumber) + elseif Decomp == "EqualArea" + CellToProc = Grids.DecomposeEqualArea(Grid,ProcNumber) + else + CellToProc = ones(Int,Grid.NumFaces) + println(" False Decomp method ") + end + + SubGrid = Grids.ConstructSubGrid(Grid,CellToProc,Proc) + + Exchange = Parallels.ExchangeStruct{FT}(backend,SubGrid,OrdPoly,CellToProc,Proc,ProcNumber,Model.HorLimit) + return SubGrid, Exchange +end diff --git a/src/Grids/JacobiSphere3GPU.jl b/src/Grids/JacobiSphere3GPU.jl index 03321f5..f6a9b2d 100644 --- a/src/Grids/JacobiSphere3GPU.jl +++ b/src/Grids/JacobiSphere3GPU.jl @@ -21,7 +21,7 @@ function JacobiSphere2GPU!(X,dXdxI,J,FE,F,Rad) KJacobiSphere2Kernel!(X,dXdxI,J,FE.xw,FE.DS,F,Rad,ndrange=ndrange) end -function JacobiSphere3GPU!(X,dXdxI,J,FE,F,z,zs,Rad) +function JacobiSphere3GPU!(AdaptGrid,X,dXdxI,J,FE,F,z,zs,Rad) backend = get_backend(X) FT = eltype(X) @@ -36,7 +36,7 @@ function JacobiSphere3GPU!(X,dXdxI,J,FE,F,z,zs,Rad) KJacobiSphere3Kernel! = JacobiSphere3Kernel!(backend,group) - KJacobiSphere3Kernel!(X,dXdxI,J,FE.xw,FE.xwZ,FE.DS,F,z,Rad,zs,ndrange=ndrange) + KJacobiSphere3Kernel!(AdaptGrid,X,dXdxI,J,FE.xw,FE.xwZ,FE.DS,F,z,Rad,zs,ndrange=ndrange) KernelAbstractions.synchronize(backend) end @@ -61,7 +61,7 @@ end end end -@kernel function JacobiSphere3Kernel!(X,dXdxI,JJ,@Const(ksi),@Const(zeta),@Const(D), +@kernel function JacobiSphere3Kernel!(AdaptGrid,X,dXdxI,JJ,@Const(ksi),@Const(zeta),@Const(D), @Const(F),@Const(z),Rad,@Const(zs)) gi, gj, gk, gz, gF = @index(Group, NTuple) @@ -98,8 +98,8 @@ end F[3,3,IF] * (eltype(X)(1)+ksi1)*(eltype(X)(1)+ksi2) + F[4,3,IF] * (eltype(X)(1)-ksi1)*(eltype(X)(1)+ksi2)) zLoc = eltype(X)(1/2) * ((eltype(X)(1)-ksi3) * z1 + (eltype(X)(1)+ksi3) * z2) - hR[I,J,K,iz], D33 = GalChen(zLoc,H,zs[I,J,IF]) -# hR, D33 = Sleve(zLoc,H,zs) +# hR[I,J,K,iz], D33 = GalChen(zLoc,H,zs[I,J,IF]) + hR[I,J,K,iz], D33 = AdaptGrid(zLoc,zs[I,J,IF]) D33 = eltype(X)(1/2) * D33*(z2-z1) r = sqrt(X1 * X1 + X2 * X2 + X3 * X3) f = Rad / r @@ -268,22 +268,47 @@ end A[1,3] * (A[2,1] * A[3,2] - A[2,2] * A[3,1]) end -@inline function GalChen(zRef,H,zs) - z = zRef + (H - zRef) * zs / H - DzDzRef = eltype(zRef)(1) - zs / H - return z, DzDzRef + +abstract type AdaptGrid end + +Base.@kwdef struct GalChen <: AdaptGrid end + +function (::GalChen)(H) + function AdaptHeight(zRef,zs) + z = zRef + (H - zRef) * zs / H + DzDzRef = eltype(zRef)(1) - zs / H + return z, DzDzRef + end + return AdaptHeight end -@inline function Sleve(zRef,H,zs) - etaH = eltype(zRef)(.7) - s = eltype(zRef)(8/10) - eta = zRef / H - if eta <= etaH - z = eta * H + zs * sinh((etaH - eta) / s / etaH) / sinh(1 / s) - DzDzRef = eltype(zRef)(1) - zs / H / s / etaH * cosh((etaH - eta) / s / etaH) / sinh(1 / s) - else - z = eta * H - DzDzRef = eltype(zRef)(1) - end - return z, DzDzRef +Base.@kwdef struct Sleve{T} <: AdaptGrid + etaH::T = .7 + s::T = 8/10 + +end + +function (F::Sleve)(H) + (;s,etaH) = F + function AdaptHeight(zRef,zs) + eta = zRef / H + if eta <= etaH + z = eta * H + zs * sinh((etaH - eta) / s / etaH) / sinh(1 / s) + DzDzRef = eltype(zRef)(1) - zs / H / s / etaH * cosh((etaH - eta) / s / etaH) / sinh(1 / s) + else + z = eta * H + DzDzRef = eltype(zRef)(1) + end + return z, DzDzRef + end + return AdaptHeight +end + +function AdaptGrid(FT,Type,H) + + if Type == "GalChen" + AdaptGridFunction = Grids.GalChen()(H) + elseif Type == "Sleve" + AdaptGridFunction = Grids.Sleve{FT}()(H) + end end diff --git a/src/Grids/SphericalGrid.jl b/src/Grids/SphericalGrid.jl new file mode 100644 index 0000000..0072b10 --- /dev/null +++ b/src/Grids/SphericalGrid.jl @@ -0,0 +1,163 @@ +function SphericalGrid(backend,FT,nLon,nLat,LatB,OrientFace,Rad,nz) + + nBar=[ 0 1 0 1 + -1 0 -1 0]; + Dim=3; + Type=Quad() + Rad=Rad; + Form="Sphere"; + + NumNodes = nLon * (nLat + 1) + Nodes = map(1:NumNodes) do i + Node() + end + dLon = 2 * pi / nLon + dLat = 2 * LatB /nLat + P=zeros(Float64,3,nLon+1,nLat+1) + Lat = -LatB + @inbounds for iLat=1:nLat+1 + Lon = 0.0 + @inbounds for iLon=1:nLon+1 + P[:,iLon,iLat] .= sphere2cart(Lon,Lat,Rad) + Lon = Lon + dLon + end + Lat = Lat + dLat + end + NodeNumber = 1 + Lon = 0.0 + Lat = -LatB + for i = 1 : nLat + 1 + Lon = 0.0 + for j = 1 : nLon + PP = sphere2cart(Lon,Lat,Rad) + Nodes[NodeNumber] = Node(Point(PP),NodeNumber) + NodeNumber += 1 + Lon += dLon + end + Lat += dLat + end + + + NumEdges=nLon*nLat+nLon*(nLat+1) + NumEdgesLon=nLon*(nLat+1) + + Edges= map(1:NumEdges) do i + Edge([1,2],Nodes,0,0,"",0) + end + + EdgeNumber=1 + EdgeNumberLon=1 + EdgeNumberLat=1 + BC="" + N1=1 + N2=nLon+1 + @inbounds for iLat=1:nLat + @inbounds for iLon=1:nLon + Edges[EdgeNumber]=Edge([N1,N2],Nodes,EdgeNumber,EdgeNumber,"Lat",EdgeNumberLat) + EdgeNumber=EdgeNumber+1 + EdgeNumberLat=EdgeNumberLat+1 + N1=N1+1 + N2=N2+1 + end + end + N1=1 + N2=2 + TypeE = "Lon" + @inbounds for iLat=1:nLat+1 + @inbounds for iLon=1:nLon + if iLon == nLon + Edges[EdgeNumber] = Edge([N1,1+(iLat-1)*nLon],Nodes,EdgeNumber,EdgeNumber,TypeE,EdgeNumberLon) + EdgeNumber = EdgeNumber + 1 + EdgeNumberLon = EdgeNumberLon + 1 + N1 += 1 + N2 += 1 + else + Edges[EdgeNumber]=Edge([N1,N2],Nodes,EdgeNumber,EdgeNumber,TypeE,EdgeNumberLon) + EdgeNumber=EdgeNumber+1 + EdgeNumberLon=EdgeNumberLon+1 + N1=N1+1 + N2=N2+1 + if iLon==nLon + N1=N1+1 + N2=N2+1 + end + end + end + end + + NumFaces=nLon*nLat + Faces=map(1:NumFaces) do i + Face() + end + E1=nLon*nLat+1 + E3=E1+nLon + E2=2 + E4=1 + FaceNumber=1 + TypeF="o" + @inbounds for iLat=1:nLat + @inbounds for iLon=1:nLon + if iLon == nLon + (Faces[FaceNumber],Edges)=Face([E1,1+(iLat-1)*nLon,E3,E4],Nodes,Edges,FaceNumber,TypeF,OrientFace, + P=[P[:,iLon,iLat] P[:,iLon+1,iLat] P[:,iLon+1,iLat+1] P[:,iLon,iLat+1]]) + FaceNumber=FaceNumber+1 + E1=E1+1 + E3=E3+1 + else + (Faces[FaceNumber],Edges)=Face([E1,E2,E3,E4],Nodes,Edges,FaceNumber,TypeF,OrientFace, + P=[P[:,iLon,iLat] P[:,iLon+1,iLat] P[:,iLon+1,iLat+1] P[:,iLon,iLat+1]]) + FaceNumber=FaceNumber+1 + E1=E1+1 + E2=E2+1 + E3=E3+1 + E4=E4+1 + end + end + E2=E2+1 + E4=E4+1 + end + + Orientation!(Edges,Faces); + Renumbering!(Edges,Faces); + FacesInNodes!(Nodes,Faces) + + zP=zeros(nz) + z=KernelAbstractions.zeros(backend,FT,nz+1) + dzeta=zeros(nz) + H=0.0 + colors=[[]] + NumGhostFaces = 0 + nBar3 = zeros(0,0) + nBar = zeros(0,0) + NumBoundaryFaces = 0 + NumEdgesI = 0 + NumEdgesB = 0 + AdaptGrid = "" + return GridStruct{FT, + typeof(z)}( + nz, + zP, + z, + dzeta, + H, + NumFaces, + NumGhostFaces, + Faces, + NumEdges, + Edges, + NumNodes, + Nodes, + Form, + Type, + Dim, + Rad, + NumEdgesI, + NumEdgesB, + nBar3, + nBar, + colors, + NumBoundaryFaces, + AdaptGrid, + ) + +end diff --git a/src/Grids/SubGrid.jl b/src/Grids/SubGrid.jl index cca0f76..ba3fe6e 100644 --- a/src/Grids/SubGrid.jl +++ b/src/Grids/SubGrid.jl @@ -189,6 +189,7 @@ function ConstructSubGrid(GlobalGrid,Proc,ProcNumber) nBar3=zeros(0,0) nBar=zeros(0,0) Type = GlobalGrid.Type + AdaptGrid = GlobalGrid.AdaptGrid return GridStruct{FT, typeof(z)}( nz, @@ -213,6 +214,7 @@ function ConstructSubGrid(GlobalGrid,Proc,ProcNumber) nBar, colors, NumBoundaryFaces, + AdaptGrid, ) end diff --git a/src/Grids/Triangular.jl b/src/Grids/Triangular.jl index 5b13ad1..725fd72 100644 --- a/src/Grids/Triangular.jl +++ b/src/Grids/Triangular.jl @@ -453,6 +453,7 @@ function TriangularGridToGrid(backend,FT,TriangularGrid,Rad,nz) nBar3 = zeros(0,0) nBar = zeros(0,0) NumBoundaryFaces = 0 + AdaptGrid = "" return GridStruct{FT, typeof(z)}( nz, @@ -477,6 +478,7 @@ function TriangularGridToGrid(backend,FT,TriangularGrid,Rad,nz) nBar, colors, NumBoundaryFaces, + AdaptGrid, ) end @@ -551,6 +553,7 @@ function DelaunayGridToPolyGrid(backend,FT,TriangularGrid,Rad,nz) nBar3 = zeros(0,0) nBar = zeros(0,0) NumBoundaryFaces = 0 + AdaptGrid = "" return GridStruct{FT, typeof(z)}( @@ -576,6 +579,7 @@ function DelaunayGridToPolyGrid(backend,FT,TriangularGrid,Rad,nz) nBar, colors, NumBoundaryFaces, + AdaptGrid, ) end diff --git a/src/Grids/geometry_circle.jl b/src/Grids/geometry_circle.jl index 8ea8ceb..4d56719 100644 --- a/src/Grids/geometry_circle.jl +++ b/src/Grids/geometry_circle.jl @@ -7,7 +7,7 @@ function sphere2cart(lam,phi,r) x=cos(lam)*cos(phi)*r y=sin(lam)*cos(phi)*r z=sin(phi)*r -return [xy;z]; +return [x;y;z]; end function sphereDeg2cart(lam,phi,r) diff --git a/src/Surfaces/Surfaces.jl b/src/Surfaces/Surfaces.jl index 1d07feb..8082a40 100644 --- a/src/Surfaces/Surfaces.jl +++ b/src/Surfaces/Surfaces.jl @@ -1,4 +1,6 @@ module Surfaces - +include("solver.jl") +include("wrfscheme.jl") +include("functions.jl") end diff --git a/src/Surfaces/functions.jl b/src/Surfaces/functions.jl new file mode 100644 index 0000000..1af85db --- /dev/null +++ b/src/Surfaces/functions.jl @@ -0,0 +1,42 @@ +abstract type AbstractTransportType end +struct MomentumTransport <: AbstractTransportType end +struct HeatTransport <: AbstractTransportType end + +abstract type UniversalFunction end + +Base.@kwdef struct Businger{FT} <: UniversalFunction + am::FT = 4.7 + ah::FT = 4.7 + Pr::FT = 0.74 +end + +# Nishizawa2018 Eq. A7 +f_momentum(uf::Businger, zeta) = sqrt(sqrt(1 - 15 * zeta)) + +# Nishizawa2018 Eq. A8 +f_heat(uf::Businger, zeta) = sqrt(1 - 9 * zeta) + +function psi(uf::Businger, zeta, ::MomentumTransport) + FT = eltype(zeta) + if zeta < 0 + # Businger1971 Eq. A3 (zeta < 0) + f_m = f_momentum(uf, zeta) + log_term = log((1 + f_m)^2 * (1 + f_m^2) / 8) + return log_term - 2 * atan(f_m) + FT(pi) / 2 + else + # Businger1971 Eq. A3 (zeta >= 0) + return -uf.am * zeta + end +end + +function psi(uf::Businger, zeta, tt::HeatTransport) + if zeta < 0 + # Businger1971 Eq. A4 (zeta < 0) + f_h = f_heat(uf, zeta) + return 2 * log((1 + f_h) / 2) + else + # Businger1971 Eq. A4 (zeta >= 0) + return -uf.ah * zeta / uf.Pr + end +end + diff --git a/src/Surfaces/solver.jl b/src/Surfaces/solver.jl index 96b6e92..3ff422f 100644 --- a/src/Surfaces/solver.jl +++ b/src/Surfaces/solver.jl @@ -1,19 +1,56 @@ function Most() -#RiB = Grav/Tgn*((AbsT-Tgn)*(dz(iz)-zRauh))/(VT*VT+Eps) +#RiB = Grav/Tgn*((AbsT-Tgn)*(dz(iz)-z0M))/(VT*VT+Eps) end -function RichardIteration(zRauh,zRauhT,dz,RiB) +function VertProfile(uf,z0,z,L,tt) + log(z / z0) - psi(uf, z / L, tt) + + psi(uf, z0 / L, tt) +end +function MOSTIteration(uf,z0M,z0h,z,U,theta,thetaS,LandClass,Phys) + FT = eltype(theta) + Karm = 0.4 + Pr = 0.74 + Ri_b = min(Phys.Grav / thetaS *(theta - thetaS) * z / (U * U), 0.2) + L = (theta - thetaS) / U^2 + for Iter = 1 : 10 + f = Ri_b - Pr * z / L * + VertProfile(uf, z0h, z, L, HeatTransport()) / + VertProfile(uf, z0M, z, L, MomentumTransport())^2 + LP = L * (FT(1) + sqrt(eps(FT))) + flipsign(sqrt(eps(FT)), L) + fP = Ri_b - Pr * z / LP * + VertProfile(uf, z0h, z, LP, HeatTransport()) / + VertProfile(uf, z0M, z, LP, MomentumTransport())^2 + df = (fP - f) /(LP - L) + LNew = L - f/df + if abs(LNew-L) < eps(FT)^FT(1/3) + L = LNew + break + else + L = LNew + end + end + +end + +function RichardIteration(z0M,z0h,dz,RiB,LandClass) FT = eltype(RiB) - RiB = MAX(RiB,-FT(10)) + RiB = max(RiB,-FT(10)) chi = RiB - Iter=0 + Iter = 0 + Karm = 0.4 for Iter = 1 : 21 - fchi = FixPoint(chi,RiB) - chiS = chi*(FT(1)+sqrt(eps(FT)))+sign(sqrt(eps(FT)),chi) - fchiS = FixPoint(chiSi,RiB) + @show Iter, chi + fchi = RiB - RiBulk(dz, z0M, z0h, chi, RiB, LandClass) + @show fchi + chiS = chi * (FT(1) + sqrt(eps(FT))) + flipsign(sqrt(eps(FT)), chi) + fchiS = RiB - RiBulk(dz, z0M, z0h, chiS, RiB, LandClass) + @show fchiS dfdchi = (fchiS - fchi) /(chiS - chi) + @show fchi, dfdchi chiNew = chi - fchi/dfdchi + @show chi, chiNew if abs(chi-chiNew) < eps(FT)^FT(1/3) + chi = chiNew break elseif Iter>20 stop @@ -21,17 +58,15 @@ function RichardIteration(zRauh,zRauhT,dz,RiB) chi = chiNew end end - LandClass = 1 - TermH = log((dz+zRauh)/zRauhT) - - PsiH(chiNew*(FT(1)One+zRauh/dz),dz,zRauh,zRauhT,RiB,LandClass) + - PsiH(chiNew*zRauhT/dz,dz,zRauh,zRauhT,RiB,LandCLass) - TermM = log((dz+zRauh)/zRauh) - - PsiM(chiNew*(FT(1) + zRauh/dz),dz,zRauh,zRauhT,RiB,LandCLass) + - PsiM(chiNew*zRauh/dz,dz,zRauh,zRauhT,RiB,LandCLass) + @show Iter + TermH = log((dz + z0M) / z0h) - + PsiH(chi * (FT(1) + z0M / dz), dz, z0M, z0h, RiB, LandClass) + + PsiH(chi * z0h / dz, dz, z0M, z0h, RiB, LandClass) + TermM = log((dz+z0M)/z0M) - + PsiM(chi*(FT(1) + z0M / dz), dz, z0M, z0h, RiB, LandClass) + + PsiM(chi * z0M / dz, dz, z0M, z0h, RiB, LandClass) Cm = Karm * Karm / TermM / TermM - Ch = Karm Karm / TermM / TermH + Ch = Karm * Karm / TermM / TermH + return Cm, Ch end -@inline function FixPoint(chi,RiB) - RiB - RiBulk(chi,RiB) -end diff --git a/src/Surfaces/wrfscheme.jl b/src/Surfaces/wrfscheme.jl index 37eb743..3bafcfc 100644 --- a/src/Surfaces/wrfscheme.jl +++ b/src/Surfaces/wrfscheme.jl @@ -1,15 +1,18 @@ function RiBulk(z,z0M,z0h,Chi,RiB,LandClass) FT = eltype(z) - PsiHeat = log((z+z0)/z0h) - - PsiH(Chi*(FT(1)+z0/z),z,z0M,z0h,RiB,LandClass) + + PsiHeat = log((z+z0M)/z0h) - + PsiH(Chi*(FT(1)+z0M/z),z,z0M,z0h,RiB,LandClass) + PsiH(Chi*z0h/z,z,z0M,z0h,RiB,LandClass) - PsiMom = log((z+z0M)/z0) - + @show PsiH + PsiMom = log((z+z0M)/z0M) - PsiM(Chi*(FT(1)+z0M/z),z,z0M,z0h,RiB,LandClass) + PsiM(Chi*z0M/z,z,z0M,z0h,RiB,LandClass) + @show PsiM return Chi * PsiHeat / (PsiMom * PsiMom) end function PsiH(Chi,z,z0M,z0h,RiB,LandClass) + FT = eltype(z) if LandClass <= 8 a = FT(-5.3) b = FT(1.1) @@ -31,7 +34,7 @@ function PsiM(Chi,z,z0M,z0h,RiB,LandClass) a = FT(-6.1) b = FT(2.5) if RiB >= FT(0) - PsiM = a*log(Chi+(FT(1)+Chi**b)**(One/b)) + PsiM = a*log(Chi+(FT(1)+Chi^b)^(FT(1)/b)) else PsiM = (PsiMK(Chi)+Chi*Chi*PsiM_C(Chi))/(FT(1)+Chi*Chi) end @@ -60,17 +63,17 @@ end function PsiH_C(Chi) FT = eltype(Chi) - beta = 34.d0 - Var2 = (FT(1)-beta*Chi)**(FT(1)/FT(3)) - PsiH_C = FT(3)/FT(2)*log((Var2**FT(2)+Var2+FT(1))/FT(3)) - + beta = FT(34) + Var2 = (FT(1)-beta*Chi)^(FT(1)/FT(3)) + PsiH_C = FT(3)/FT(2)*log((Var2^FT(2)+Var2+FT(1))/FT(3)) - sqrt(FT(3))*atan((FT(2)*Var2+FT(1))/sqrt(FT(3))) + Pi/sqrt(FT(3)) end function PsiM_C(Chi) - beta = 10.d0 - Var2 = (FT(1)-beta*Chi)**(FT(1)/FT(3)) - PsiM_C = FT(3)/FT(2)*log((Var2**FT(2)+Var2+FT(1))/FT(3)) - + beta = FT(10) + Var2 = (FT(1)-beta*Chi)^(FT(1)/FT(3)) + PsiM_C = FT(3)/FT(2)*log((Var2^FT(2)+Var2+FT(1))/FT(3)) - sqrt(FT(3))*atan((FT(2)*Var2+FT(1))/sqrt(FT(3))) + pi/sqrt(FT(3)) end @@ -80,11 +83,12 @@ end function PsiH_SHEBA(Chi) FT = eltype(Chi) - ah = 5.d0 - bh = 5.d0 - ch = 3.d0 + ah = FT(5) + bh = FT(5) + ch = FT(3) grbh = sqrt(ch^FT(2)-FT(4)) + @show "PsiH_SHEBA",Chi PsiH_SHEBA = - bh/FT(2) * log(FT(1)+ch*Chi+Chi) + (-ah/grbh + bh*ch/(FT(2)*grbh)) * (log((FT(2)*Chi+ch-grbh)/(FT(2)*Chi+ch+grbh)) - @@ -95,8 +99,8 @@ end function PsiM_SHEBA(Chi) FT = eltype(Chi) - am = 5.d0 - bm = am/6.5d0 + am = FT(5) + bm = am/FT(6.5) grbm = ((FT(1)-bm)/bm)^(FT(1)/FT(3)) Var3 = (FT(1)+Chi)^(FT(1)/FT(3))