Merge branch 'main' into sk/embed_domain_into_mesh

CliMA · Sep 23, 2021 · 8c22aa0 · 8c22aa0
2 parents a121b41 + 10923dd
commit 8c22aa0
Show file tree

Hide file tree

Showing 40 changed files with 766 additions and 496 deletions.
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -177,6 +177,17 @@ steps:
       queue: central
       slurm_ntasks: 1
 
+  - label: ":computer: Bickley jet CG unstructured mesh"
+    key: "cpu_bickleyjet_cg_unsmesh"
+    command:
+      - "julia --color=yes --project=examples/bickleyjet examples/bickleyjet/bickleyjet_cg_unsmesh.jl"
+    artifact_paths:
+      - "examples/bickleyjet/output/cg_unsmesh/*"
+    agents:
+      config: cpu
+      queue: central
+      slurm_ntasks: 1
+
   - label: ":computer: Bickley jet CG vector invariant"
     key: "cpu_bickleyjet_cg_invariant"
     command:
@@ -237,7 +248,7 @@ steps:
     command:
       - "julia --color=yes --project=examples/hybrid examples/hybrid/bubble.jl"
     artifact_paths:
-      - "examples/hybrid/output/*"
+      - "examples/hybrid/output/bubble/*"
     agents:
       config: cpu
       queue: central

