From cd2ed5eb34257a13db6d049da43a597e2ea1ef8e Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Fri, 19 Jul 2024 17:19:39 +0200 Subject: [PATCH 1/5] Update README.md Specifying ?branch=main in the CI badge. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 986afe6..41c165f 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ --> [![Docs-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://trixi-framework.github.io/TrixiAtmo.jl/dev) [![Slack](https://img.shields.io/badge/chat-slack-e01e5a)](https://join.slack.com/t/trixi-framework/shared_invite/zt-sgkc6ppw-6OXJqZAD5SPjBYqLd8MU~g) -[![Build Status](https://github.com/trixi-framework/TrixiAtmo.jl/workflows/CI/badge.svg)](https://github.com/trixi-framework/TrixiAtmo.jl/actions?query=workflow%3ACI) +[![Build Status](https://github.com/trixi-framework/TrixiAtmo.jl/workflows/CI/badge.svg?branch=main)](https://github.com/trixi-framework/TrixiAtmo.jl/actions?query=workflow%3ACI) [![Codecov](https://codecov.io/gh/trixi-framework/TrixiAtmo.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/trixi-framework/TrixiAtmo.jl) [![Coveralls](https://coveralls.io/repos/github/trixi-framework/TrixiAtmo.jl/badge.svg?branch=main)](https://coveralls.io/github/trixi-framework/TrixiAtmo.jl?branch=main) [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT) From c9404d7375d8838f337eb1dcd4e48d27268fe9af Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Fri, 19 Jul 2024 17:32:00 +0200 Subject: [PATCH 2/5] Update README.md Ok, apparently specifying ?branch=main in the CI badge doesn't work. Now I think failing PRs will cause the badge to say "failing." --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 41c165f..986afe6 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ --> [![Docs-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://trixi-framework.github.io/TrixiAtmo.jl/dev) [![Slack](https://img.shields.io/badge/chat-slack-e01e5a)](https://join.slack.com/t/trixi-framework/shared_invite/zt-sgkc6ppw-6OXJqZAD5SPjBYqLd8MU~g) -[![Build Status](https://github.com/trixi-framework/TrixiAtmo.jl/workflows/CI/badge.svg?branch=main)](https://github.com/trixi-framework/TrixiAtmo.jl/actions?query=workflow%3ACI) +[![Build Status](https://github.com/trixi-framework/TrixiAtmo.jl/workflows/CI/badge.svg)](https://github.com/trixi-framework/TrixiAtmo.jl/actions?query=workflow%3ACI) [![Codecov](https://codecov.io/gh/trixi-framework/TrixiAtmo.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/trixi-framework/TrixiAtmo.jl) [![Coveralls](https://coveralls.io/repos/github/trixi-framework/TrixiAtmo.jl/badge.svg?branch=main)](https://coveralls.io/github/trixi-framework/TrixiAtmo.jl?branch=main) [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT) From a9f91cacfb98677c74de496802745523b51c0ba1 Mon Sep 17 00:00:00 2001 From: Benedict <135045760+benegee@users.noreply.github.com> Date: Tue, 23 Jul 2024 11:36:32 +0200 Subject: [PATCH 3/5] adapt URL scheme according to github docs (#22) --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 986afe6..d5a0287 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ --> [![Docs-dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://trixi-framework.github.io/TrixiAtmo.jl/dev) [![Slack](https://img.shields.io/badge/chat-slack-e01e5a)](https://join.slack.com/t/trixi-framework/shared_invite/zt-sgkc6ppw-6OXJqZAD5SPjBYqLd8MU~g) -[![Build Status](https://github.com/trixi-framework/TrixiAtmo.jl/workflows/CI/badge.svg)](https://github.com/trixi-framework/TrixiAtmo.jl/actions?query=workflow%3ACI) +[![Build Status](https://github.com/trixi-framework/TrixiAtmo.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/trixi-framework/TrixiAtmo.jl/actions?query=workflow%3ACI) [![Codecov](https://codecov.io/gh/trixi-framework/TrixiAtmo.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/trixi-framework/TrixiAtmo.jl) [![Coveralls](https://coveralls.io/repos/github/trixi-framework/TrixiAtmo.jl/badge.svg?branch=main)](https://coveralls.io/github/trixi-framework/TrixiAtmo.jl?branch=main) [![License: MIT](https://img.shields.io/badge/License-MIT-success.svg)](https://opensource.org/licenses/MIT) From 25592b60d5699a8e90319ddf4817551d12fe10ab Mon Sep 17 00:00:00 2001 From: Benedict <135045760+benegee@users.noreply.github.com> Date: Tue, 23 Jul 2024 22:31:41 +0200 Subject: [PATCH 4/5] Do coverage testing only on ubuntu (#23) --- .github/workflows/ci.yml | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index b9eb5d3..e0b76fe 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -48,21 +48,31 @@ jobs: show-versioninfo: true - uses: julia-actions/cache@v2 - - name: Run tests + - name: Run tests without coverage + if: ${{ matrix.os != 'ubuntu-latest' }} + uses: julia-actions/julia-runtest@v1 + with: + coverage: false + + - name: Run tests with coverage + if: ${{ matrix.os == 'ubuntu-latest' }} uses: julia-actions/julia-runtest@v1 with: coverage: true - name: Process Julia coverage data + if: ${{ matrix.os == 'ubuntu-latest' }} uses: julia-actions/julia-processcoverage@v1 - name: Upload coverage data (Codecov) + if: ${{ matrix.os == 'ubuntu-latest' }} uses: codecov/codecov-action@v4 with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} - name: Upload coverage data (Coveralls) + if: ${{ matrix.os == 'ubuntu-latest' }} uses: coverallsapp/github-action@v2 with: github-token: ${{ secrets.GITHUB_TOKEN }} From 32fa7c76e28408c65525682cce6d4d1d6b8a445c Mon Sep 17 00:00:00 2001 From: Benedict <135045760+benegee@users.noreply.github.com> Date: Wed, 24 Jul 2024 00:55:09 +0200 Subject: [PATCH 5/5] Add equations by Lucas (#21) add equations by Lucas and consistency check --------- Co-authored-by: Tristan Montoya --- Project.toml | 2 + src/TrixiAtmo.jl | 15 +- src/auxiliary/auxiliary.jl | 14 + .../compressible_moist_euler_2d_lucas.jl | 1032 +++++++++++++++++ src/equations/equations.jl | 4 + test/Project.toml | 6 +- test/runtests.jl | 27 +- test/test_trixi_consistency.jl | 64 + test/test_trixiatmo.jl | 51 + 9 files changed, 1190 insertions(+), 25 deletions(-) create mode 100644 src/auxiliary/auxiliary.jl create mode 100644 src/equations/compressible_moist_euler_2d_lucas.jl create mode 100644 test/test_trixi_consistency.jl create mode 100644 test/test_trixiatmo.jl diff --git a/Project.toml b/Project.toml index d468d75..1b2a381 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.0-DEV" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" @@ -13,6 +14,7 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] LinearAlgebra = "1" MuladdMacro = "0.2.2" +Reexport = "1.0" Static = "0.8" StaticArrays = "1" Trixi = "0.7, 0.8" diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 84eac18..61b3702 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -9,17 +9,20 @@ See also: [trixi-framework/TrixiAtmo.jl](https://github.com/trixi-framework/Trix module TrixiAtmo using Trixi -# Import additional symbols that are not exported by Trixi.jl -# using Trixi: using MuladdMacro: @muladd using StaticArrays: SVector using Static: True, False using LinearAlgebra: norm +using Reexport: @reexport +@reexport using StaticArrays: SVector -foo() = true -bar() = false -baz() = Trixi.examples_dir() - +include("auxiliary/auxiliary.jl") include("equations/equations.jl") +export CompressibleMoistEulerEquations2D + +export flux_chandrashekar, flux_LMARS + +export examples_dir + end # module TrixiAtmo diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl new file mode 100644 index 0000000..316a687 --- /dev/null +++ b/src/auxiliary/auxiliary.jl @@ -0,0 +1,14 @@ +""" + examples_dir() + +Return the directory where the example files provided with TrixiAtmo.jl are located. +If TrixiAtmo.jl is installed as a regular package (with `]add Trixi`), these files are +read-only and should *not* be modified. +To find out which files are available, use, e.g., `readdir`: + +# Examples +```@example +readdir(examples_dir()) +``` +""" +examples_dir() = pkgdir(TrixiAtmo, "examples") diff --git a/src/equations/compressible_moist_euler_2d_lucas.jl b/src/equations/compressible_moist_euler_2d_lucas.jl new file mode 100644 index 0000000..f9941a3 --- /dev/null +++ b/src/equations/compressible_moist_euler_2d_lucas.jl @@ -0,0 +1,1032 @@ +# Implemented by Lucas Gemein +# https://github.com/NichtLucas/Trixi.jl/tree/thesis_gemein_2022 + +using Trixi +using Trixi: ln_mean, inv_ln_mean +import Trixi: varnames, flux_chandrashekar, boundary_condition_slip_wall, + cons2prim, cons2entropy, max_abs_speed_naive, max_abs_speeds, + entropy, energy_total, flux + +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + struct CompressibleMoistEulerEquations2D{RealT <: Real} <: + AbstractCompressibleMoistEulerEquations{2, 6} + p_0::RealT # constant reference pressure 1000 hPa(100000 Pa) + c_pd::RealT # dry air constant + c_vd::RealT # dry air constant + R_d::RealT # dry air gas constant + c_pv::RealT # moist air constant + c_vv::RealT # moist air constant + R_v::RealT # moist air gas constant + c_pl::RealT # liqid water constant + g::RealT # gravitation constant + kappa::RealT # ratio of the gas constant R_d + gamma::RealT # = inv(kappa- 1); can be used to write slow divisions as fast multiplications + L_00::RealT # latent heat of evaporation at 0 K + a::RealT + end + + function CompressibleMoistEulerEquations2D(; g = 9.81, RealT = Float64) + p_0 = 100000.0 + c_pd = 1004.0 + c_vd = 717.0 + R_d = c_pd - c_vd + c_pv = 1885.0 + c_vv = 1424.0 + R_v = c_pv - c_vv + c_pl = 4186.0 + gamma = c_pd / c_vd # = 1/(1 - kappa) + kappa = 1 - inv(gamma) + L_00 = 3147620.0 + a = 360.0 + return CompressibleMoistEulerEquations2D{RealT}(p_0, c_pd, c_vd, R_d, c_pv, c_vv, + R_v, c_pl, g, kappa, gamma, L_00, a) + end + + # Calculate 1D flux for a single point. + @inline function flux(u, orientation::Integer, + equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + qv = rho_qv / rho + ql = rho_ql / rho + p, T = get_current_condition(u, equations) + if orientation == 1 + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + f4 = (rho_E + p) * v1 + f5 = rho_v1 * qv + f6 = rho_v1 * ql + else + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + f4 = (rho_E + p) * v2 + f5 = rho_v2 * qv + f6 = rho_v2 * ql + end + return SVector(f1, f2, f3, f4, f5, f6) + end + + # Calculate 1D flux for a single point in the normal direction. + # Note, this directional vector is not normalized. + @inline function flux(u, normal_direction::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_e, rho_qv, rho_ql = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + qv = rho_qv / rho + ql = rho_ql / rho + p, T = get_current_condition(u, equations) + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + rho_v_normal = rho * v_normal + f1 = rho_v_normal + f2 = (rho_v_normal) * v1 + p * normal_direction[1] + f3 = (rho_v_normal) * v2 + p * normal_direction[2] + f4 = (rho_e + p) * v_normal + f5 = rho_v_normal * qv + f6 = (rho_v_normal) * ql + return SVector(f1, f2, f3, f4, f5, f6) + end + + # Slip-wall boundary condition + # Determine the boundary numerical surface flux for a slip wall condition. + # Imposes a zero normal velocity at the wall. + @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_pv, c_pl, c_vd, c_vv = equations + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal = normal_direction / norm_ + + # rotate the internal solution state + u_local = rotate_to_x(u_inner, normal, equations) + + # compute the primitive variables + rho_local, v_normal, v_tangent, p_local, qv_local, ql_local = cons2prim(u_local, + equations) + qd_local = 1 - qv_local - ql_local + gamma = (qd_local * c_pd + qv_local * c_pv + ql_local * c_pl) * + inv(qd_local * c_vd + qv_local * c_vv + ql_local * c_pl) + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0.0 + sound_speed = sqrt(gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1.0 + 0.5 * (gamma - 1) * v_normal / sound_speed)^(2.0 * gamma * + inv(gamma - 1)) + else # v_normal > 0.0 + A = 2.0 / ((gamma + 1) * rho_local) + B = p_local * (gamma - 1) / (gamma + 1) + p_star = p_local + + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4.0 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return SVector(zero(eltype(u_inner)), + p_star * normal[1], + p_star * normal[2], + zero(eltype(u_inner)), + zero(eltype(u_inner)), + zero(eltype(u_inner))) * norm_ + end + + # Fix sign for structured mesh. + @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::CompressibleMoistEulerEquations2D) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x and -y sides due to the orientation convention used by StructuredMesh + if isodd(direction) + boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, + equations) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, + equations) + end + + return boundary_flux + end + + # Rotate momentum flux. The same as in compressible Euler. + @inline function rotate_to_x(u, normal_vector::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + # cos and sin of the angle between the x-axis and the normalized normal_vector are + # the normalized vector's x and y coordinates respectively (see unit circle). + c = normal_vector[1] + s = normal_vector[2] + + # Apply the 2D rotation matrix with normal and tangent directions of the form + # [ 1 0 0 0; + # 0 n_1 n_2 0; + # 0 t_1 t_2 0; + # 0 0 0 1 ] + # where t_1 = -n_2 and t_2 = n_1 + + return SVector(u[1], + c * u[2] + s * u[3], + -s * u[2] + c * u[3], + u[4], u[5], u[6]) + end + + # Recreates the convergence test initial condition from compressible euler 2D. + function initial_condition_convergence_test_dry(x, t, + equations::CompressibleMoistEulerEquations2D) + @unpack R_d, R_v, c_vd, c_vv, c_pl, L_00 = equations + + c = 2 + A = 0.1 + L = 2 + f = 1 / L + ω = 2 * pi * f + ini = c + A * sin(ω * (x[1] + x[2] - t)) + + qv = 0 + ql = 0 + + mu = ((1 - qv - ql) * c_vd + qv * c_vv + ql * c_pl) + + T = (ini - 1) / c_vd + E = (mu * T + qv * L_00 + 1) + + rho = ini + rho_v1 = ini + rho_v2 = ini + rho_e = E * ini + rho_qv = qv * ini + rho_ql = ql * ini + + return SVector(rho, rho_v1, rho_v2, rho_e, rho_qv, rho_ql) + end + + # Recreates the convergence test source term from compressible euler 2D. + @inline function source_terms_convergence_test_dry(u, x, t, + equations::CompressibleMoistEulerEquations2D) + # Same settings as in `initial_condition` + @unpack R_d, R_v, c_vd, c_vv, c_pl, L_00 = equations + c = 2 + A = 0.1 + L = 2 + f = 1 / L + ω = 2 * pi * f + + x1, x2 = x + si, co = sincos(ω * (x1 + x2 - t)) + rho = c + A * si + rho_x = ω * A * co + + qv = 0 + ql = 0 + mu = ((1 - qv - ql) * c_vd + qv * c_vv + ql * c_pl) + xi = ((1 - qv - ql) * R_d + qv * R_v) + + T = (rho - 1) / c_vd + dT = rho_x / c_vd + E = (mu * T + qv * L_00 + 1) + dE = E * rho_x + rho * mu * dT + dp = xi * (T * rho_x + rho * dT) + # Note that d/dt rho = -d/dx rho = -d/dy rho. + + du1, du2, du3, du4, du5, du6 = source_terms_moist_bubble(u, x, t, equations) + + du1 += rho_x + du2 += rho_x + dp + du3 += du2 + du4 += dE + 2 * dp + du5 += qv * du1 + du6 += ql * du1 + + return SVector(du1, du2, du3, du4, du5, du6) + end + + # Initial condition for the convergence analysis. + # Extends the convergence test to the system with gravity and phase change. + function initial_condition_convergence_test_moist(x, t, + equations::CompressibleMoistEulerEquations2D) + @unpack R_d, R_v, c_vd, c_vv, c_pl, L_00 = equations + + c = 2 + A = 0.1 + L = 2 + f = 1 / L + ω = 2 * pi * f + ini = c + A * sin(ω * (x[1] + x[2] - t)) + + qv = 100 / L_00 + ql = qv / 10 + + mu = ((1 - qv - ql) * c_vd + qv * c_vv + ql * c_pl) + + T = (ini - 1) / mu + 10 / c_vd + 40 + E = (mu * T + qv * L_00 + 1) + + rho = ini + rho_v1 = ini + rho_v2 = ini + rho_e = E * ini + rho_qv = qv * ini + rho_ql = ql * ini + + return SVector(rho, rho_v1, rho_v2, rho_e, rho_qv, rho_ql) + end + + # Source term for the convergence analysis. + # Extends the convergence test to the system with gravity and phase change. + @inline function source_terms_convergence_test_moist(u, x, t, + equations::CompressibleMoistEulerEquations2D) + # Same settings as in `initial_condition` + @unpack R_d, R_v, c_vd, c_vv, c_pl, L_00 = equations + c = 2 + A = 0.1 + L = 2 + f = 1 / L + ω = 2 * pi * f + + x1, x2 = x + si, co = sincos(ω * (x1 + x2 - t)) + rho = c + A * si + rho_x = ω * A * co + + qv = 100 / L_00 + ql = qv / 10 + mu = ((1 - qv - ql) * c_vd + qv * c_vv + ql * c_pl) + xi = ((1 - qv - ql) * R_d + qv * R_v) + + T = (rho - 1) / mu + 10 / c_vd + 40 + dT = rho_x / c_vd + E = (mu * T + qv * L_00 + 1) + dE = E * rho_x + rho * mu * dT + dp = xi * (T * rho_x + rho * dT) + + # Calculate Error in Sources with exact solution and u + u_exact = SVector(rho, rho, rho, rho * E, rho * qv, rho * ql) + + du1, du2, du3, du4, du5, du6 = (source_terms_moist_bubble(u, x, t, equations) - + source_terms_moist_bubble(u_exact, x, t, equations)) + # du1, du2, du3, du4, du5, du6 = zeros(Float64, 6) + # Note that d/dt rho = -d/dx rho = -d/dy rho. + + du1 += rho_x + du2 += rho_x + dp + du3 += rho_x + dp + du4 += dE + 2 * dp + du5 += qv * rho_x + du6 += ql * rho_x + + return SVector(du1, du2, du3, du4, du5, du6) + end + + # Gravity source term + @inline function source_terms_geopotential(u, x, t, + equations::CompressibleMoistEulerEquations2D) + @unpack g = equations + rho, rho_v1, rho_v2, rho_e, rho_qv, rho_ql = u + tmp = rho_v2 + + return SVector(zero(eltype(u)), zero(eltype(u)), + -g * rho, -g * tmp, + zero(eltype(u)), zero(eltype(u))) + end + + @inline function source_terms_geopotential(u, + equations::CompressibleMoistEulerEquations2D) + @unpack g = equations + rho, rho_v1, rho_v2, rho_e, rho_qv, rho_ql = u + tmp = rho_v2 + + return SVector(zero(eltype(u)), zero(eltype(u)), + -g * rho, -g * tmp, + zero(eltype(u)), zero(eltype(u))) + end + + # Rayleigh damping sponge source term form A. Sridhar et al., + # Large-eddy simulations with ClimateMachine: a new open-sourcecode for + # atmospheric simulations on GPUs and CPUs, 2 Oct 2021, doi: 10.5194/gmd-15-6259-2022, + # https://arxiv.org/abs/2110.00853 [physics.ao-ph] . + @inline function source_terms_nonhydrostatic_rayleigh_sponge(u, x, t, + equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_e, rho_qv, rho_ql = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + z = x[2] + + # relaxed background velocity + vr1, vr2 = (10.0, 0.0) + # damping threshold + z_s = 9000.0 + # boundary top + z_top = 16000.0 + # positive even power with default value 2 + gamma = 2.0 + # relaxation coefficient > 0 + alpha = 0.5 + + tau_s = zero(eltype(u)) + if z > z_s + tau_s = alpha * sin(0.5 * (z - z_s) * inv(z_top - z_s))^(gamma) + end + + return SVector(zero(eltype(u)), + -tau_s * rho * (v1 - vr1), + -tau_s * rho * (v2 - vr2), + zero(eltype(u)), zero(eltype(u)), zero(eltype(u))) + end + + # Source term with gravity and phase change + @inline function source_terms_moist_bubble(u, x, t, + equations::CompressibleMoistEulerEquations2D) + return source_terms_geopotential(u, equations) + + source_terms_phase_change(u, equations) + end + + # Calculate pressure and temperature from a state u. + @inline function get_current_condition(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_vd, R_d, c_vv, c_pv, R_v, c_pl, L_00 = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + rho_qd = rho - rho_qv - rho_ql + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + # Inner energy + rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + + # Absolute temperature + T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) + + # Pressure + p = (rho_qd * R_d + rho_qv * R_v) * T + + return SVector(p, T) + end + + # Calculate Q_ph for a state u. + # This source term models the phase chance between could water and vapor. + @inline function phase_change_term(u, equations::CompressibleMoistEulerEquations2D) + @unpack R_v = equations + rho, _, _, _, rho_qv, rho_ql = u + _, T = get_current_condition(u, equations) + rho_qd = rho - rho_qv - rho_ql + + T_C = T - 273.15 + # saturation vapor pressure + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + + # saturation density of vapor + rho_star_qv = p_vs / (R_v * T) + + # Fisher-Burgmeister-Function + a = rho_star_qv - rho_qv + b = rho - rho_qv - rho_qd + + # saturation control factor + # < 1: stronger saturation effect + # > 1: weaker saturation effect + C = 1 + + return (a + b - sqrt(a^2 + b^2)) * C + end + + # Add the source containing Q_ph + @inline function source_terms_phase_change(u, + equations::CompressibleMoistEulerEquations2D) + Q_ph = phase_change_term(u, equations) + + return SVector(zero(eltype(u)), zero(eltype(u)), zero(eltype(u)), + zero(eltype(u)), Q_ph, -Q_ph) + end + + # Low Mach number approximate Riemann solver (LMARS) from + # X. Chen, N. Andronova, B. Van Leer, J. E. Penner, J. P. Boyd, C. Jablonowski, S. + # Lin, A Control-Volume Model of the Compressible Euler Equations with a Vertical Lagrangian + # Coordinate Monthly Weather Review Vol. 141.7, pages 2526–2544, 2013, + # https://journals.ametsoc.org/view/journals/mwre/141/7/mwr-d-12-00129.1.xml. + @inline function flux_LMARS(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + @unpack a = equations + # Unpack left and right state + rho_ll, rho_v1_ll, rho_v2_ll, rho_e_ll, rho_qv_ll, rho_ql_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_e_rr, rho_qv_rr, rho_ql_rr = u_rr + p_ll, T_ll = get_current_condition(u_ll, equations) + p_rr, T_rr = get_current_condition(u_rr, equations) + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_rr = rho_v2_rr / rho_rr + + v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + # diffusion parameter <= 1 + beta = 1 + + # Compute the necessary interface flux components + norm_ = norm(normal_direction) + + rho = 0.5 * (rho_ll + rho_rr) + p_interface = 0.5 * (p_ll + p_rr) - beta * 0.5 * a * rho * (v_rr - v_ll) / norm_ + v_interface = 0.5 * (v_ll + v_rr) - beta * 1 / (2 * a * rho) * (p_rr - p_ll) * norm_ + + if (v_interface > 0) + f1, f2, f3, f4, f5, f6 = u_ll * v_interface + f4 += p_ll * v_interface + else + f1, f2, f3, f4, f5, f6 = u_rr * v_interface + f4 += p_rr * v_interface + end + + return SVector(f1, + f2 + p_interface * normal_direction[1], + f3 + p_interface * normal_direction[2], + f4, f5, f6) + end + + # Convert conservative variables to primitive. + @inline function cons2prim(u, equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p = get_current_condition(u, equations)[1] + qv = rho_qv / rho + ql = rho_ql / rho + + return SVector(rho, v1, v2, p, qv, ql) + end + + # Convert conservative variables to primitive with + # temperature instead of pressure. + @inline function cons2temp(u, equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + v1 = rho_v1 / rho + v2 = rho_v2 / rho + T = get_current_condition(u, equations)[2] + qv = rho_qv / rho + ql = rho_ql / rho + + return SVector(rho, v1, v2, T, qv, ql) + end + + # Convert conservative variables to entropy + @inline function cons2entropy(u, equations::CompressibleMoistEulerEquations2D) + @unpack R_d, R_v, c_pd, c_pv, c_pl, L_00 = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + p, T = get_current_condition(u, equations) + v_square = v1^2 + v2^2 + rho_qd = rho - rho_qv - rho_ql + + # Workaround if an individual density is zero + # Thermodynamic entropy + s_d = 0 + s_v = 0 + s_l = 0 + + # Thermodynamic entropy + if (rho_qd > 0.0) + s_d = c_pd * log(T) - R_d * log(rho_qd * R_d * T) + end + if (rho_qv > 0.0) + s_v = c_pv * log(T) - R_v * log(rho_qv * R_v * T) + end + if (rho_ql > 0.0) + s_l = c_pl * log(T) + end + + g_d = (c_pd - s_d) * T + g_v = L_00 + (c_pv - s_v) * T + g_l = (c_pl - s_l) * T + + w1 = g_d - 0.5 * v_square + w2 = v1 + w3 = v2 + w4 = -1 + w5 = g_v - g_d + w6 = g_l - g_d + + return inv(T) * SVector(w1, w2, w3, w4, w5, w6) + end + + # Convert primitive to conservative variables. + @inline function prim2cons(prim, equations::CompressibleMoistEulerEquations2D) + rho, v1, v2, p, qv, ql = prim + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_qv = rho * qv + rho_ql = rho * ql + rho_E = p * equations.inv_gamma_minus_one + 0.5 * (rho_v1 * v1 + rho_v2 * v2) + return SVector(rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql) + end + + # Convert conservative variables to primitive with + # potential temperature instead of pressure. + @inline function cons2drypot(u, equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + qv = rho_qv / rho + ql = rho_ql / rho + + pot1 = rho + pot2 = v1 + pot3 = v2 + pot4 = dry_pottemp_thermodynamic(u, equations) + pot5 = qv + pot6 = ql + + return SVector(pot1, pot2, pot3, pot4, pot5, pot6) + end + + # Convert conservative variables to primitive with + # moist potential temperature instead of pressure. + @inline function cons2moistpot(u, equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + qv = rho_qv / rho + ql = rho_ql / rho + + pot1 = rho + pot2 = v1 + pot3 = v2 + pot4 = moist_pottemp_thermodynamic(u, equations) + pot5 = qv + pot6 = ql + + return SVector(pot1, pot2, pot3, pot4, pot5, pot6) + end + + # Convert conservative variables to moisture related variables. + @inline function cons2moist(u, equations::CompressibleMoistEulerEquations2D) + @unpack R_d, R_v, c_pd, c_pv, c_pl, p_0 = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + qv = rho_qv / rho + ql = rho_ql / rho + p, T = get_current_condition(u, equations) + + p_v = rho_qv * R_v * T + T_C = T - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + H = p_v * inv(p_vs) + + rho_d = rho - (rho_qv + rho_ql) + r_v = inv(rho_d) * rho_qv + r_l = inv(rho_d) * rho_ql + + # Potential temperature + R_m = R_d + r_v * R_v + c_pml = c_pd + r_v * c_pv + r_l * c_pl + kappa_m = R_m * inv(c_pml) + pot = T * (p_0 / p)^(kappa_m) + + pot1 = qv + pot2 = ql + pot3 = r_v + r_l + pot4 = T + pot5 = H + pot6 = equivalent_pottemp_thermodynamic(u, equations) + + return SVector(pot1, pot2, pot3, pot4, pot5, pot6) + end + + @inline function density(u, equations::CompressibleMoistEulerEquations2D) + rho = u[1] + return rho + end + + @inline function density_dry(u, equations::CompressibleMoistEulerEquations2D) + rho_qd = u[1] - (u[5] + u[6]) + return rho_qd + end + + @inline function density_vapor(u, equations::CompressibleMoistEulerEquations2D) + rho_qv = u[5] + return rho_qv + end + + @inline function temperature(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_vd, R_d, c_vv, c_pv, R_v, c_pl, L_00 = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + rho_qd = rho - rho_qv - rho_ql + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + # inner energy + rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + + # Absolute Temperature + T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) + return T + end + + @inline function saturation_pressure(u, equations::CompressibleMoistEulerEquations2D) + @unpack R_v = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + T = get_current_condition(u, equations)[2] + p_v = rho_qv * R_v * T + T_C = T - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + H = p_v / p_vs + return H + end + + @inline function density_liquid(u, equations::CompressibleMoistEulerEquations2D) + rho_ql = u[6] + return rho_ql + end + + @inline function ratio_liquid(u, equations::CompressibleMoistEulerEquations2D) + rho = u[1] + rho_ql = u[6] + ql = rho_ql / rho + return ql + end + + @inline function ratio_vapor(u, equations::CompressibleMoistEulerEquations2D) + rho = u[1] + rho_qv = u[5] + qv = rho_qv / rho + return qv + end + + @inline function pressure(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_vd, R_d, c_vv, c_pv, R_v, c_pl, L_00 = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + rho_qd = rho - rho_qv - rho_ql + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + # inner energy + rho_e = (rho_E - 0.5 * (rho_v1 * v1 + rho_v2 * v2)) + + # Absolute Temperature + T = (rho_e - L_00 * rho_qv) / (rho_qd * c_vd + rho_qv * c_vv + rho_ql * c_pl) + + # Pressure + p = (rho_qd * R_d + rho_qv * R_v) * T + return p + end + + @inline function density_pressure(u, equations::CompressibleMoistEulerEquations2D) + rho = u[1] + rho_times_p = rho * get_current_condition(u, equations)[1] + return rho_times_p + end + + # Calculate thermodynamic entropy for a conservative state `cons`. + # This is the dryspecific entropy multiplied by the dry air density. + @inline function entropy_thermodynamic(cons, + equations::CompressibleMoistEulerEquations2D) + @unpack c_vd, c_vv, c_pl, R_d, R_v = equations + # Pressure + p, T = get_current_condition(cons, equations) + rho_qd = cons[1] - cons[5] - cons[6] + rho_qv = cons[5] + rho_ql = cons[6] + # Thermodynamic entropy + s_d = c_vd * log(T) - R_d * log(rho_qd * R_d) + + s_v = c_vv * log(T) - R_v * log(rho_qv * R_v) + + s_l = c_pl * log(T) + + return rho_qd * s_d + rho_qv * s_v + rho_ql * s_l + end + + # Calculate mathematical entropy for a conservative state `cons`. + @inline function entropy_math(cons, equations::CompressibleMoistEulerEquations2D) + # Mathematical entropy + S = -entropy_thermodynamic(cons, equations) + + return S + end + + # Default entropy is the mathematical entropy. + @inline entropy(cons, equations::CompressibleMoistEulerEquations2D) = entropy_math(cons, + equations) + + # Calculate total energy for a conservative state `cons`. + @inline energy_total(cons, ::CompressibleMoistEulerEquations2D) = cons[4] + + # Calculate kinetic energy for a conservative state `cons`. + @inline function energy_kinetic(u, equations::CompressibleMoistEulerEquations2D) + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = u + return (rho_v1^2 + rho_v2^2) / (2 * rho) + end + + # Calculate internal energy for a conservative state `cons`. + @inline function energy_internal(cons, equations::CompressibleMoistEulerEquations2D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) + end + + # Calculate the dry potential temperature for a conservative state `cons`. + @inline function dry_pottemp_thermodynamic(cons, + equations::CompressibleMoistEulerEquations2D) + @unpack R_d, p_0, kappa = equations + # Pressure + p = get_current_condition(cons, equations)[1] + # Potential temperature + pot = p_0 * (p / p_0)^(1 - kappa) / (R_d * cons[1]) + + return pot + end + + # Calculate the moist potential temperature for a conservative state `cons`. + @inline function moist_pottemp_thermodynamic(cons, + equations::CompressibleMoistEulerEquations2D) + @unpack R_d, R_v, c_pd, c_pv, c_pl, p_0 = equations + # Pressure + p, T = get_current_condition(cons, equations) + rho_d = cons[1] - (cons[5] + cons[6]) + r_v = inv(rho_d) * cons[5] + r_l = inv(rho_d) * cons[6] + + # Potential temperature + R_m = R_d + r_v * R_v + c_pml = c_pd + r_v * c_pv + r_l * c_pl + kappa_m = R_m * inv(c_pml) + pot = T * (p_0 / p)^(kappa_m) + return pot + end + + # Calculate the equivalent potential temperature for a conservative state `cons`. + @inline function equivalent_pottemp_thermodynamic(cons, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_pv, c_pl, R_d, R_v, p_0, kappa, L_00 = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = cons + rho_d = rho - rho_qv - rho_ql + p, T = get_current_condition(cons, equations) + p_v = rho_qv * R_v * T + p_d = p - p_v + T_C = T - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + H = p_v / p_vs + r_v = rho_qv / rho_d + r_l = rho_ql / rho_d + r_t = r_v + r_l + L_v = L_00 + (c_pv - c_pl) * T + c_p = c_pd + r_t * c_pl + + # equivalent potential temperature + aeq_pot = (T * (p_0 / p_d)^(R_d / c_p) * H^(-r_v * R_v / c_p) * + exp(L_v * r_v * inv(c_p * T))) + + return aeq_pot + end + + # Convert conservative variables to primitive varuables with + # equivalent potential temperature instead of pressure + # and mixing ratios innstead of specific contend. + @inline function cons2aeqpot(cons, equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_pv, c_pl, R_d, R_v, p_0, L_00 = equations + rho, rho_v1, rho_v2, rho_E, rho_qv, rho_ql = cons + rho_d = rho - rho_qv - rho_ql + p, T = get_current_condition(cons, equations) + p_v = rho_qv * R_v * T + p_d = p - p_v + T_C = T - 273.15 + p_vs = 611.2 * exp(17.62 * T_C / (243.12 + T_C)) + H = p_v / p_vs + r_v = rho_qv / rho_d + r_l = rho_ql / rho_d + r_t = r_v + r_l + L_v = L_00 + (c_pv - c_pl) * T + c_p = c_pd + r_t * c_pl + + # equivalent potential temperature + aeq_pot = (T * (p_0 / p_d)^(R_d / c_p) * H^(-r_v * R_v / c_p) * + exp(L_v * r_v * inv(c_p * T))) + + v1 = rho_v1 / rho + v2 = rho_v2 / rho + + pot1 = rho + pot2 = v1 + pot3 = v2 + pot4 = aeq_pot + pot5 = r_v + pot6 = r_t + return SVector(pot1, pot2, pot3, pot4, pot5, pot6) + end + + @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_pv, c_pl, c_vd, c_vv = equations + rho_ll, v1_ll, v2_ll, p_ll, qv_ll, ql_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, qv_rr, ql_rr = cons2prim(u_rr, equations) + qd_ll = 1 - qv_ll - ql_ll + qd_rr = 1 - qv_rr - ql_rr + # Get the density and gas gamma + gamma_ll = (qd_ll * c_pd + qv_ll * c_pv + ql_ll * c_pl) * + inv(qd_ll * c_vd + qv_ll * c_vv + ql_ll * c_pl) + gamma_rr = (qd_rr * c_pd + qv_rr * c_pv + ql_rr * c_pl) * + inv(qd_rr * c_vd + qv_rr * c_vv + ql_rr * c_pl) + + # Compute the sound speeds on the left and right + v_mag_ll = sqrt(v1_ll^2 + v2_ll^2) + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + v_mag_rr = sqrt(v1_rr^2 + v2_rr^2) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr) + end + + # Adjusted version of LLF dissipation from compressible euler. + @inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_pv, c_pl, c_vd, c_vv = equations + rho_ll, v1_ll, v2_ll, p_ll, qv_ll, ql_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, p_rr, qv_rr, ql_rr = cons2prim(u_rr, equations) + qd_ll = 1 - qv_ll - ql_ll + qd_rr = 1 - qv_rr - ql_rr + # Get the density and gas gamma + gamma_ll = (qd_ll * c_pd + qv_ll * c_pv + ql_ll * c_pl) * + inv(qd_ll * c_vd + qv_ll * c_vv + ql_ll * c_pl) + gamma_rr = (qd_rr * c_pd + qv_rr * c_pv + ql_rr * c_pl) * + inv(qd_rr * c_vd + qv_rr * c_vv + ql_rr * c_pl) + # Calculate normal velocities and sound speed + # left + v_ll = (v1_ll * normal_direction[1] + + + v2_ll * normal_direction[2]) + c_ll = sqrt(gamma_ll * p_ll / rho_ll) + # right + v_rr = (v1_rr * normal_direction[1] + + + v2_rr * normal_direction[2]) + c_rr = sqrt(gamma_rr * p_rr / rho_rr) + + return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction) + end + + # Adjusted version of lambda_max from compressible euler. + @inline function max_abs_speeds(u, equations::CompressibleMoistEulerEquations2D) + @unpack c_pd, c_pv, c_pl, c_vd, c_vv = equations + rho, v1, v2, p, qv, ql = cons2prim(u, equations) + qd = 1 - qv - ql + + gamma = (qd * c_pd + qv * c_pv + ql * c_pl) * inv(qd * c_vd + qv * c_vv + ql * c_pl) + c = sqrt(gamma * p / rho) + + return (abs(v1) + c, abs(v2) + c) + end + + # Adjusted EC flux in a normal direction with R_q=0. This is based on + # A. Gouasmi, K. Duraisamy, S. M. Murman, Formulation of Entropy-Stable schemes for the + # multicomponent compressible Euler equations, 4 Feb 2020, doi:10.1016/j.cma.2020.112912, + # https://arxiv.org/abs/1904.00972 [math.NA]. + @inline function flux_chandrashekar(u_ll, u_rr, normal_direction::AbstractVector, + equations::CompressibleMoistEulerEquations2D) + @unpack R_d, R_v, c_vd, c_vv, c_pl, L_00 = equations + R_q = 0 + # Unpack left and right state + rho_ll, rho_v1_ll, rho_v2_ll, rho_E_ll, rho_qv_ll, rho_ql_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_E_rr, rho_qv_rr, rho_ql_rr = u_rr + + rho_qd_ll = rho_ll - rho_qv_ll - rho_ql_ll + rho_qd_rr = rho_rr - rho_qv_rr - rho_ql_rr + v1_ll = rho_v1_ll / rho_ll + v1_rr = rho_v1_rr / rho_rr + v2_ll = rho_v2_ll / rho_ll + v2_rr = rho_v2_rr / rho_rr + + # inner energy + rho_e_ll = (rho_E_ll - 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll)) + rho_e_rr = (rho_E_rr - 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr)) + + # Absolute Temperature + T_ll = (rho_e_ll - L_00 * rho_qv_ll) / + (rho_qd_ll * c_vd + rho_qv_ll * c_vv + rho_ql_ll * c_pl) + T_rr = (rho_e_rr - L_00 * rho_qv_rr) / + (rho_qd_rr * c_vd + rho_qv_rr * c_vv + rho_ql_rr * c_pl) + + # Compute the necessary mean values + rho_qd_mean = 0 + rho_qv_mean = 0 + rho_ql_mean = 0 + inv_T_mean = 0 + if (!(rho_qd_ll == 0.0) && !(rho_qd_rr == 0.0)) + rho_qd_mean = ln_mean(rho_qd_ll, rho_qd_rr) + end + if (!(rho_qv_ll == 0.0) && !(rho_qv_rr == 0.0)) + rho_qv_mean = ln_mean(rho_qv_ll, rho_qv_rr) + end + if (!(rho_ql_ll == 0.0) && !(rho_ql_rr == 0.0)) + rho_ql_mean = ln_mean(rho_ql_ll, rho_ql_rr) + end + if (!(inv(T_ll) == 0.0) && !(inv(T_rr) == 0.0)) + inv_T_mean = inv_ln_mean(inv(T_ll), inv(T_rr)) + end + + v1_avg = 0.5 * (v1_ll + v1_rr) + v2_avg = 0.5 * (v2_ll + v2_rr) + v1_square_avg = 0.5 * (v1_ll^2 + v1_rr^2) + v2_square_avg = 0.5 * (v2_ll^2 + v2_rr^2) + rho_qd_avg = 0.5 * (rho_qd_ll + rho_qd_rr) + rho_qv_avg = 0.5 * (rho_qv_ll + rho_qv_rr) + rho_ql_avg = 0.5 * (rho_ql_ll + rho_ql_rr) + inv_T_avg = 0.5 * (inv(T_ll) + inv(T_rr)) + v_dot_n_avg = normal_direction[1] * v1_avg + normal_direction[2] * v2_avg + + p_int = inv(inv_T_avg) * (R_d * rho_qd_avg + R_v * rho_qv_avg + R_q * rho_ql_avg) + K_avg = 0.5 * (v1_square_avg + v2_square_avg) + + f_1d = rho_qd_mean * v_dot_n_avg + f_1v = rho_qv_mean * v_dot_n_avg + f_1l = rho_ql_mean * v_dot_n_avg + f1 = f_1d + f_1v + f_1l + f2 = f1 * v1_avg + normal_direction[1] * p_int + f3 = f1 * v2_avg + normal_direction[2] * p_int + f4 = ((c_vd * inv_T_mean - K_avg) * f_1d + + (L_00 + c_vv * inv_T_mean - K_avg) * f_1v + + (c_pl * inv_T_mean - K_avg) * f_1l + v1_avg * f2 + v2_avg * f3) + + return SVector(f1, f2, f3, f4, f_1v, f_1l) + end + + varnames(::typeof(cons2cons), ::CompressibleMoistEulerEquations2D) = ("rho", "rho_v1", + "rho_v2", "rho_E", + "rho_qv", + "rho_ql") + varnames(::typeof(cons2prim), ::CompressibleMoistEulerEquations2D) = ("rho", "v1", "v2", + "p", "qv", "ql") + varnames(::typeof(cons2temp), ::CompressibleMoistEulerEquations2D) = ("rho", "v1", "v2", + "T", "qv", "ql") + varnames(::typeof(cons2drypot), ::CompressibleMoistEulerEquations2D) = ("rho", "v1", + "v2", + "drypottemp", + "qv", "ql") + varnames(::typeof(cons2moistpot), ::CompressibleMoistEulerEquations2D) = ("rho", "v1", + "v2", + "moistpottemp", + "qv", "ql") + varnames(::typeof(cons2moist), ::CompressibleMoistEulerEquations2D) = ("qv", "ql", "rt", + "T", "H", + "aeqpottemp") + varnames(::typeof(cons2aeqpot), ::CompressibleMoistEulerEquations2D) = ("rho", "v1", + "v2", + "aeqpottemp", + "rv", "rt") +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 8b13789..a833960 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -1 +1,5 @@ +using Trixi: AbstractEquations +abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: + AbstractEquations{NDIMS, NVARS} end +include("compressible_moist_euler_2d_lucas.jl") diff --git a/test/Project.toml b/test/Project.toml index b3405d2..640c292 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,9 @@ [deps] +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] -Test = "1" \ No newline at end of file +OrdinaryDiffEq = "6.53.2" +Test = "1" +Trixi = "0.7, 0.8" \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index a8e58d3..95801f0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,24 +1,15 @@ -using TrixiAtmo using Test -# We run tests in parallel with CI jobs setting the `TRIXI_TEST` environment +# Trixi.jl runs tests in parallel with CI jobs setting the `TRIXI_TEST` environment # variable to determine the subset of tests to execute. -# By default, we just run the threaded tests since they are relatively cheap -# and test a good amount of different functionality. -const TRIXI_TEST = get(ENV, "TRIXI_TEST", "all") -const TRIXI_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3) -const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) +# +# We could do the same once we have a lot of tests +const TRIXIATMO_TEST = get(ENV, "TRIXIATMO_TEST", "all") +const TRIXIATMO_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3) +const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) -@time @testset "TrixiAtmo.jl tests" begin - @time if TRIXI_TEST == "all" - @test TrixiAtmo.foo() == true - @test TrixiAtmo.bar() == false - @test TrixiAtmo.baz() isa String - end - - @time if TRIXI_TEST == "all" || TRIXI_TEST == "upstream" - @testset "baz()" begin - @test TrixiAtmo.baz() isa String - end +@time @testset verbose=true showtiming=true "TrixiAtmo.jl tests" begin + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "trixi_consistency" + include("test_trixi_consistency.jl") end end diff --git a/test/test_trixi_consistency.jl b/test/test_trixi_consistency.jl new file mode 100644 index 0000000..b13b5d3 --- /dev/null +++ b/test/test_trixi_consistency.jl @@ -0,0 +1,64 @@ +module TestTrixiConsistency + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = TrixiAtmo.examples_dir() +TRIXI_EXAMPLES_DIR = Trixi.examples_dir() + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@trixiatmo_testset "Dry air consistency check" begin + # Dry air warm bubble test case in Trixi.jl + maxiters = 100 + trixi_elixir = joinpath(TRIXI_EXAMPLES_DIR, "structured_2d_dgsem", + "elixir_euler_warm_bubble.jl") + + # Override maxiter and fluxes + trixi_include(trixi_elixir; + volume_flux = Trixi.flux_chandrashekar, + surface_flux = Trixi.FluxLMARS(360.0), + maxiters = maxiters) + + # Save errors + errors_trixi = Main.analysis_callback(Main.sol) + + # Create an instance of Trixi's equations, just used for dispatch below + equations_trixi = Trixi.CompressibleEulerEquations2D(Main.warm_bubble_setup.gamma) + add_zeros = SVector(zero(eltype(Main.sol)), zero(eltype(Main.sol))) + + # Now use moist equations instead + equations_moist = CompressibleMoistEulerEquations2D() + + # Redefine source terms for CompressibleMoistEulerEquations2D + @inline function (setup::Main.WarmBubbleSetup)(u, x, t, + equations_moist::CompressibleMoistEulerEquations2D) + ret_trixi = setup(u, x, t, equations_trixi) + return vcat(ret_trixi, add_zeros) + end + + # Redefine initial condition for CompressibleMoistEulerEquations2D + @inline function (setup::Main.WarmBubbleSetup)(x, t, + equations_moist::CompressibleMoistEulerEquations2D) + ret_trixi = setup(x, t, equations_trixi) + return vcat(ret_trixi, add_zeros) + end + + # Run again with overrides + trixi_include(trixi_elixir, + equations = equations_moist, + volume_flux = flux_chandrashekar, + surface_flux = flux_LMARS, + maxiters = maxiters) + + errors_atmo = Main.analysis_callback(Main.sol) + + for (error_trixi, error_atmo) in zip(errors_trixi.l2, errors_atmo.l2) + @test isapprox(error_trixi, error_atmo, rtol = 1e-12) + end + for (error_trixi, error_atmo) in zip(errors_trixi.linf, errors_atmo.linf) + @test isapprox(error_trixi, error_atmo, rtol = 3e-12) + end +end +end diff --git a/test/test_trixiatmo.jl b/test/test_trixiatmo.jl new file mode 100644 index 0000000..8aa71a5 --- /dev/null +++ b/test/test_trixiatmo.jl @@ -0,0 +1,51 @@ +# TODO - This is a repetition from Trixi.jl + +using Test: @test +using Trixi: Trixi, examples_dir, trixi_include +import TrixiAtmo + +""" + @trixiatmo_testset "name of the testset" #= code to test #= + +Similar to `@testset`, but wraps the code inside a temporary module to avoid +namespace pollution. It also `include`s this file again to provide the +definition of `@test_trixi_include`. Moreover, it records the execution time +of the testset similarly to [`timed_testset`](@ref). +""" +macro trixiatmo_testset(name, expr) + @assert name isa String + # TODO: `@eval` is evil + # We would like to use + # mod = gensym(name) + # ... + # module $mod + # to create new module names for every test set. However, this is not + # compatible with the dirty hack using `@eval` to get the mapping when + # loading structured, curvilinear meshes. Thus, we need to use a plain + # module name here. + quote + local time_start = time_ns() + @eval module TrixiAtmoTestModule + using Test + using TrixiAtmo + include(@__FILE__) + # We define `EXAMPLES_DIR` in (nearly) all test modules and use it to + # get the path to the elixirs to be tested. However, that's not required + # and we want to fail gracefully if it's not defined. + try + import ..EXAMPLES_DIR + import ..TRIXI_EXAMPLES_DIR + catch + nothing + end + @testset $name $expr + end + local time_stop = time_ns() + if TrixiAtmo.Trixi.mpi_isroot() + flush(stdout) + @info("Testset "*$name*" finished in " + *string(1.0e-9 * (time_stop - time_start))*" seconds.\n") + end + nothing + end +end