Potential Temperature Euler Eq.

examples/elixir_potential_euler_dry_bubble.jl

@@ -0,0 +1,114 @@
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+using TrixiAtmo: flux_LMARS, source_terms_bubble, flux_theta
+using Plots
+
+function initial_condition_warm_bubble(x, t, equations::CompressiblePotentialEulerEquations2D)
+	g = equations.g
+	c_p = equations.c_p
+	c_v = equations.c_v
+	# center of perturbation
+	center_x = 10000.0
+	center_z = 2000.0
+	# radius of perturbation
+	radius = 2000.0
+	# distance of current x to center of perturbation
+	r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2)
+
+	# perturbation in potential temperature
+	potential_temperature_ref = 300.0
+	potential_temperature_perturbation = 0.0
+	if r <= radius
+		potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2
+	end
+	potential_temperature = potential_temperature_ref + potential_temperature_perturbation
+
+	# Exner pressure, solves hydrostatic equation for x[2]
+	exner = 1 - g / (c_p * potential_temperature_ref) * x[2]
+
+	# pressure
+	p_0 = 100_000.0  # reference pressure
+	R = c_p - c_v    # gas constant (dry air)
+	p = p_0 * exner^(c_p / R)
+	T = potential_temperature * exner
+
+	# density
+	rho = p / (R * T)
+	v1 = 20.0
+	v2 = 0.0
+
+	return SVector(rho, rho * v1, rho * v2, rho * potential_temperature)
+end
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+
+equations = CompressiblePotentialEulerEquations2D()
+
+boundary_conditions = (x_neg = boundary_condition_periodic,
+	x_pos = boundary_condition_periodic,
+	y_neg = boundary_condition_slip_wall,
+	y_pos = boundary_condition_slip_wall)
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+
+surface_flux = flux_LMARS
+
+volume_flux = flux_theta
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (20_000.0, 10_000.0)
+
+cells_per_dimension = (32, 16)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+	periodicity = (true, false))
+initial_condition = initial_condition_warm_bubble
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+	source_terms = source_terms_bubble,
+	boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1000.0)  # 1000 seconds final time
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+	extra_analysis_errors = (:entropy_conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = analysis_interval,
+	save_initial_solution = true,
+	save_final_solution = true,
+	output_directory = "out",
+	solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback,
+	analysis_callback,
+	alive_callback,
+	save_solution,
+	stepsize_callback)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+	maxiters = 1.0e7,
+	dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+	save_everystep = false, callback = callbacks);
+
+summary_callback()

src/TrixiAtmo.jl

@@ -24,8 +24,9 @@ include("solvers/solvers.jl")
 include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
 
 export CompressibleMoistEulerEquations2D
+export CompressiblePotentialEulerEquations2D
 
-export flux_chandrashekar, flux_LMARS
+export flux_chandrashekar, flux_LMARS, flux_theta
 
 export examples_dir