diff --git a/benchmarks/bickleyjet/bickleyjet_dg.jl b/benchmarks/bickleyjet/bickleyjet_dg.jl
@@ -21,8 +21,8 @@ const parameters = (
 
 
 domain = Domains.RectangleDomain(
-    -2π..2π,
-    -2π..2π,
+    Geometry.XPoint(-2π)..Geometry.XPoint(2π),
+    Geometry.YPoint(-2π)..Geometry.YPoint(2π),
     x1periodic = true,
     x2periodic = true,
 )

diff --git a/examples/bickleyjet/bickleyjet_cg.jl b/examples/bickleyjet/bickleyjet_cg.jl
@@ -24,8 +24,8 @@ const parameters = (
 
 
 domain = Domains.RectangleDomain(
-    -2π..2π,
-    -2π..2π,
+    Geometry.XPoint(-2π)..Geometry.XPoint(2π),
+    Geometry.YPoint(-2π)..Geometry.YPoint(2π),
     x1periodic = true,
     x2periodic = true,
 )

diff --git a/examples/bickleyjet/bickleyjet_cg_invariant.jl b/examples/bickleyjet/bickleyjet_cg_invariant.jl
@@ -22,8 +22,8 @@ const parameters = (
 )
 
 domain = Domains.RectangleDomain(
-    -2π..2π,
-    -2π..2π,
+    Geometry.XPoint(-2π)..Geometry.XPoint(2π),
+    Geometry.YPoint(-2π)..Geometry.YPoint(2π),
     x1periodic = true,
     x2periodic = true,
 )

diff --git a/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl b/examples/bickleyjet/bickleyjet_cg_invariant_hypervisc.jl
@@ -23,8 +23,8 @@ const parameters = (
 )
 
 domain = Domains.RectangleDomain(
-    -2π..2π,
-    -2π..2π,
+    Geometry.XPoint(-2π)..Geometry.XPoint(2π),
+    Geometry.YPoint(-2π)..Geometry.YPoint(2π),
     x1periodic = true,
     x2periodic = true,
 )

diff --git a/examples/bickleyjet/bickleyjet_cg_unsmesh.jl b/examples/bickleyjet/bickleyjet_cg_unsmesh.jl
@@ -24,8 +24,8 @@ const parameters = (
 )
 
 domain = Domains.RectangleDomain(
-    -2π..2π,
-    -2π..2π,
+    Geometry.XPoint(-2π)..Geometry.XPoint(2π),
+    Geometry.YPoint(-2π)..Geometry.YPoint(2π),
     x1periodic = true,
     x2periodic = true,
 )
@@ -124,7 +124,7 @@ ENV["GKSwstype"] = "nul"
 import Plots
 Plots.GRBackend()
 
-dirname = "cg_Uns2DMesh"
+dirname = "cg_unsmesh"
 path = joinpath(@__DIR__, "output", dirname)
 mkpath(path)
 

diff --git a/examples/bickleyjet/bickleyjet_dg.jl b/examples/bickleyjet/bickleyjet_dg.jl
@@ -30,8 +30,8 @@ numflux_name = get(ARGS, 1, "rusanov")
 boundary_name = get(ARGS, 2, "")
 
 domain = Domains.RectangleDomain(
-    -2π..2π,
-    -2π..2π,
+    Geometry.XPoint(-2π)..Geometry.XPoint(2π),
+    Geometry.YPoint(-2π)..Geometry.YPoint(2π),
     x1periodic = true,
     x2periodic = boundary_name != "noslip",
     x2boundary = boundary_name != "noslip" ? nothing : (:south, :north),

diff --git a/examples/column/advect.jl b/examples/column/advect.jl
@@ -24,14 +24,18 @@ b = FT(4pi)
 n = 128
 α = FT(0.1)
 
-domain = Domains.IntervalDomain(a, b, x3boundary = (:left, :right))
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(a),
+    Geometry.ZPoint{FT}(b),
+    boundary_tags = (:left, :right),
+)
 mesh = Meshes.IntervalMesh(domain, nelems = n)
 
 cs = Spaces.CenterFiniteDifferenceSpace(mesh)
 fs = Spaces.FaceFiniteDifferenceSpace(cs)
 
 V = Geometry.Cartesian3Vector.(ones(FT, fs))
-θ = sin.(Fields.coordinate_field(cs))
+θ = sin.(Fields.coordinate_field(cs).z)
 
 # Solve advection Equation: ∂θ/dt = -∂(vθ)
 

diff --git a/examples/column/advect_diffusion.jl b/examples/column/advect_diffusion.jl
@@ -31,7 +31,11 @@ t₁ = FT(10)
 𝓌 = FT(1)
 δ = FT(1)
 
-domain = Domains.IntervalDomain(z₀, z₁, x3boundary = (:bottom, :top))
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(z₀),
+    Geometry.ZPoint{FT}(z₁),
+    boundary_tags = (:bottom, :top),
+)
 mesh = Meshes.IntervalMesh(domain, nelems = n)
 
 cs = Spaces.CenterFiniteDifferenceSpace(mesh)
@@ -48,7 +52,7 @@ function ∇gaussian(z, t; μ = -1 // 2, ν = 1, 𝓌 = 1, δ = 1)
            exp(-(z - μ - 𝓌 * t)^2 / (4 * ν * (δ + t))) / sqrt(1 + t / δ)
 end
 
-T = gaussian.(zc, -0; μ = μ, δ = δ, ν = ν, 𝓌 = 𝓌)
+T = gaussian.(zc.z, -0; μ = μ, δ = δ, ν = ν, 𝓌 = 𝓌)
 V = Geometry.Cartesian3Vector.(ones(FT, fs))
 
 # Solve Adv-Diff Equation: ∂_t T = α ∇²T

diff --git a/examples/column/ekman.jl b/examples/column/ekman.jl
@@ -36,7 +36,11 @@ const Cd = ν / (L / nelems)
 const ug = 1.0
 const vg = 0.0
 const d = sqrt(2 * ν / f)
-domain = Domains.IntervalDomain(0.0, L; x3boundary = (:bottom, :top))
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(0.0),
+    Geometry.ZPoint{FT}(L);
+    boundary_tags = (:bottom, :top),
+)
 #mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
 mesh = Meshes.IntervalMesh(domain; nelems = nelems)
 
@@ -48,8 +52,6 @@ fspace = Spaces.FaceFiniteDifferenceSpace(cspace)
 # https://github.com/CliMA/Thermodynamics.jl/blob/main/src/TemperatureProfiles.jl#L115-L155
 # https://clima.github.io/Thermodynamics.jl/dev/TemperatureProfiles/#DecayingTemperatureProfile
 function adiabatic_temperature_profile(z; T_surf = 300.0, T_min_ref = 230.0)
-
-
     u = FT(ug)
     v = FT(vg)
     return (u = u, v = v)
@@ -58,7 +60,7 @@ end
 
 
 zc = Fields.coordinate_field(cspace)
-Yc = adiabatic_temperature_profile.(zc)
+Yc = adiabatic_temperature_profile.(zc.z)
 w = Geometry.Cartesian3Vector.(zeros(Float64, fspace))
 
 Y_init = copy(Yc)

diff --git a/examples/column/heat.jl b/examples/column/heat.jl
@@ -26,7 +26,11 @@ b = FT(1.0)
 n = 10
 α = FT(0.1)
 
-domain = Domains.IntervalDomain(a, b, x3boundary = (:bottom, :top))
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(a),
+    Geometry.ZPoint{FT}(b),
+    boundary_tags = (:bottom, :top),
+)
 mesh = Meshes.IntervalMesh(domain, nelems = n)
 
 cs = Spaces.CenterFiniteDifferenceSpace(mesh)

