Add Lucas' thesis

examples/elixir_moist_euler_EC_bubble.jl

@@ -0,0 +1,275 @@
+
+using OrdinaryDiffEq
+using Trixi, TrixiAtmo
+using TrixiAtmo: cons2aeqpot, source_terms_geopotential
+using NLsolve: nlsolve
+
+###############################################################################
+# semidiscretization of the compressible moist Euler equations
+
+equations = CompressibleMoistEulerEquations2D()
+
+function moist_state(y, dz, y0, r_t0, theta_e0, equations::CompressibleMoistEulerEquations2D)
+  @unpack p_0, g, c_pd, c_pv, c_vd, c_vv, R_d, R_v, c_pl, L_00 = equations
+  (p, rho, T, r_t, r_v, rho_qv, theta_e) = y
+  p0 = y0[1]
+
+  F = zeros(7,1)
+  rho_d = rho / (1 + r_t)
+  p_d = R_d * rho_d * T
+  T_C = T - 273.15
+  p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C))
+  L = L_00 - (c_pl - c_pv) * T
+
+  F[1] = (p - p0) / dz + g * rho 
+  F[2] = p - (R_d * rho_d + R_v * rho_qv) * T
+  # H = 1 is assumed
+  F[3] = (theta_e - T * (p_d / p_0)^(-R_d / (c_pd + c_pl * r_t)) *
+         exp(L * r_v / ((c_pd + c_pl * r_t) * T)))
+  F[4] = r_t - r_t0
+  F[5] = rho_qv - rho_d * r_v
+  F[6] = theta_e - theta_e0
+  a = p_vs / (R_v * T) - rho_qv
+  b = rho - rho_qv - rho_d
+  # H=1 => phi=0
+  F[7] = a+b-sqrt(a*a+b*b)
+
+  return F
+end
+
+function generate_function_of_y(dz, y0, r_t0, theta_e0, equations::CompressibleMoistEulerEquations2D)
+  function function_of_y(y)
+    return moist_state(y, dz, y0, r_t0, theta_e0, equations)
+  end
+end
+#Create Initial atmosphere by generating a layer data set
+struct AtmossphereLayers{RealT<:Real}
+  equations::CompressibleMoistEulerEquations2D
+  # structure:  1--> i-layer (z = total_hight/precision *(i-1)),  2--> rho, rho_theta, rho_qv, rho_ql
+  LayerData::Matrix{RealT} # Contains the layer data for each height
+  total_hight::RealT # Total height of the atmosphere
+  preciseness::Int # Space between each layer data (dz)
+  layers::Int # Amount of layers (total height / dz)
+  ground_state::NTuple{2, RealT} # (rho_0, p_tilde_0) to define the initial values at height z=0
+  equivalentpotential_temperature::RealT  # Value for theta_e since we have a constant temperature theta_e0=theta_e.
+  mixing_ratios::NTuple{2, RealT} # Constant mixing ratio. Also defines initial guess for rho_qv0 = r_v0 * rho_0.
+end
+
+
+function AtmossphereLayers(equations ; total_hight=10010.0, preciseness=10, ground_state=(1.4, 100000.0), equivalentpotential_temperature=320, mixing_ratios=(0.02, 0.02), RealT=Float64)
+  @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl = equations
+  rho0, p0 = ground_state
+  r_t0, r_v0 = mixing_ratios
+  theta_e0 = equivalentpotential_temperature
+
+  rho_qv0 = rho0 * r_v0
+  T0 = theta_e0
+  y0 = [p0, rho0, T0, r_t0, r_v0, rho_qv0, theta_e0]
+
+  n = convert(Int, total_hight / preciseness)
+  dz = 0.01
+  LayerData = zeros(RealT, n+1, 4)
+
+  F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations)
+  sol = nlsolve(F, y0)
+  p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero
+
+  rho_d = rho / (1 + r_t)
+  rho_ql = rho - rho_d - rho_qv
+  kappa_M=(R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql)
+  rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) *r_v) / (1 + r_t) 
+
+  LayerData[1, :] = [rho, rho_theta, rho_qv, rho_ql]
+   for i in (1:n)
+    y0 = deepcopy(sol.zero)
+    dz = preciseness
+    F = generate_function_of_y(dz, y0, r_t0, theta_e0, equations)
+    sol = nlsolve(F, y0)
+    p, rho, T, r_t, r_v, rho_qv, theta_e = sol.zero
+
+    rho_d = rho / (1 + r_t)
+    rho_ql = rho - rho_d - rho_qv
+    kappa_M=(R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql)
+    rho_theta = rho * (p0 / p)^kappa_M * T * (1 + (R_v / R_d) *r_v) / (1 + r_t) 
+
+    LayerData[i+1, :] = [rho, rho_theta, rho_qv, rho_ql]
+   end
+
+  return AtmossphereLayers{RealT}(equations, LayerData, total_hight, dz, n, ground_state, theta_e0, mixing_ratios)
+end
+
+# Generate background state from the Layer data set by linearely interpolating the layers
+function initial_condition_moist_bubble(x, t, equations::CompressibleMoistEulerEquations2D, AtmosphereLayers::AtmossphereLayers)
+  @unpack LayerData, preciseness, total_hight = AtmosphereLayers
+  dz = preciseness
+  z = x[2] 
+  if (z > total_hight && !(isapprox(z, total_hight)))
+    error("The atmossphere does not match the simulation domain")
+  end
+  n = convert(Int, floor(z/dz)) + 1
+  z_l = (n-1) * dz
+  (rho_l, rho_theta_l, rho_qv_l, rho_ql_l) = LayerData[n, :]
+  z_r = n * dz
+  if (z_l == total_hight)
+    z_r = z_l + dz 
+    n = n-1
+  end
+  (rho_r, rho_theta_r, rho_qv_r, rho_ql_r) = LayerData[n+1, :]
+  rho = (rho_r * (z - z_l) + rho_l * (z_r - z)) / dz
+  rho_theta = rho * (rho_theta_r / rho_r * (z - z_l) + rho_theta_l / rho_l * (z_r - z)) / dz
+  rho_qv = rho * (rho_qv_r / rho_r * (z - z_l) + rho_qv_l / rho_l * (z_r - z)) / dz
+  rho_ql = rho * (rho_ql_r / rho_r * (z - z_l) + rho_ql_l / rho_l * (z_r - z)) / dz
+
+  rho, rho_e, rho_qv, rho_ql = PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D)
+
+  v1 = 60.0
+  v2 = 60.0
+  rho_v1 = rho * v1
+  rho_v2 = rho * v2
+  rho_E = rho_e + 1/2 * rho *(v1^2 + v2^2)
+
+
+  return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql)
+end
+
+# Add perturbation to the profile
+function PerturbMoistProfile(x, rho, rho_theta, rho_qv, rho_ql, equations::CompressibleMoistEulerEquations2D)
+  @unpack kappa, p_0, c_pd, c_vd, c_pv, c_vv, R_d, R_v, c_pl, L_00 = equations
+  xc = 2000
+  zc = 2000
+  rc = 2000
+  # Peak perturbation at center
+  Δθ = 10
+
+  r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2)
+  rho_d = rho - rho_qv - rho_ql
+  kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql)
+  p_loc = p_0 *(R_d * rho_theta / p_0)^(1/(1-kappa_M))
+  T_loc = p_loc / (R_d * rho_d + R_v * rho_qv)
+  rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv
+
+
+  # Assume pressure stays constant
+  if (r < rc && Δθ > 0) 
+    θ_dens = rho_theta / rho * (p_loc / p_0)^(kappa_M - kappa)
+    θ_dens_new = θ_dens * (1 + Δθ * cospi(0.5*r/rc)^2 / 300)
+    rt =(rho_qv + rho_ql) / rho_d 
+    rv = rho_qv / rho_d
+    θ_loc = θ_dens_new * (1 + rt)/(1 + (R_v / R_d) * rv)
+    if rt > 0 
+      while true 
+        T_loc = θ_loc * (p_loc / p_0)^kappa
+        T_C = T_loc - 273.15
+        # SaturVapor
+        pvs = 611.2 * exp(17.62 * T_C / (243.12 + T_C))
+        rho_d_new = (p_loc - pvs) / (R_d * T_loc)
+        rvs = pvs / (R_v * rho_d_new * T_loc)
+        θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs)
+        if abs(θ_new-θ_loc) <= θ_loc * 1.0e-12
+          break
+        else
+          θ_loc=θ_new
+        end
+      end
+    else
+      rvs = 0
+      T_loc = θ_loc * (p_loc / p_0)^kappa
+      rho_d_new = p_loc / (R_d * T_loc)
+      θ_new = θ_dens_new * (1 + rt) / (1 + (R_v / R_d) * rvs)
+    end
+    rho_qv = rvs * rho_d_new
+    rho_ql = (rt - rvs) * rho_d_new
+    rho = rho_d_new * (1 + rt)
+    rho_d = rho - rho_qv - rho_ql
+    kappa_M = (R_d * rho_d + R_v * rho_qv) / (c_pd * rho_d + c_pv * rho_qv + c_pl * rho_ql)
+    rho_theta = rho * θ_dens_new * (p_loc / p_0)^(kappa - kappa_M)
+    rho_e = (c_vd * rho_d + c_vv * rho_qv + c_pl * rho_ql) * T_loc + L_00 * rho_qv
+  end
+
+  return SVector(rho, rho_e, rho_qv, rho_ql)
+end
+
+AtmossphereData = AtmossphereLayers(equations)
+
+function initial_condition_moist(x, t, equations)
+  return initial_condition_moist_bubble(x, t, equations, AtmossphereData)
+end
+
+initial_condition = initial_condition_moist
+
+boundary_condition = (x_neg=boundary_condition_periodic,
+                      x_pos=boundary_condition_periodic,
+                      y_neg=boundary_condition_periodic,
+                      y_pos=boundary_condition_periodic)
+
+source_term = source_terms_geopotential
+
+###############################################################################
+# Get the DG approximation space
+
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+
+
+surface_flux = flux_chandrashekar
+volume_flux = flux_chandrashekar
+
+
+volume_integral=VolumeIntegralFluxDifferencing(volume_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (4000.0, 4000.0)
+
+cells_per_dimension = (16, 16)
+
+# Create curved mesh with 16 x 16 elements
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_condition,
+                                    source_terms=source_term)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 30.0)
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+solution_variables = cons2aeqpot
+
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, extra_analysis_errors=(:entropy_conservation_error, ))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=1000,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=solution_variables)
+
+stepsize_callback = StepsizeCallback(cfl=0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+summary_callback() # print the timer summary