format

Author: MarcoArtiano

examples/elixir_euler_density_current.jl

@@ -2,52 +2,51 @@ using Trixi
 using OrdinaryDiffEq
 
 function initial_condition_density_current(x, t,
-	equations::CompressibleEulerEquations2D)
-	g = 9.81
-	c_p = 1004.0
-	c_v = 717.0
-	# center of perturbation
-	center_x = -2000.0
-	center_z = 3000.0
-	# radius of perturbation
-	radius = 1.0
-	x_r = 4000.0
-	z_r = 2000.0
-
-	# distance of current x to center of perturbation
-	r = sqrt((x[1] - center_x)^2 / x_r^2 + (x[2] - center_z)^2 / z_r^2)
-
-	# perturbation in potential temperature
-	potential_temperature_ref = 300.0
-	potential_temperature_perturbation = 0.0
-	if r <= radius
-		potential_temperature_perturbation = -15 / 2 * (1 + cospi(r))
-	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 prim2cons(SVector(rho, v1, v2, p), equations)
+                                           equations::CompressibleEulerEquations2D)
+    g = 9.81
+    c_p = 1004.0
+    c_v = 717.0
+    # center of perturbation
+    center_x = -2000.0
+    center_z = 3000.0
+    # radius of perturbation
+    radius = 1.0
+    x_r = 4000.0
+    z_r = 2000.0
+
+    # distance of current x to center of perturbation
+    r = sqrt((x[1] - center_x)^2 / x_r^2 + (x[2] - center_z)^2 / z_r^2)
+
+    # perturbation in potential temperature
+    potential_temperature_ref = 300.0
+    potential_temperature_perturbation = 0.0
+    if r <= radius
+        potential_temperature_perturbation = -15 / 2 * (1 + cospi(r))
+    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 prim2cons(SVector(rho, v1, v2, p), equations)
 end
 
 @inline function source_terms_gravity(u, x, t,
-	equations::CompressibleEulerEquations2D)
-	rho, _, rho_v2, _ = u
-	return SVector(zero(eltype(u)), zero(eltype(u)), -9.81 * rho,
-		-9.81 * rho_v2)
+                                      equations::CompressibleEulerEquations2D)
+    rho, _, rho_v2, _ = u
+    return SVector(zero(eltype(u)), zero(eltype(u)), -9.81 * rho,
+                   -9.81 * rho_v2)
 end
 
 equations = CompressibleEulerEquations2D(1.4)
@@ -59,7 +58,7 @@ surface_flux = FluxLMARS(340.0)
 volume_flux = flux_kennedy_gruber
 
 solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
-	volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ## Agnesi Profile
 L = 25_600  # Length of the box domain
@@ -75,15 +74,15 @@ f4(s) = SVector((s + 1 - 1) * L / 2, H)
 
 trees_per_dimension = (32, 32)
 mesh = P4estMesh(trees_per_dimension, polydeg = 3,
-	faces = (f1, f2, f3, f4),
-	initial_refinement_level = 1, periodicity = (true, false))
+                 faces = (f1, f2, f3, f4),
+                 initial_refinement_level = 1, periodicity = (true, false))
 
-boundary_conditions = Dict(
-	:y_neg => boundary_condition_slip_wall,
-	:y_pos => boundary_condition_slip_wall)
+boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall,
+                           :y_pos => boundary_condition_slip_wall)
 
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_current, solver, source_terms = source_terms_gravity,
-	boundary_conditions = boundary_conditions)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_current,
+                                    solver, source_terms = source_terms_gravity,
+                                    boundary_conditions = boundary_conditions)
 
 tspan = (0.0, 900.0)
 ode = semidiscretize(semi, tspan)
@@ -96,11 +95,11 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
 callbacks = CallbackSet(summary_callback,
-	analysis_callback,
-	alive_callback, visualization_callback)
+                        analysis_callback,
+                        alive_callback, visualization_callback)
 
 sol = solve(ode,
-SSPRK43(thread = OrdinaryDiffEq.True()),
-	maxiters = 1.0e7,
-	dt = 1e-1,
-	save_everystep = false, callback = callbacks)
+            SSPRK43(thread = OrdinaryDiffEq.True()),
+            maxiters = 1.0e7,
+            dt = 1e-1,
+            save_everystep = false, callback = callbacks)