diff --git a/examples/column/hydrostatic.jl b/examples/column/hydrostatic.jl
@@ -29,7 +29,11 @@ const C_v = FT(R_d / (γ - 1)) # heat capacit at constant volume
 const R_m = FT(R_d) # moist R, assumed to be dry
 
 
-domain = Domains.IntervalDomain(FT(0.0), FT(30e3), x3boundary = (:bottom, :top))
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(0.0),
+    Geometry.ZPoint{FT}(30e3),
+    boundary_tags = (:bottom, :top),
+)
 #mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
 mesh = Meshes.IntervalMesh(domain; nelems = 30)
 
@@ -89,7 +93,7 @@ function discrete_hydrostatic_balance!(ρ, w, ρθ, Δz::FT, _grav::FT, Π::Func
 end
 
 zc = Fields.coordinate_field(cspace)
-Yc = decaying_temperature_profile.(zc)
+Yc = decaying_temperature_profile.(zc.z)
 w = Geometry.Cartesian3Vector.(zeros(FT, fspace))
 
 Y_init = copy(Yc)
@@ -116,7 +120,7 @@ function tendency!(dY, Y, _, t)
     @. dYc.ρθ = -(∂(w * If(Yc.ρθ)))
     @. dw = B(
         Geometry.CartesianVector(
-            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc)),
+            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc.z)),
         ),
     )
     return dY

diff --git a/examples/column/hydrostatic_ekman.jl b/examples/column/hydrostatic_ekman.jl
@@ -39,7 +39,11 @@ const ug = 1.0
 const vg = 0.0
 const uvg = Geometry.Cartesian12Vector(ug, vg)
 const d = sqrt(2 * ν / f)
-domain = Domains.IntervalDomain(0.0, L; x3boundary = (:bottom, :top))
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(0.0),
+    Geometry.ZPoint{FT}(L);
+    boundary_tags = (:bottom, :top),
+)
 #mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
 mesh = Meshes.IntervalMesh(domain; nelems = nelems)
 
@@ -86,7 +90,7 @@ end
 Φ(z) = grav * z
 
 zc = Fields.coordinate_field(cspace)
-Yc = adiabatic_temperature_profile.(zc)
+Yc = adiabatic_temperature_profile.(zc.z)
 w = Geometry.Cartesian3Vector.(zeros(FT, fspace))
 
 Y_init = copy(Yc)
@@ -154,7 +158,7 @@ function tendency!(dY, Y, _, t)
     )
     @. dw = B(
         Geometry.CartesianVector(
-            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc)),
+            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc.z)),
         ) + divf(ν * ∂c(w)) - Af(w, w),
     )
 

diff --git a/examples/column/hydrostatic_implicit.jl b/examples/column/hydrostatic_implicit.jl
@@ -37,7 +37,11 @@ const C_v = R_d / (γ - 1) # heat capacit at constant volume
 const R_m = R_d # moist R, assumed to be dry
 
 
-domain = Domains.IntervalDomain(0.0, 30e3, x3boundary = (:bottom, :top))
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(0.0),
+    Geometry.ZPoint{FT}(30e3),
+    boundary_tags = (:bottom, :top),
+)
 #mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
 mesh = Meshes.IntervalMesh(domain; nelems = 30)
 
