diff --git a/BatchScripts/Mac/runHillAgnesi.sh b/BatchScripts/Mac/HillAgnesi.sh similarity index 100% rename from BatchScripts/Mac/runHillAgnesi.sh rename to BatchScripts/Mac/HillAgnesi.sh diff --git a/Examples/testFEMSeiConsShallow.jl b/Examples/ConsNonLinShallow.jl similarity index 83% rename from Examples/testFEMSeiConsShallow.jl rename to Examples/ConsNonLinShallow.jl index 02ffda6..620e1b0 100644 --- a/Examples/testFEMSeiConsShallow.jl +++ b/Examples/ConsNonLinShallow.jl @@ -152,10 +152,15 @@ LatB = 0.0 #Quad GridType = "CubedSphere" -nPanel = 40 +nPanel = 80 #GridType = "HealPix" ns = 57 +#Grid construction +RadEarth = Phys.RadEarth +Grid, Exchange = Grids.InitGridSphere(backend,FTB,OrdPoly,nz,nPanel,RefineLevel,ns, + nLat,nLon,LatB,GridType,Decomp,RadEarth,Model,ParallelCom) + print("Which Problem do you want so solve? \n") print("1 - GalewskiSphere\n\ 2 - HaurwitzSphere\n") @@ -165,38 +170,45 @@ a = parse(Int,text) if a == 1 Problem = "GalewskiSphere" Param = Examples.Parameters(FTB,Problem) - RadEarth = Phys.RadEarth - GridLength = Grids.GridLength(Grid) + GridLengthMin,GridLengthMax = Grids.GridLength(Grid) cS = sqrt(Phys.Grav * Param.H0G) - dtau = GridLength / cS * 0.1 + dtau = GridLengthMin / cS / sqrt(2) * .5 EndTime = 24 * 3600 # One day + PrintTime = 3600 nAdveVel = round(EndTime / dtau) dtau = EndTime / nAdveVel - nprint = ceil(nAdveVel/50) + nprint = ceil(PrintTime/dtau) GridTypeOut = GridType*"NonLinShallowGal" + @show GridLengthMin,GridLengthMax @show nAdveVel @show dtau @show nprint elseif a == 2 Problem = "HaurwitzSphere" Param = Examples.Parameters(FTB,Problem) - RadEarth = Phys.RadEarth - dtau = 30 #g=9.81, H=8000 - nAdveVel = ceil(Int,6*24*3600/dtau) + GridLengthMin,GridLengthMax = Grids.GridLength(Grid) + cS = sqrt(Phys.Grav * Param.h0) + dtau = GridLengthMin / cS / sqrt(2) * .3 + EndTime = 24 * 3600 # One day + PrintTime = 3600 + nAdveVel = round(EndTime / dtau) + dtau = EndTime / nAdveVel + nprint = ceil(PrintTime/dtau) GridTypeOut = GridType*"NonLinShallowHaurwitz" + @show GridLengthMin,GridLengthMax @show nAdveVel + @show dtau + @show nprint else print("Error") end println("The chosen Problem is ") Examples.InitialProfile!(Model,Problem,Param,Phys) -#Grid construction -Grid, Exchange = Grids.InitGridSphere(backend,FTB,OrdPoly,nz,nPanel,RefineLevel,ns,nLat,nLon,LatB,GridType,Decomp,RadEarth, - Model,ParallelCom) +#Output vtkSkeletonMesh = Outputs.vtkStruct{Float64}(backend,Grid,Grid.NumFaces,Flat) -#finite elements +#Finite elements DG = FEMSei.DG0Struct{FTB}(Grids.Quad(),backend,Grid) VecDG = FEMSei.VecDG0Struct{FTB}(Grids.Quad(),backend,Grid) RT = FEMSei.RT0Struct{FTB}(Grids.Quad(),backend,Grid) @@ -235,9 +247,8 @@ FEMSei.Project!(backend,FTB,Uh,DG,Grid,nQuad,FEMSei.Jacobi!,Model.InitialProfile FileNumber=0 VelSp = zeros(Grid.NumFaces,2) FEMSei.ConvertScalarVelocitySp!(backend,FTB,VelSp,Uhu,RT,Uh,DG,Grid,FEMSei.Jacobi!) -Outputs.vtkSkeleton!(vtkSkeletonMesh, "ConsShallow", Proc, ProcNumber, [Uh VelSp] ,FileNumber) +Outputs.vtkSkeleton!(vtkSkeletonMesh, GridTypeOut, Proc, ProcNumber, [Uh VelSp] ,FileNumber) -nprint = 10 for i = 1 : nAdveVel @show i,(i-1)*dtau/3600 @. F = 0 @@ -251,6 +262,18 @@ for i = 1 : nAdveVel FEMSei.CrossRhs!(backend,FTB,Fhu,RT,Uhu,RT,Grid,Grids.Quad(),nQuad,FEMSei.Jacobi!) FEMSei.GradHeightSquared!(backend,FTB,Fhu,RT,Uh,DG,Grid,Grids.Quad(),nQuad,FEMSei.Jacobi!) ldiv!(RT.LUM,Fhu) + @. UNew = U + 1 / 3 * dtau * F + @. F = 0 + # Tendency h + FEMSei.DivRhs!(backend,FTB,Fh,DG,UNewhu,RT,Grid,Grids.Quad(),nQuad,FEMSei.Jacobi!) + ldiv!(DG.LUM,Fh) + # Tendency hu + FEMSei.ProjectScalarHDivVecDG1!(backend,FTB,uRec,VecDG,UNewh,DG,UNewhu,RT,Grid, + Grids.Quad(),nQuad,FEMSei.Jacobi!) + FEMSei.DivMomentumVector!(backend,FTB,Fhu,RT,UNewhu,RT,uRec,VecDG,Grid,Grids.Quad(),nQuad,FEMSei.Jacobi!) + FEMSei.CrossRhs!(backend,FTB,Fhu,RT,UNewhu,RT,Grid,Grids.Quad(),nQuad,FEMSei.Jacobi!) + FEMSei.GradHeightSquared!(backend,FTB,Fhu,RT,UNewh,DG,Grid,Grids.Quad(),nQuad,FEMSei.Jacobi!) + ldiv!(RT.LUM,Fhu) @. UNew = U + 0.5 * dtau * F @. F = 0 # Tendency h @@ -268,7 +291,7 @@ for i = 1 : nAdveVel if mod(i,nprint) == 0 global FileNumber += 1 FEMSei.ConvertScalarVelocitySp!(backend,FTB,VelSp,Uhu,RT,Uh,DG,Grid,FEMSei.Jacobi!) - Outputs.vtkSkeleton!(vtkSkeletonMesh, "ConsShallow", Proc, ProcNumber, [Uh VelSp] ,FileNumber) + Outputs.vtkSkeleton!(vtkSkeletonMesh, GridTypeOut, Proc, ProcNumber, [Uh VelSp] ,FileNumber) end end @show "finished" diff --git a/Examples/testFEMSeiProjectionVectorTri.jl b/Examples/testFEMSeiProjectionVectorTri.jl new file mode 100644 index 0000000..846cca5 --- /dev/null +++ b/Examples/testFEMSeiProjectionVectorTri.jl @@ -0,0 +1,280 @@ +import CGDycore: + Examples, Parallels, Models, Grids, Outputs, Integration, GPU, DyCore, FEMSei, FiniteVolumes +using MPI +using Base +using CUDA +using AMDGPU +using Metal +using KernelAbstractions +using StaticArrays +using ArgParse +using LinearAlgebra + +# 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() +Flat = false #testen +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 + +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 + +# Physical parameters +Phys = DyCore.PhysParameters{FTB}() + +#ModelParameters +Model = DyCore.ModelStruct{FTB}() + +RefineLevel = 6 +nz = 1 +nQuad = 3 +nQuadM = 3 #2 +nQuadS = 3 #3 +Decomp = "EqualArea" +nLat = 0 +nLon = 0 +LatB = 0.0 + +#Quad +GridType = "TriangularSphere" +nPanel = 80 +#GridType = "HealPix" +ns = 57 + +print("Which Problem do you want so solve? \n") +print("1 - GalewskiSphere\n\ + 2 - HaurwitzSphere\n\ + 3 - LinearBlob\n\ + 4 - AdvectionSpherical\n") +text = readline() +a = parse(Int,text) +if a == 1 + Problem = "GalewskiSphere" + Param = Examples.Parameters(FTB,Problem) + RadEarth = Phys.RadEarth + dtau = 30 + nAdveVel = 100 #ceil(Int,6*24*3600/dtau) + GridTypeOut = GridType*"NonLinShallowGal" + @show nAdveVel +elseif a == 2 + Problem = "HaurwitzSphere" + Param = Examples.Parameters(FTB,Problem) + RadEarth = Phys.RadEarth + dtau = 30 #g=9.81, H=8000 + nAdveVel = ceil(Int,6*24*3600/dtau) + GridTypeOut = GridType*"NonLinShallowHaurwitz" + @show nAdveVel +elseif a == 3 + Problem = "LinearBlob" + Param = Examples.Parameters(FTB,Problem) + RadEarth = 1.0 + dtau = 0.00025 + nAdveVel = 100 + GridTypeOut = GridType*"NonLinShallowBlob" + @show nAdveVel +elseif a == 4 + Problem = "AdvectionSphereSpherical" + Param = Examples.Parameters(FTB,Problem) + RadEarth = 1.0 + dtau = 2*pi*RadEarth/4/nPanel/Param.uMax*0.1 + @show dtau # 0.0004581489286485114 #in s = 2*pi*Rad / 4*nPanel / param.uMax * cFL (ca. 0.7) bei RK2 (RK3 1.7) + nAdveVel = 100 + nprint = 10 + GridTypeOut = GridType*"Advec" + @show nAdveVel +else + print("Error") +end +println("The chosen Problem is ") +Examples.InitialProfile!(Model,Problem,Param,Phys) +ProfileUTR = Examples.AdvectionVelocity()(Param,Phys) + +#Grid construction +Grid, Exchange = Grids.InitGridSphere(backend,FTB,OrdPoly,nz,nPanel,RefineLevel,ns,nLat,nLon,LatB,GridType,Decomp,RadEarth, + Model,ParallelCom) +for iE = 1 : Grid.NumEdges + Grids.PosEdgeInFace!(Grid.Edges[iE],Grid.Edges,Grid.Faces) +end +vtkSkeletonMesh = Outputs.vtkStruct{Float64}(backend,Grid,Grid.NumFaces,Flat) + +#finite elements +VecDG = FEMSei.VecDG0Struct{FTB}(Grids.Tri(),backend,Grid) +RTTR = FEMSei.RT0Struct{FTB}(Grids.Tri(),backend,Grid) +RT = FEMSei.RT0Struct{FTB}(Grids.Tri(),backend,Grid) + +#massmatrix und LU-decomposition +VecDG.M = FEMSei.MassMatrix(backend,FTB,VecDG,Grid,nQuadM,FEMSei.Jacobi!) +VecDG.LUM = lu(VecDG.M) +RTTR.M = FEMSei.MassMatrix(backend,FTB,RTTR,Grid,nQuadM,FEMSei.Jacobi!) +RTTR.LUM = lu(RTTR.M) +RT.M = FEMSei.MassMatrix(backend,FTB,RT,Grid,nQuadM,FEMSei.Jacobi!) +RT.LUM = lu(RT.M) + +#stiffmatrix +cVecDG = zeros(FTB,VecDG.NumG) +cRT = zeros(FTB,RTTR.NumG) + +QuadOrd=3 + +FEMSei.Interpolate1!(backend,FTB,cRT,RTTR,Grid.Type,Grid,nQuad,FEMSei.Jacobi!,ProfileUTR) +#projection from +FEMSei.ProjectHDivVecDG!(backend,FTB,cVecDG,VecDG,cRT,RTTR,Grid,Grids.Quad(),nQuad,FEMSei.Jacobi!) + +u = zeros(FTB,RT.NumG) +VelSp = zeros(Grid.NumFaces,2) + +#Fortin-Interpolation +FEMSei.Interpolate1!(backend,FTB,u,RT,Grid.Type,Grid,QuadOrd,FEMSei.Jacobi!,Model.InitialProfile) + +FEMSei.ConvertVelocitySp!(backend,FTB,VelSp,cRT,RTTR,Grid,FEMSei.Jacobi!) + +FileNumber=0 +Outputs.vtkSkeleton!(vtkSkeletonMesh, "Proj_Recov_Tri", Proc, ProcNumber, VelSp, FileNumber) + +#runge-kutta steps +time = 0.0 +cRTNew = similar(cRT) +Rhs = zeros(FTB,RT.NumG) + +for i = 1 : nAdveVel + @show i + @. Rhs = 0 + FEMSei.ProjectHDivVecDG!(backend,FTB,cVecDG,VecDG,cRT,RTTR,Grid,Grids.Tri(),nQuad,FEMSei.Jacobi!) + #FEMSei.DivMomentumVector!(backend,FTB,Rhs,u,RT,cVecDG,VecDG,RT,Grid,Grids.Tri(),nQuad,FEMSei.Jacobi!) + FEMSei.DivMomentumVector!(backend,FTB,Rhs,RT,u,RT,cVecDG,VecDG,Grid,Grids.Tri(),nQuad,FEMSei.Jacobi!) + ldiv!(RTTR.LUM,Rhs) + ################TEST + # FEMSei.ConvertVelocitySp!(backend,FTB,VelSp,Rhs,RTTR,Grid,FEMSei.Jacobi!) + # FileNumber=4000 + # Outputs.vtkSkeleton!(vtkSkeletonMesh, "Rhs", Proc, ProcNumber, VelSp, FileNumber) + ################ + @. cRTNew = cRT + 0.5 * dtau * Rhs + @. Rhs = 0 + #CG nach vorn nach RHs sezten + FEMSei.ProjectHDivVecDG!(backend,FTB,cVecDG,VecDG,cRTNew,RTTR,Grid,Grids.Tri(),nQuad,FEMSei.Jacobi!) + #FEMSei.DivMomentumVector!(backend,FTB,Rhs,u,RT,cVecDG,VecDG,RT,Grid,Grids.Tri(),nQuad,FEMSei.Jacobi!) + FEMSei.DivMomentumVector!(backend,FTB,Rhs,RT,u,RT,cVecDG,VecDG,Grid,Grids.Tri(),nQuad,FEMSei.Jacobi!) + ldiv!(RTTR.LUM,Rhs) + @. cRT = cRT + dtau * Rhs + if mod(i,nprint) == 0 + FEMSei.ConvertVelocitySp!(backend,FTB,VelSp,cRT,RTTR,Grid,FEMSei.Jacobi!) + global FileNumber += 1 + Outputs.vtkSkeleton!(vtkSkeletonMesh, "Proj_Recov_Tri", Proc, ProcNumber,VelSp, FileNumber) + end +end +@show "finished" diff --git a/Jobs/NHCart/MountFujiCart b/Jobs/NHCart/MountFujiCart new file mode 100755 index 0000000..235edc4 --- /dev/null +++ b/Jobs/NHCart/MountFujiCart @@ -0,0 +1,50 @@ +#!/bin/bash +julia --project Examples/testNHCart.jl \ + --Problem="HillAgnesiXCart" \ + --NumberThreadGPU=512 \ + --FloatTypeBackend="Float32" \ + --NumV=5 \ + --NumTr=0 \ + --RefProfile=false \ + --Equation="CompressibleShallow" \ + --State="Dry" \ + --Buoyancy=true \ + --Coriolis=false \ + --Curl=true \ + --ModelType="VectorInvariant" \ + --SurfaceFluxMom=false \ + --VerticalDiffusionMom=false \ + --Damping=true \ + --StrideDamp=10000 \ + --Relax=1.0e-2 \ + --Upwind=true \ + --BoundaryWE="Period" \ + --BoundarySN="Period" \ + --BoundaryBT="" \ + --Thermo="" \ + --SimMinutes=30 \ + --PrintSeconds=60 \ + --PrintTime=0 \ + --dtau=.3 \ + --IntMethod="Rosenbrock" \ + --Table="SSP-Knoth" \ + --TopoS="MountFuji" \ + `#Topography ` \ + --P1=0.0\ + --P2=1000.0\ + --P3=400.0\ + --Stretch=true \ + --StretchType="ICON" \ + --H=19600.0 \ + --Lx=40000.0 \ + --Ly=4000.0 \ + --x0=-20000.0 \ + --y0=0.0 \ + --nz=80 \ + --nx=40 \ + --ny=2 \ + --OrdPoly=4 \ + --HyperVisc=true \ + --HyperDCurl=1.e6 \ + --HyperDGrad=1.e6 \ + --HyperDDiv=1.e6 diff --git a/src/Examples/parameters.jl b/src/Examples/parameters.jl index 662273d..8d02ab3 100755 --- a/src/Examples/parameters.jl +++ b/src/Examples/parameters.jl @@ -1,6 +1,6 @@ Base.@kwdef struct ParamGalewskiSphere H0G = 10000.0 - hH = 0.0 #120.0 OSWALD + hH = 120.0 alphaG = 1.0/3.0 betaG = 1.0/15.0 lat0G = pi/7.0 diff --git a/src/Examples/topography.jl b/src/Examples/topography.jl index 990f2d7..8e4bd67 100644 --- a/src/Examples/topography.jl +++ b/src/Examples/topography.jl @@ -1,9 +1,29 @@ abstract type Topography end +Base.@kwdef struct MountFuji{T} <: Topography + x0C::T = 0 + y0C::T = 0 + axC::T = 40000 + ayC::T = 40000 + hC::T = 3776 +end + +function (profile::MountFuji)() + (; x0C, y0C, axC, ayC, hC) = profile + function local_profile(x) + FT = eltype(x) + h = hC / (((x[1] -x0C) / axC)^2 + ((x[2] -y0C) / ayC)^2 + FT(1)); + return h + end + return local_profile +end + + + Base.@kwdef struct AgnesiHill{T} <: Topography x0C::T = 0 aC::T = 1000 - hC::T = 0 + hC::T = 400 end function (profile::AgnesiHill)() diff --git a/src/FEMSei/FEMSei.jl b/src/FEMSei/FEMSei.jl index e26711c..e90cb89 100755 --- a/src/FEMSei/FEMSei.jl +++ b/src/FEMSei/FEMSei.jl @@ -49,12 +49,12 @@ include("Jacobi.jl") include("StiffMatrix.jl") include("Stiff1.jl") include("Project.jl") -include("fp.jl") include("ConvertVelocity.jl") include("ModelFEM.jl") include("Fcn.jl") include("TimestepperFEM.jl") include("ConstructFEM.jl") +include("Interpolate.jl") diff --git a/src/FEMSei/Interpolate.jl b/src/FEMSei/Interpolate.jl new file mode 100644 index 0000000..71240c6 --- /dev/null +++ b/src/FEMSei/Interpolate.jl @@ -0,0 +1,152 @@ +function Interpolate1!(backend,FTB,uN,Fe::HDivConfElement,ElemType,Grid,QuadOrd,Jacobi,F) + NumQuadL, WeightsL, PointsL = QuadRule(Grids.Line(),QuadOrd) + X = zeros(3) + VelSp = zeros(3) + VelCa = zeros(3) + nR = zeros(3) + nL = zeros(3) + DF = zeros(3,2) + detDF = zeros(1) + pinvDF = zeros(3,2) + X = zeros(3) + + if ElemType == Grids.Quad() + PointsE = zeros(2,NumQuadL,4) + for iQ = 1 : NumQuadL + PointsE[1,iQ,1] = PointsL[iQ] + PointsE[2,iQ,1] = -1.0 + PointsE[1,iQ,2] = 1.0 + PointsE[2,iQ,2] = PointsL[iQ] + PointsE[1,iQ,3] = PointsL[iQ] + PointsE[2,iQ,3] = 1.0 + PointsE[1,iQ,4] = -1.0 + PointsE[2,iQ,4] = PointsL[iQ] + end + elseif ElemType == Grids.Tri() + PointsE = zeros(2,NumQuadL,3) + for iQ = 1 : NumQuadL + PointsE[1,iQ,1] = PointsL[iQ] + PointsE[2,iQ,1] = -1.0 + PointsE[1,iQ,2] = -PointsL[iQ] + PointsE[2,iQ,2] = PointsL[iQ] + PointsE[1,iQ,3] = -1.0 + PointsE[2,iQ,3] = PointsL[iQ] + end + end + @. uN = 0 + for iE = 1 : Grid.NumEdges + Edge = Grid.Edges[iE] + if length(Edge.F) > 1 + iFL = Edge.F[1] + EdgeTypeL = Edge.FE[1] + for iQ = 1 : NumQuadL + Jacobi!(DF,detDF,pinvDF,X,ElemType,PointsE[1,iQ,EdgeTypeL], + PointsE[2,iQ,EdgeTypeL],Grid.Faces[iFL], Grid) + nBarL = Grid.nBar[:, EdgeTypeL] + detDFL = detDF[1] * Grid.Faces[iFL].Orientation + nL[1] = (pinvDF[1,1] * nBarL[1] + pinvDF[1,2] * nBarL[2]) * detDFL + nL[2] = (pinvDF[2,1] * nBarL[1] + pinvDF[2,2] * nBarL[2]) * detDFL + nL[3] = (pinvDF[3,1] * nBarL[1] + pinvDF[3,2] * nBarL[2]) * detDFL + _,VelSp[1],VelSp[2],VelSp[3], = F(X,0.0) + lon,lat,r = Grids.cart2sphere(X[1],X[2],X[3]) + VelCa = VelSphere2Cart(VelSp,lon,lat) + ind = Fe.Glob[EdgeTypeL,iFL] + uN[iE] -= 0.5 * nL' * VelCa * WeightsL[iQ] + end + end + end +end + +function Interpolate!(backend,FTB,uN,Fe::HDivConfElement,ElemType,Grid,QuadOrd,Jacobi,F) + NumQuadL, WeightsL, PointsL = QuadRule(Grids.Line(),QuadOrd) + X = zeros(3) + VelSp = zeros(3) + VelCa = zeros(3) + n = zeros(3) + DF = zeros(3,2) + detDF = zeros(1) + pinvDF = zeros(3,2) + X = zeros(3) + + if ElemType == Grids.Quad() + PointsE = zeros(2,NumQuadL,4) + for iQ = 1 : NumQuadL + PointsE[1,iQ,1] = PointsL[iQ] + PointsE[2,iQ,1] = -1.0 + PointsE[1,iQ,2] = 1.0 + PointsE[2,iQ,2] = PointsL[iQ] + PointsE[1,iQ,3] = PointsL[iQ] + PointsE[2,iQ,3] = 1.0 + PointsE[1,iQ,4] = -1.0 + PointsE[2,iQ,4] = PointsL[iQ] + end + end + @. uN = 0 + for iE = 1 : Grid.NumEdges + Edge = Grid.Edges[iE] + if length(Edge.F) > 1 + iF = Edge.F[1] + EdgeType = Edge.FE[1] + for iQ = 1 : NumQuadL + Jacobi!(DF,detDF,pinvDF,X,ElemType,PointsE[1,iQ,EdgeType], + PointsE[2,iQ,EdgeType],Grid.Faces[iF], Grid) + nBar = Grid.nBar[:, EdgeType] + n[1] = (pinvDF[1,1] * nBar[1] + pinvDF[1,2] * nBar[2]) * detDF[1] + n[2] = (pinvDF[2,1] * nBar[1] + pinvDF[2,2] * nBar[2]) * detDF[1] + n[3] = (pinvDF[3,1] * nBar[1] + pinvDF[3,2] * nBar[2]) * detDF[1] + _,VelSp[1],VelSp[2],VelSp[3], = F(X,0.0) + lon,lat,r = Grids.cart2sphere(X[1],X[2],X[3]) + VelCa = VelSphere2Cart(VelSp,lon,lat) + ind = Fe.Glob[EdgeType,iF] + uN[iE] += - 0.5 * n' * VelCa * WeightsL[iQ] + end + end + end +end + +function InterpolateCons!(backend,FTB,uN,Fe::HDivConfElement,Grid,QuadOrd,Jacobi,F) + NumQuadL, WeightsL, PointsL = QuadRule(Grids.Line(),QuadOrd) + X = zeros(3) + VelSp = zeros(3) + for iE = 1 : Grid.NumEdges + X[1] = Grid.Edges[iE].Mid.x + X[2] = Grid.Edges[iE].Mid.y + X[3] = Grid.Edges[iE].Mid.z + h,VelSp[1],VelSp[2],VelSp[3], = F(X,0.0) + lon,lat,r = Grids.cart2sphere(X[1],X[2],X[3]) + VelCa = VelSphere2Cart(VelSp,lon,lat) * h + uN[iE] = -0.5 * Grid.Edges[iE].a * (Grid.Edges[iE].n.x * VelCa[1] + + Grid.Edges[iE].n.y * VelCa[2] + Grid.Edges[iE].n.z * VelCa[3]) + end +end + +function Project(backend,FTB,Fe::ScalarElement,Grid,QuadOrd,Jacobi,F) + NumQuad,Weights,Points = QuadRule(Fe.Type,QuadOrd) + fRef = zeros(Fe.Comp,Fe.DoF,length(Weights)) + + p=zeros(Fe.NumG) + for i = 1 : length(Weights) + for iComp = 1 : Fe.Comp + for iD = 1 : Fe.DoF + fRef[iComp,iD,i] = Fe.phi[iD,iComp](Points[i,1],Points[i,2]) + end + end + end + DF = zeros(3,2) + detDF = zeros(1) + pinvDF = zeros(3,2) + X = zeros(3) + for iF = 1 : Grid.NumFaces + pLoc = zeros(Fe.Comp,Fe.DoF) + for i = 1 : length(Weights) + Jacobi!(DF,detDF,pinvDF,X,Grid.Type,Points[i,1],Points[i,2],Grid.Faces[iF], Grid) + detDFLoc = detDF[1] + f, = F(X,0.0) + pLoc += abs(detDFLoc)*Weights[i]*(fRef[:,:,i]*f) + end + @. p[Fe.Glob[:,iF]] += pLoc[Fe.Comp,:] + end + ldiv!(Fe.LUM,p) + return p +end + diff --git a/src/FEMSei/Project.jl b/src/FEMSei/Project.jl index e0c6b5a..c1c4a79 100644 --- a/src/FEMSei/Project.jl +++ b/src/FEMSei/Project.jl @@ -1,4 +1,4 @@ - +#= function Interpolate1!(backend,FTB,uN,Fe::HDivConfElement,ElemType,Grid,QuadOrd,Jacobi,F) NumQuadL, WeightsL, PointsL = QuadRule(Grids.Line(),QuadOrd) X = zeros(3) @@ -150,6 +150,7 @@ function Project(backend,FTB,Fe::ScalarElement,Grid,QuadOrd,Jacobi,F) ldiv!(Fe.LUM,p) return p end +=# function Project!(backend,FTB,p,Fe::ScalarElement,Grid,QuadOrd,Jacobi,F) NumQuad,Weights,Points = QuadRule(Fe.Type,QuadOrd) diff --git a/src/FEMSei/StiffMatrix.jl b/src/FEMSei/StiffMatrix.jl index 26185f2..c7f48b1 100644 --- a/src/FEMSei/StiffMatrix.jl +++ b/src/FEMSei/StiffMatrix.jl @@ -1903,3 +1903,277 @@ function DivMomentumScalar!(backend,FTB,Rhs,uHDiv,FeHDiv::HDivElement, end end end + +#vector for Tri +function DivMomentumVector!(backend,FTB,Rhs,FeTHDiv::HDivElement,uHDiv,FeHDiv::HDivElement, + uVecDG,FeVecDG::VectorElement,Grid,ElemType::Grids.ElementType,QuadOrd,Jacobi) + + NumQuad, Weights, Points = QuadRule(ElemType,QuadOrd) + fuHDiv = zeros(FeHDiv.DoF,FeHDiv.Comp,NumQuad) + fuVecDG = zeros(FeVecDG.DoF,FeVecDG.DoF,NumQuad) + fuTHDiv = zeros(FeTHDiv.DoF,FeTHDiv.Comp,NumQuad) + fuDivHDiv = zeros(FeHDiv.DoF,NumQuad) + fuGradVecDG = zeros(FeVecDG.DoF,FeVecDG.DoF,2,NumQuad) + + #computation of the ansatz functions in the quadrature points + for iQ = 1 : NumQuad + for iComp = 1 : FeTHDiv.Comp + for iD = 1 : FeTHDiv.DoF + fuTHDiv[iD,iComp,iQ] = FeTHDiv.phi[iD,iComp](Points[iQ,1],Points[iQ,2]) + end + end + for iComp = 1 : FeTHDiv.Comp + for iD = 1 : FeHDiv.DoF + fuHDiv[iD,iComp,iQ] = FeHDiv.phi[iD,iComp](Points[iQ,1],Points[iQ,2]) + end + end + for iD = 1 : FeVecDG.DoF + fuVecDG[iD,1,iQ] = FeVecDG.phi[iD,1](Points[iQ,1],Points[iQ,2]) + fuVecDG[iD,2,iQ] = FeVecDG.phi[iD,2](Points[iQ,1],Points[iQ,2]) + fuVecDG[iD,3,iQ] = FeVecDG.phi[iD,3](Points[iQ,1],Points[iQ,2]) + end + for iD = 1 : FeHDiv.DoF + fuDivHDiv[iD,iQ] = FeHDiv.Divphi[iD,1](Points[iQ,1],Points[iQ,2]) + end + + for iD = 1 : FeVecDG.DoF + fuGradVecDG[iD,1,1,iQ] = FeVecDG.Gradphi[iD,1,1](Points[iQ,1],Points[iQ,2]) + fuGradVecDG[iD,1,2,iQ] = FeVecDG.Gradphi[iD,1,2](Points[iQ,1],Points[iQ,2]) + fuGradVecDG[iD,2,1,iQ] = FeVecDG.Gradphi[iD,2,1](Points[iQ,1],Points[iQ,2]) + fuGradVecDG[iD,2,2,iQ] = FeVecDG.Gradphi[iD,2,2](Points[iQ,1],Points[iQ,2]) + fuGradVecDG[iD,3,1,iQ] = FeVecDG.Gradphi[iD,3,1](Points[iQ,1],Points[iQ,2]) + fuGradVecDG[iD,3,2,iQ] = FeVecDG.Gradphi[iD,3,2](Points[iQ,1],Points[iQ,2]) + end + end + + + #local variables + uVecDGLoc = zeros(FeVecDG.DoF) + uuVecDGLoc = zeros(3) + uuGradVecDGLoc = zeros(3,2) + uHDivLoc = zeros(FeHDiv.DoF) + uuHDivLoc = zeros(2) + RhsLoc = zeros(FeTHDiv.DoF) + + GradVecDGHDiv = zeros(3) + VecDGDivHDiv = zeros(3) + + #Jacobi-Allocation + DF = zeros(3,2) + detDF = zeros(1) + pinvDF = zeros(3,2) + X = zeros(3) + + DFL = zeros(3,2) + detDFL = zeros(1) + pinvDFL = zeros(3,2) + XL = zeros(3) + + DFR = zeros(3,2) + detDFR = zeros(1) + pinvDFR = zeros(3,2) + XR = zeros(3) + innersum = zeros(3) + + #copy from global variables to local variables + @inbounds for iF = 1 : Grid.NumFaces + @inbounds for iD = 1 : FeVecDG.DoF + ind = FeVecDG.Glob[iD,iF] + uVecDGLoc[iD] = uVecDG[ind] + end + @inbounds for iD = 1 : FeHDiv.DoF + ind = FeHDiv.Glob[iD,iF] + uHDivLoc[iD] = uHDiv[ind] + end + @. RhsLoc = 0.0 + + @inbounds for iQ = 1 : NumQuad + @. uuVecDGLoc = 0.0 + @. uuGradVecDGLoc = 0.0 + #computation of local variables in a quadrature point + @inbounds for iD = 1 : FeVecDG.DoF + @inbounds for i = 1 : 3 + uuVecDGLoc[i] += uVecDGLoc[iD] * fuVecDG[iD,i,iQ] + @inbounds for j = 1 : 2 + uuGradVecDGLoc[i,j] += uVecDGLoc[iD] * fuGradVecDG[iD,i,j,iQ] + end + end + end + uuDivHDivLoc = 0.0 + @. uuHDivLoc = 0.0 + @inbounds for iD = 1 : FeHDiv.DoF + uuDivHDivLoc += uHDivLoc[iD] * fuDivHDiv[iD,iQ] + @inbounds for i = 1 : 2 + uuHDivLoc[i] += uHDivLoc[iD] * fuHDiv[iD,i,iQ] + end + end + + @inbounds for i = 1 : 3 + GradVecDGHDiv[i] = 0.0 + VecDGDivHDiv[i] = uuDivHDivLoc * uuVecDGLoc[i] + @inbounds for j = 1 : 2 + GradVecDGHDiv[i] += uuGradVecDGLoc[i,j] * uuHDivLoc[j] + end + end + @. innersum = Grid.Faces[iF].Orientation * (VecDGDivHDiv + GradVecDGHDiv) + #computation of Jacobi + Jacobi!(DF,detDF,pinvDF,X,Grid.Type,Points[iQ,1],Points[iQ,2],Grid.Faces[iF], Grid) + detDFLoc = detDF[1] + product = 0.0 + @inbounds for iD = 1 : FeTHDiv.DoF + @inbounds for i = 1 : 3 + temp = 0.0 + @inbounds for j = 1 : 2 + temp += DF[i,j] * fuTHDiv[iD,j,iQ] + end + product += innersum[i] * temp + end + end + end + @inbounds for iD = 1 : FeTHDiv.DoF + ind = FeTHDiv.Glob[iD,iF] + Rhs[ind] += RhsLoc[iD] + end + end + + #UPWIND on Edges + NumQuadL, WeightsL, PointsL = QuadRule(Grids.Line(),QuadOrd) + uFFRef = zeros(FeHDiv.DoF, FeHDiv.Comp,NumQuadL,3) + uVecDGRef = zeros(FeVecDG.DoF,FeVecDG.Comp,NumQuadL,3) + fTRef = zeros(FeTHDiv.DoF,FeTHDiv.Comp,NumQuadL,3) + uFFLocL = zeros(FeHDiv.DoF) + uFFLocR = zeros(FeHDiv.DoF) + uHDivLocLeft = zeros(FeHDiv.DoF) + uHDivLocRight = zeros(FeHDiv.DoF) + uVecDGLocLeft = zeros(FeVecDG.DoF) + uVecDGLocRight = zeros(FeVecDG.DoF) + + for iQ = 1 : NumQuadL + for iComp = 1 : FeTHDiv.Comp + for iD = 1 : FeTHDiv.DoF + fTRef[iD,iComp,iQ,1] = FeTHDiv.phi[iD,iComp](PointsL[iQ],-1.0) + fTRef[iD,iComp,iQ,2] = FeTHDiv.phi[iD,iComp](-PointsL[iQ],PointsL[iQ]) + fTRef[iD,iComp,iQ,3] = FeTHDiv.phi[iD,iComp](-1.0,PointsL[iQ]) + end + end + #uf RT UHDiv + for iComp = 1 : FeHDiv.Comp + for iD = 1 : FeHDiv.DoF + uFFRef[iD,iComp,iQ,1] = FeHDiv.phi[iD,iComp](PointsL[iQ],-1.0) + uFFRef[iD,iComp,iQ,2] = FeHDiv.phi[iD,iComp](-PointsL[iQ],PointsL[iQ]) + uFFRef[iD,iComp,iQ,3] = FeHDiv.phi[iD,iComp](-1.0,PointsL[iQ]) + end + end + #hf + for iComp = 1 : FeVecDG.Comp + for iD = 1 : FeVecDG.DoF + uVecDGRef[iD,iComp,iQ,1] = FeVecDG.phi[iD,iComp](PointsL[iQ],-1.0) + uVecDGRef[iD,iComp,iQ,2] = FeVecDG.phi[iD,iComp](-PointsL[iQ],PointsL[iQ]) + uVecDGRef[iD,iComp,iQ,3] = FeVecDG.phi[iD,iComp](-1.0,PointsL[iQ]) + end + end + end + + #allocation Mloc + MLoc11 = zeros(FeTHDiv.DoF, FeHDiv.DoF) + MLoc12 = zeros(FeTHDiv.DoF, FeHDiv.DoF) + MLoc21 = zeros(FeTHDiv.DoF, FeHDiv.DoF) + MLoc22 = zeros(FeTHDiv.DoF, FeHDiv.DoF) + fTLocL = zeros(3,FeTHDiv.DoF) + fTLocR = zeros(3,FeTHDiv.DoF) + + uL = zeros(3) + uR = zeros(3) + + @inbounds for iE in 1:Grid.NumEdges + Edge = Grid.Edges[iE] + if length(Edge.F) > 1 + iFL = Edge.F[1] + EdgeTypeL = Edge.FE[1] + iFR = Edge.F[2] + EdgeTypeR = Edge.FE[2] + #computation normales of edges + @views nBarLocL = Grid.nBar[:, EdgeTypeL] + @views nBarLocR = Grid.nBar[:, EdgeTypeR] + #gamma upwind value + uE = 0.0 + @inbounds for iD = 1 : FeHDiv.DoF + ind = FeHDiv.Glob[iD,iFL] + uHDivLocLeft[iD] = uHDiv[ind] + ind = FeHDiv.Glob[iD,iFR] + uHDivLocRight[iD] = uHDiv[ind] + end + @inbounds for iD = 1 : FeVecDG.DoF + ind = FeVecDG.Glob[iD,iFL] + uVecDGLocLeft[iD] = uVecDG[ind] + ind = FeVecDG.Glob[iD,iFR] + uVecDGLocRight[iD] = uVecDG[ind] + end + @inbounds for iQ = 1:length(WeightsL) + @inbounds for iD = 1 : FeHDiv.DoF + @views uE += WeightsL[iQ]*(uHDivLocLeft[iD]*(uFFRef[iD,:,iQ,EdgeTypeL]'*nBarLocL)+uHDivLocRight[iD]*(uFFRef[iD,:,iQ,EdgeTypeR]'*nBarLocR)) + end + end + #upwind value + gammaU = 0.5 + gammaLoc = uE > 0 ? gammaU : -gammaU + + @. MLoc11 = 0 + @. MLoc12 = 0 + @. MLoc21 = 0 + @. MLoc22 = 0 + + @inbounds for iQ in 1:length(WeightsL) + #computation of Jacobi EdgetypeL + Jacobi!(DFL,detDFL,pinvDFL,XL,Grid.Type,Points[iQ,1],Points[iQ,2],Grid.Faces[iFL], Grid) + detDFLLoc = detDFL[1] * Grid.Faces[iFL].Orientation + #EdgeTypeR + Jacobi!(DFR,detDFR,pinvDFR,XR,Grid.Type,Points[iQ,1],Points[iQ,2],Grid.Faces[iFR], Grid) + detDFRLoc = detDFR[1] * Grid.Faces[iFR].Orientation + + @inbounds for iDoF = 1 : FeTHDiv.DoF + @inbounds for j = 1 : 3 + fTLocL[j,iDoF] = (1/detDFLLoc) * (DFL[j,1] * fTRef[iDoF,1,iQ,EdgeTypeL] + + DFL[j,2] * fTRef[iDoF,2,iQ,EdgeTypeL]) + fTLocR[j,iDoF] = (1/detDFRLoc) * (DFR[j,1] * fTRef[iDoF,1,iQ,EdgeTypeR] + + DFR[j,2] * fTRef[iDoF,2,iQ,EdgeTypeR]) + end + end + @inbounds for iD in 1:FeHDiv.DoF + @views uFFLocL[iD] = uFFRef[iD,:,iQ,EdgeTypeL]' * nBarLocL + @views uFFLocR[iD] = uFFRef[iD,:,iQ,EdgeTypeR]' * nBarLocR + end + @. uL = 0.0 + @. uR = 0.0 + + @inbounds for iD in 1:FeVecDG.DoF + uL[1] += uVecDGRef[iD,1,iQ,EdgeTypeL] * uVecDGLocLeft[iD] + uL[2] += uVecDGRef[iD,2,iQ,EdgeTypeL] * uVecDGLocLeft[iD] + uL[3] += uVecDGRef[iD,3,iQ,EdgeTypeL] * uVecDGLocLeft[iD] + uR[1] += uVecDGRef[iD,1,iQ,EdgeTypeR] * uVecDGLocRight[iD] + uR[2] += uVecDGRef[iD,2,iQ,EdgeTypeR] * uVecDGLocRight[iD] + uR[3] += uVecDGRef[iD,3,iQ,EdgeTypeR] * uVecDGLocRight[iD] + end + @inbounds for iDoF = 1 : FeTHDiv.DoF + @inbounds for jDoF = 1 : FeHDiv.DoF + @inbounds for k = 1 : 3 + MLoc11[iDoF,jDoF] += WeightsL[iQ] * (fTLocL[k,iDoF]' * uL[k]) * uFFLocL[jDoF] + MLoc12[iDoF,jDoF] += WeightsL[iQ] * (fTLocL[k,iDoF]' * uR[k]) * uFFLocR[jDoF] + MLoc21[iDoF,jDoF] += WeightsL[iQ] * (fTLocR[k,iDoF]' * uL[k]) * uFFLocL[jDoF] + MLoc22[iDoF,jDoF] += WeightsL[iQ] * (fTLocR[k,iDoF]' * uR[k]) * uFFLocR[jDoF] + end + end + end + end + #all together over edges + @inbounds for iD = 1 : FeTHDiv.DoF + indL = FeTHDiv.Glob[iD,iFL] + @views Rhs[indL] += (-0.5 - gammaLoc) * MLoc11[iD,:]' * uHDivLocLeft + @views Rhs[indL] += (-0.5 + gammaLoc) * MLoc12[iD,:]' * uHDivLocRight + indR = FeTHDiv.Glob[iD,iFR] + @views Rhs[indR] += (+0.5 + gammaLoc) * MLoc21[iD,:]' * uHDivLocLeft + @views Rhs[indR] += (+0.5 - gammaLoc) * MLoc22[iD,:]' * uHDivLocRight + end + end + end +end diff --git a/src/FEMSei/VecDG (1).jl b/src/FEMSei/VecDG (1).jl new file mode 100644 index 0000000..1916722 --- /dev/null +++ b/src/FEMSei/VecDG (1).jl @@ -0,0 +1,313 @@ +mutable struct VecDG0Struct{FT<:AbstractFloat, + IT2<:AbstractArray} <: VectorElement + Glob::IT2 + DoF::Int + Comp::Int + phi::Array{Polynomial,2} + Gradphi::Array{Polynomial,3} + NumG::Int + NumI::Int + Type::Grids.ElementType + M::AbstractSparseMatrix + LUM::SparseArrays.UMFPACK.UmfpackLU{Float64, Int64} +end + +#VecDG0 Quad + +function VecDG0Struct{FT}(::Grids.Quad,backend,Grid) where FT<:AbstractFloat + Glob = KernelAbstractions.zeros(backend,Int,0,0) + Type = Grids.Quad() + DoF = 3 #vorn + Comp = 3 + @polyvar x1 x2 + phi = Array{Polynomial,2}(undef,DoF,Comp) + phi[1,1] = 1.0 + 0.0*x1 + 0.0*x2 + phi[1,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[1,3] = 0.0 + 0.0*x1 + 0.0*x2 + phi[2,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[2,2] = 1.0 + 0.0*x1 + 0.0*x2 + phi[2,3] = 0.0 + 0.0*x1 + 0.0*x2 + phi[3,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[3,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[3,3] = 1.0 + 0.0*x1 + 0.0*x2 + + Gradphi = Array{Polynomial,3}(undef,DoF,Comp,2) + for i = 1 : DoF + for j = 1 : Comp + Gradphi[i,j,1] = differentiate(phi[i,j],x1) + Gradphi[i,j,2] = differentiate(phi[i,j],x2) + end + end + + Glob = KernelAbstractions.zeros(backend,Int,DoF,Grid.NumFaces) + GlobCPU = zeros(Int,DoF,Grid.NumFaces) + NumG = DoF * Grid.NumFaces + NumI = DoF * Grid.NumFaces + for iF = 1 : Grid.NumFaces + for iD = 1 : DoF + GlobCPU[iD,iF] = (Grid.Faces[iF].F - 1 ) * DoF + iD + end + end + copyto!(Glob,GlobCPU) + M = sparse([1],[1],[1.0]) + LUM = lu(M) + return VecDG0Struct{FT, + typeof(Glob)}( + Glob, + DoF, + Comp, + phi, + Gradphi, + NumG, + NumI, + Type, + M, + LUM, + ) +end + +#VecDG0 Tri + +function VecDG0Struct{FT}(::Grids.Tri,backend,Grid) where FT<:AbstractFloat + Glob = KernelAbstractions.zeros(backend,Int,0,0) + Type = Grids.Tri() + DoF = 3 + Comp = 3 + @polyvar x1 x2 + phi = Array{Polynomial,2}(undef,DoF,Comp) + phi[1,1] = 1.0 + 0.0*x1 + 0.0*x2 + phi[1,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[1,3] = 0.0 + 0.0*x1 + 0.0*x2 + phi[2,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[2,2] = 1.0 + 0.0*x1 + 0.0*x2 + phi[2,3] = 0.0 + 0.0*x1 + 0.0*x2 + phi[3,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[3,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[3,3] = 1.0 + 0.0*x1 + 0.0*x2 + + Gradphi = Array{Polynomial,3}(undef,DoF,Comp,2) + for i = 1 : DoF + for j = 1 : Comp + Gradphi[i,j,1] = differentiate(phi[i,j],x1) + Gradphi[i,j,2] = differentiate(phi[i,j],x2) + end + end + + Glob = KernelAbstractions.zeros(backend,Int,DoF,Grid.NumFaces) + GlobCPU = zeros(Int,DoF,Grid.NumFaces) + NumG = DoF * Grid.NumFaces + NumI = DoF * Grid.NumFaces + for iF = 1 : Grid.NumFaces + for iD = 1 : DoF + GlobCPU[iD,iF] = (Grid.Faces[iF].F - 1 ) * DoF + iD + end + end + copyto!(Glob,GlobCPU) + M = sparse([1],[1],[1.0]) + LUM = lu(M) + return VecDG0Struct{FT, + typeof(Glob)}( + Glob, + DoF, + Comp, + phi, + Gradphi, + NumG, + NumI, + Type, + M, + LUM, + ) +end + +mutable struct VecDG1Struct{FT<:AbstractFloat, + IT2<:AbstractArray} <: VectorElement + Glob::IT2 + DoF::Int + Comp::Int + phi::Array{Polynomial,2} + Gradphi::Array{Polynomial,3} + NumG::Int + NumI::Int + Type::Grids.ElementType + M::AbstractSparseMatrix + LUM::SparseArrays.UMFPACK.UmfpackLU{Float64, Int64} +end + +#VecDG1 Quad + +function VecDG1Struct{FT}(::Grids.Quad,backend,Grid) where FT<:AbstractFloat + Glob = KernelAbstractions.zeros(backend,Int,0,0) + Type = Grids.Quad() + DoF = 12 #vorn + Comp = 3 + @polyvar x1 x2 + phi = Array{Polynomial,2}(undef,DoF,Comp) + phi[1,1] = (1.0-1.0*x1) * (1.0-1.0*x2) + phi[1,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[1,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[2,1] = (1.0+1.0*x1) * (1.0-1.0*x2) + phi[2,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[2,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[3,1] = (1.0+1.0*x1) * (1.0+1.0*x2) + phi[3,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[3,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[4,1] = (1.0-1.0*x1) * (1.0+1.0*x2) + phi[4,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[4,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[5,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[5,2] = (1.0-1.0*x1) * (1.0-1.0*x2) + phi[5,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[6,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[6,2] = (1.0+1.0*x1) * (1.0-1.0*x2) + phi[6,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[7,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[7,2] = (1.0+1.0*x1) * (1.0+1.0*x2) + phi[7,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[8,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[8,2] = (1.0-1.0*x1) * (1.0+1.0*x2) + phi[8,3] = 0.0 + 0.0*x1 + 0.0*x2 + + phi[9,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[9,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[9,3] = (1.0-1.0*x1) * (1.0-1.0*x2) + + phi[10,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[10,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[10,3] = (1.0+1.0*x1) * (1.0-1.0*x2) + + phi[11,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[11,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[11,3] = (1.0+1.0*x1) * (1.0+1.0*x2) + + phi[12,1] = 0.0 + 0.0*x1 + 0.0*x2 + phi[12,2] = 0.0 + 0.0*x1 + 0.0*x2 + phi[12,3] = (1.0-1.0*x1) * (1.0+1.0*x2) + + Gradphi = Array{Polynomial,3}(undef,DoF,Comp,2) + for i = 1 : DoF + for j = 1 : Comp + Gradphi[i,j,1] = differentiate(phi[i,j],x1) + Gradphi[i,j,2] = differentiate(phi[i,j],x2) + end + end + + Glob = KernelAbstractions.zeros(backend,Int,DoF,Grid.NumFaces) + GlobCPU = zeros(Int,DoF,Grid.NumFaces) + NumG = DoF * Grid.NumFaces + NumI = DoF * Grid.NumFaces + for iF = 1 : Grid.NumFaces + for iD = 1 : DoF + GlobCPU[iD,iF] = (Grid.Faces[iF].F - 1 ) * DoF + iD + end + end + copyto!(Glob,GlobCPU) + M = sparse([1],[1],[1.0]) + LUM = lu(M) + return VecDG1Struct{FT, + typeof(Glob)}( + Glob, + DoF, + Comp, + phi, + Gradphi, + NumG, + NumI, + Type, + M, + LUM, + ) +end + +#VecDG1 Tri + +function VecDG1Struct{FT}(::Grids.Tri,backend,Grid) where FT<:AbstractFloat + Glob = KernelAbstractions.zeros(backend,Int,0,0) + Type = Grids.Tri() + DoF = 9 #vorn + Comp = 3 + @polyvar ksi1 ksi2 + phi = Array{Polynomial,2}(undef,DoF,Comp) + phi[1,1] = 0.0*ksi1 + 0.0*ksi2 + 2.0 + phi[1,2] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[1,3] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + + phi[2,1] = 0.0*ksi1 + 0.0*ksi2 + 2.0 + phi[2,2] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[2,3] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + + phi[3,1] = 0.0*ksi1 + 0.0*ksi2 + 2.0 + phi[3,2] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[3,3] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + + phi[4,1] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[4,2] = 0.0*ksi1 + 1.0*ksi2 - 1.0/3.0 + phi[4,3] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + + phi[5,1] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[5,2] = 0.0*ksi1 + 1.0*ksi2 - 1.0/3.0 + phi[5,3] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + + phi[6,1] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[6,2] = 0.0*ksi1 + 1.0*ksi2 - 1.0/3.0 + phi[6,3] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + + phi[7,1] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[7,2] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[7,3] = 1.0*ksi1 + 0.0*ksi2 - 1.0/3.0 + + phi[8,1] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[8,2] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[8,3] = 1.0*ksi1 + 0.0*ksi2 - 1.0/3.0 + + phi[9,1] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[9,2] = 0.0 + 0.0*ksi1 + 0.0*ksi2 + phi[9,3] = 1.0*ksi1 + 0.0*ksi2 - 1.0/3.0 + + for s = 1 : DoF + for t = 1 : Comp + phi[s,t] = subs(nu[s,t], ksi1 => (x1+1)/2, ksi2 => (x2+1)/2) + end + end + + Gradphi = Array{Polynomial,3}(undef,DoF,Comp,2) + for i = 1 : DoF + for j = 1 : Comp + Gradphi[i,j,1] = differentiate(phi[i,j],x1) + Gradphi[i,j,2] = differentiate(phi[i,j],x2) + end + end + + Glob = KernelAbstractions.zeros(backend,Int,DoF,Grid.NumFaces) + GlobCPU = zeros(Int,DoF,Grid.NumFaces) + NumG = DoF * Grid.NumFaces + NumI = DoF * Grid.NumFaces + for iF = 1 : Grid.NumFaces + for iD = 1 : DoF + GlobCPU[iD,iF] = (Grid.Faces[iF].F - 1 ) * DoF + iD + end + end + copyto!(Glob,GlobCPU) + M = sparse([1],[1],[1.0]) + LUM = lu(M) + return VecDG1Struct{FT, + typeof(Glob)}( + Glob, + DoF, + Comp, + phi, + Gradphi, + NumG, + NumI, + Type, + M, + LUM, + ) +end \ No newline at end of file diff --git a/src/Grids/CartGrid.jl b/src/Grids/CartGrid.jl index 7106d7f..370d46a 100644 --- a/src/Grids/CartGrid.jl +++ b/src/Grids/CartGrid.jl @@ -226,6 +226,7 @@ function CartGrid(backend,FT,nx::Int,ny::Int,lx::Float64,ly::Float64,x0::Float64 nBar3 = zeros(0,0) nBar = zeros(0,0) AdaptGrid = "" + Rad = 1.0 return GridStruct{FT, typeof(z)}( nz, diff --git a/src/Grids/GridLength.jl b/src/Grids/GridLength.jl index d54c249..2c94ac8 100644 --- a/src/Grids/GridLength.jl +++ b/src/Grids/GridLength.jl @@ -1,9 +1,12 @@ function GridLength(Grid) Faces = Grid.Faces - LengthLoc = 0 + LengthMaxLoc = 0 + LengthMinLoc = 4 * pi * Grid.Rad^2 for iF = 1 : Grid.NumFaces - LengthLoc = max(sqrt(Faces[iF].Area), LengthLoc) + LengthMaxLoc = max(sqrt(Faces[iF].Area), LengthMaxLoc) + LengthMinLoc = min(sqrt(Faces[iF].Area), LengthMinLoc) end - Length = MPI.Allreduce(LengthLoc, MPI.MAX, MPI.COMM_WORLD) - return Length + LengthMax = MPI.Allreduce(LengthMaxLoc, MPI.MAX, MPI.COMM_WORLD) + LengthMin = MPI.Allreduce(LengthMinLoc, MPI.MIN, MPI.COMM_WORLD) + return LengthMin,LengthMax end diff --git a/src/Grids/Renumbering.jl b/src/Grids/Renumbering.jl index 2052e80..c2b6d9a 100644 --- a/src/Grids/Renumbering.jl +++ b/src/Grids/Renumbering.jl @@ -31,17 +31,19 @@ function PosEdgeInFace!(Edge,Edges,Faces) end if size(Edge.F,1)>1 iF = Edge.F[1] - EdgeType = Edge.FE[1] - if Faces[iF].Orientation * Faces[iF].OrientE[EdgeType] == -1 - iTemp = Edge.FE[1] - Edge.FE[1] = Edge.FE[2] - Edge.FE[2] = iTemp - iTemp = Edge.F[1] - Edge.F[1] = Edge.F[2] - Edge.F[2] = iTemp - iTemp = Edge.FG[1] - Edge.FG[1] = Edge.FG[2] - Edge.FG[2] = iTemp + if iF > 0 + EdgeType = Edge.FE[1] + if Faces[iF].Orientation * Faces[iF].OrientE[EdgeType] == -1 + iTemp = Edge.FE[1] + Edge.FE[1] = Edge.FE[2] + Edge.FE[2] = iTemp + iTemp = Edge.F[1] + Edge.F[1] = Edge.F[2] + Edge.F[2] = iTemp + iTemp = Edge.FG[1] + Edge.FG[1] = Edge.FG[2] + Edge.FG[2] = iTemp + end end end end diff --git a/src/Grids/SubGridGhost.jl b/src/Grids/SubGridGhost.jl index 697426d..e4ba58c 100644 --- a/src/Grids/SubGridGhost.jl +++ b/src/Grids/SubGridGhost.jl @@ -1,5 +1,4 @@ function ConstructSubGridGhost(GlobalGrid,Proc,ProcNumber;order=true) - @show "ConstructSubGridGhost" backend = get_backend(GlobalGrid.z) FT = eltype(GlobalGrid.z) diff --git a/src/Outputs/vtkSphere.jl b/src/Outputs/vtkSphere.jl old mode 100755 new mode 100644 index 6437eb3..03e1e4b --- a/src/Outputs/vtkSphere.jl +++ b/src/Outputs/vtkSphere.jl @@ -324,7 +324,7 @@ function vtkSkeleton!(vtkCache,filename, part::Int, nparts::Int, c, FileNumber) step = FileNumber stepS = "$step" - vtk_filename_noext = filename * stepS; + vtk_filename_noext = pwd()*"/output/VTK/" * filename * stepS; vtk = pvtk_grid(vtk_filename_noext, pts, cells; compress=3, part = part, nparts = nparts) cName=["Height","uC","vC","wC","uS","vS"] for iC = 1 : size(c,2)