examples/elixir_euler_linear_hydrostatic.jl

@@ -2,97 +2,99 @@ using Trixi
 using OrdinaryDiffEq
 
 struct HydrostaticSetup
-	# Physical constants
-	g::Float64       # gravity of earth
-	c_p::Float64     # heat capacity for constant pressure (dry air)
-	c_v::Float64     # heat capacity for constant volume (dry air)
-	gamma::Float64   # heat capacity ratio (dry air)
-	p_0::Float64     # atmospheric pressure
-	T_0::Float64     # 
-	u0::Float64      #
-	Nf::Float64      # 
-	z_B::Float64     #
-	z_T::Float64     #
-	function HydrostaticSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v, p_0 = 100_000.0, T_0 = 250.0, u0 = 20.0, z_B = 15000.0, z_T = 30000.0)
-		Nf = g / sqrt(c_p * T_0)
-		new(g, c_p, c_v, gamma, p_0, T_0, u0, Nf, z_B, z_T)
-	end
+    # Physical constants
+    g::Float64       # gravity of earth
+    c_p::Float64     # heat capacity for constant pressure (dry air)
+    c_v::Float64     # heat capacity for constant volume (dry air)
+    gamma::Float64   # heat capacity ratio (dry air)
+    p_0::Float64     # atmospheric pressure
+    T_0::Float64     # 
+    u0::Float64      #
+    Nf::Float64      # 
+    z_B::Float64     #
+    z_T::Float64     #
+    function HydrostaticSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v,
+                              p_0 = 100_000.0, T_0 = 250.0, u0 = 20.0, z_B = 15000.0,
+                              z_T = 30000.0)
+        Nf = g / sqrt(c_p * T_0)
+        new(g, c_p, c_v, gamma, p_0, T_0, u0, Nf, z_B, z_T)
+    end
 end
 
 function (setup::HydrostaticSetup)(u, x, t, equations::CompressibleEulerEquations2D)
-	@unpack g, c_p, c_v, gamma, p_0, T_0, z_B, z_T, Nf, u0 = setup
+    @unpack g, c_p, c_v, gamma, p_0, T_0, z_B, z_T, Nf, u0 = setup
 
-	rho, rho_v1, rho_v2, rho_e = u
+    rho, rho_v1, rho_v2, rho_e = u
 
-	R = c_p - c_v
+    R = c_p - c_v
 
-	v1 = rho_v1 / rho
-	v2 = rho_v2 / rho
-	p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
-	rho_theta = (p / p_0)^(c_v / c_p) * p_0 / R
-	theta = rho_theta / rho
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1 * v1 + rho_v2 * v2))
+    rho_theta = (p / p_0)^(c_v / c_p) * p_0 / R
+    theta = rho_theta / rho
 
-	alfa = 0.1
+    alfa = 0.1
 
-	xr_B = 80000.0
-	xr_T = 120000.0
+    xr_B = 80000.0
+    xr_T = 120000.0
 
-	if x[2] <= z_B
-		S_v = 0.0
-	else
-		S_v = -alfa * sinpi(0.5 * (x[2] - z_B) / (z_T - z_B))^2
-	end
-	if x[1] < xr_B
-		S_h1 = 0.0
-	else
-		S_h1 = -alfa * sinpi(0.5 * (x[1] - xr_B) / (xr_T - xr_B))^2
-	end
+    if x[2] <= z_B
+        S_v = 0.0
+    else
+        S_v = -alfa * sinpi(0.5 * (x[2] - z_B) / (z_T - z_B))^2
+    end
+    if x[1] < xr_B
+        S_h1 = 0.0
+    else
+        S_h1 = -alfa * sinpi(0.5 * (x[1] - xr_B) / (xr_T - xr_B))^2
+    end
 