@@ -100,9 +104,9 @@ function discrete_hydrostatic_balance!(
 end
 
 zc = Fields.coordinate_field(cspace)
-Yc = decaying_temperature_profile.(zc)
+Yc = decaying_temperature_profile.(zc.z)
 w = Geometry.Cartesian3Vector.(zeros(FT, fspace))
-zf = parent(Fields.coordinate_field(fspace))
+zf = parent(Fields.coordinate_field(fspace).z)
 Δz = zf[2:end] - zf[1:(end - 1)]
 Y_init = copy(Yc)
 w_init = copy(w)
@@ -129,7 +133,7 @@ function tendency!(dY, Y, _, t)
     @. dYc.ρθ = -(∂(w * If(Yc.ρθ)))
     @. dw = B(
         Geometry.CartesianVector(
-            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc)),
+            -(If(Yc.ρθ / Yc.ρ) * ∂f(Π(Yc.ρθ))) - ∂f(Φ(zc.z)),
         ),
     )
     return dY
@@ -188,7 +192,7 @@ function jacobian!(J, Y, p, t)
     # A_W = diagm(0=>-ones(N-1)./ρh/2, 1=>-ones(N-1)./ρh/2)[1:N-1, 1:N]
 
     # P = ([zeros(N,N)     zeros(N,N)      D_ρ;
-    #       zeros(N,N)     zeros(N,N)      D_Θ   
+    #       zeros(N,N)     zeros(N,N)      D_Θ
     #       A_W*_grav        G_W          zeros(N+1,N+1)])
 
 end

diff --git a/examples/column/wave.jl b/examples/column/wave.jl
@@ -19,25 +19,18 @@ global_logger(TerminalLogger())
 
 const FT = Float64
 
-# https://github.com/CliMA/CLIMAParameters.jl/blob/master/src/Planet/planet_parameters.jl#L5
-const MSLP = 1e5 # mean sea level pressure
-const grav = 9.8 # gravitational constant
-const R_d = 287.058 # R dry (gas constant / mol mass dry air)
-const γ = 1.4 # heat capacity ratio
-const C_p = R_d * γ / (γ - 1) # heat capacity at constant pressure
-const C_v = R_d / (γ - 1) # heat capacit at constant volume
-const R_m = R_d # moist R, assumed to be dry
-
-
-domain = Domains.IntervalDomain(0.0, 4pi, x3boundary = (:left, :right))
-#mesh = Meshes.IntervalMesh(domain, Meshes.ExponentialStretching(7.5e3); nelems = 30)
+domain = Domains.IntervalDomain(
+    Geometry.ZPoint{FT}(0.0),
+    Geometry.ZPoint{FT}(4pi),
+    boundary_tags = (:left, :right),
+)
 mesh = Meshes.IntervalMesh(domain; nelems = 30)
 
 cspace = Spaces.CenterFiniteDifferenceSpace(mesh)
 fspace = Spaces.FaceFiniteDifferenceSpace(cspace)
 
 zc = Fields.coordinate_field(cspace)
-u = sin.(zc)
+u = sin.(zc.z)
 p = Geometry.Cartesian3Vector.(zeros(Float64, fspace))
 
 using RecursiveArrayTools

diff --git a/examples/hybrid/bubble.jl b/examples/hybrid/bubble.jl
@@ -31,15 +31,15 @@ function hvspace_2D(
 )
     FT = Float64
     vertdomain = Domains.IntervalDomain(
-        FT(zlim[1]),
-        FT(zlim[2]);
-        x3boundary = (:bottom, :top),
+        Geometry.ZPoint{FT}(zlim[1]),
+        Geometry.ZPoint{FT}(zlim[2]);
+        boundary_tags = (:bottom, :top),
     )
     vertmesh = Meshes.IntervalMesh(vertdomain, nelems = velem)
     vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)
     horzdomain = Domains.RectangleDomain(
-        xlim[1]..xlim[2],
-        -0..0,
+        Geometry.XPoint{FT}(xlim[1])..Geometry.XPoint{FT}(xlim[2]),
+        Geometry.YPoint{FT}(-0)..Geometry.YPoint{FT}(0),
         x1periodic = true,
         x2boundary = (:a, :b),
     )