-	if x[1] > -xr_B
-		S_h2 = 0.0
-	else
-		S_h2 = -alfa * sinpi(0.5 * (x[1] + xr_B) / (-xr_T + xr_B))^2
-	end
+    if x[1] > -xr_B
+        S_h2 = 0.0
+    else
+        S_h2 = -alfa * sinpi(0.5 * (x[1] + xr_B) / (-xr_T + xr_B))^2
+    end
 
-	exner = exp(-Nf^2 / g * x[2])
+    exner = exp(-Nf^2 / g * x[2])
 
-	theta_0 = T_0 / exner
+    theta_0 = T_0 / exner
 
-	K = p_0 * (R / p_0)^gamma
-	du2 = rho * (v1 - u0) * (S_v + S_h1 + S_h2)
-	du3 = rho_v2 * (S_v + S_h1 + S_h2)
-	du4 = rho * (theta - theta_0) * (S_v + S_h1 + S_h2) * K * gamma / (gamma - 1.0) * (rho_theta)^(gamma - 1.0) + du2 * v1 + du3 * v2
-
-	return SVector(zero(eltype(u)), du2, du3 - g * rho, du4 - g * rho_v2)
+    K = p_0 * (R / p_0)^gamma
+    du2 = rho * (v1 - u0) * (S_v + S_h1 + S_h2)
+    du3 = rho_v2 * (S_v + S_h1 + S_h2)
+    du4 = rho * (theta - theta_0) * (S_v + S_h1 + S_h2) * K * gamma / (gamma - 1.0) *
+          (rho_theta)^(gamma - 1.0) + du2 * v1 + du3 * v2
 
+    return SVector(zero(eltype(u)), du2, du3 - g * rho, du4 - g * rho_v2)
 end
 
 function (setup::HydrostaticSetup)(x, t, equations::CompressibleEulerEquations2D)
-	@unpack g, c_p, c_v, p_0, T_0, u0, Nf = setup
+    @unpack g, c_p, c_v, p_0, T_0, u0, Nf = setup
 
-	# Exner pressure, solves hydrostatic equation for x[2]
-	exner = exp(-Nf^2 / g * 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)
+    # Exner pressure, solves hydrostatic equation for x[2]
+    exner = exp(-Nf^2 / g * 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)
 
-	# density
-	rho = p / (R * T_0)
-	v1 = u0
-	v2 = 0.0
+    # density
+    rho = p / (R * T_0)
+    v1 = u0
+    v2 = 0.0
 
-	return prim2cons(SVector(rho, v1, v2, p), equations)
+    return prim2cons(SVector(rho, v1, v2, p), equations)
 end
 
 linear_hydrostatic_setup = HydrostaticSetup()
 
 equations = CompressibleEulerEquations2D(linear_hydrostatic_setup.gamma)
 
 boundary_conditions = (x_neg = boundary_condition_periodic,
-	x_pos = boundary_condition_periodic,
-	y_neg = boundary_condition_slip_wall,
-	y_pos = boundary_condition_slip_wall)
+                       x_pos = boundary_condition_periodic,
+                       y_neg = boundary_condition_slip_wall,
+                       y_pos = boundary_condition_slip_wall)
 
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
@@ -118,8 +120,9 @@ f4(s) = SVector((s + 1 - 1) * L / 2, H)
 cells_per_dimension = (100, 50)
 mesh = StructuredMesh(cells_per_dimension, (f1, f2, f3, f4), periodicity = (true, false))
 
-semi = SemidiscretizationHyperbolic(mesh, equations, linear_hydrostatic_setup, solver, source_terms = linear_hydrostatic_setup,
-	boundary_conditions = boundary_conditions)
+semi = SemidiscretizationHyperbolic(mesh, equations, linear_hydrostatic_setup, solver,
+                                    source_terms = linear_hydrostatic_setup,
+                                    boundary_conditions = boundary_conditions)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
@@ -136,15 +139,15 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
 callbacks = CallbackSet(summary_callback,
-	analysis_callback,
-	alive_callback)
+                        analysis_callback,
+                        alive_callback)
 
 ###############################################################################
 # run the simulation
 sol = solve(ode,
-	SSPRK43(thread = OrdinaryDiffEq.True()),
-	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)
+            SSPRK43(thread = OrdinaryDiffEq.True()),
+            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()