diff --git a/examples/tree_3d_dgsem/elixir_mhd_diffusive_alfven_wave.jl b/examples/tree_3d_dgsem/elixir_mhd_diffusive_alfven_wave.jl new file mode 100644 index 00000000000..a1e55fe367b --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_mhd_diffusive_alfven_wave.jl @@ -0,0 +1,101 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the visco-resistive compressible MHD equations + +prandtl_number() = 0.72 +mu_const = 2e-2 +eta_const = 2e-2 + +equations = IdealGlmMhdEquations3D(5 / 3) +equations_parabolic = ViscoResistiveMhd3D(equations, mu = mu_const, + Prandtl = prandtl_number(), + eta = eta_const, + gradient_variables = GradientVariablesPrimitive()) + +volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# Create a uniformly refined mesh +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + n_cells_max = 200_000) # set maximum capacity of tree data structure + +function initial_condition_constant_alfven(x, t, equations) + # Alfvén wave in three space dimensions modified by a periodic density variation. + # For the system without the density variations see: Altmann thesis http://dx.doi.org/10.18419/opus-3895. + # Domain must be set to [-1, 1]^3, γ = 5/3. + omega = 2.0 * pi # may be multiplied by frequency + # r = length-variable = length of computational domain + r = 2.0 + # e = epsilon + e = 0.02 + nx = 1 / sqrt(r^2 + 1.0) + ny = r / sqrt(r^2 + 1.0) + sqr = 1.0 + Va = omega / (ny * sqr) + phi_alv = omega / ny * (nx * (x[1] - 0.5 * r) + ny * (x[2] - 0.5 * r)) - Va * t + + rho = 1.0 + e * cos(phi_alv + 1.0) + v1 = -e * ny * cos(phi_alv) / rho + v2 = e * nx * cos(phi_alv) / rho + v3 = e * sin(phi_alv) / rho + p = 1.0 + B1 = nx - rho * v1 * sqr + B2 = ny - rho * v2 * sqr + B3 = -rho * v3 * sqr + psi = 0.0 + + return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) +end + +initial_condition = initial_condition_constant_alfven + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 0.1) +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 100, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +cfl = 0.5 +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + save_solution, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1e-5, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary. +summary_callback() diff --git a/examples/tree_3d_dgsem/elixir_mhd_diffusive_convergence.jl b/examples/tree_3d_dgsem/elixir_mhd_diffusive_convergence.jl new file mode 100644 index 00000000000..c6a56b94783 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_mhd_diffusive_convergence.jl @@ -0,0 +1,133 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the visco-resistive compressible MHD equations + +prandtl_number() = 0.72 +mu_const = 1e-3 +eta_const = 1e-3 +prandtl_const = prandtl_number() + +equations = IdealGlmMhdEquations3D(5 / 3) +equations_parabolic = ViscoResistiveMhd3D(equations, mu = mu_const, + Prandtl = prandtl_number(), + eta = eta_const, + gradient_variables = GradientVariablesPrimitive()) + +# volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) +# solver = DGSEM(polydeg = 3, +# surface_flux = (flux_hindenlang_gassner, flux_nonconservative_powell), +# volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) +volume_flux = (flux_central, flux_nonconservative_powell) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +# Create a uniformly refined mesh +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 1, + n_cells_max = 150_000) # set maximum capacity of tree data structure + +function initial_condition_constant_alfven_3d(x, t, equations) + # Alfvén wave in three space dimensions modified by a periodic density variation. + # For the system without the density variations see: Altmann thesis http://dx.doi.org/10.18419/opus-3895. + # Domain must be set to [-1, 1]^3, γ = 5/3. + omega = 2.0 * pi # may be multiplied by frequency + # r = length-variable = length of computational domain + r = 2.0 + # e = epsilon + e = 0.02 + nx = 1 / sqrt(r^2 + 1.0) + ny = r / sqrt(r^2 + 1.0) + sqr = 1.0 + Va = omega / (ny * sqr) + phi_alv = omega / ny * (nx * (x[1] - 0.5 * r) + ny * (x[2] - 0.5 * r)) - Va * t + + # 3d Alfven wave + rho = 1.0 + e * cos(phi_alv + 1.0) + v1 = -e * ny * cos(phi_alv) / rho + v2 = e * nx * cos(phi_alv) / rho + v3 = e * sin(phi_alv) / rho + p = 1.0# + e*cos(phi_alv + 1.0) + B1 = nx - rho * v1 * sqr + B2 = ny - rho * v2 * sqr + B3 = -rho * v3 * sqr + psi = 0.0 + + return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) +end + +@inline function source_terms_mhd_convergence_test_3d(u, x, t, equations) + r_1 = 0.02*sqrt(5)*pi*sin(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + + r_2 = sqrt(5)*pi^2*mu_const*(-0.04*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0)) + (0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)*(0.0012*cos(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 1) - 0.0004*cos(1)) + 3.2e-5*sin(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0)))/(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^3 + + r_3 = -sqrt(5)*pi^2*mu_const*(-0.02*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0)) + (0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)*(0.0006*cos(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 1) - 0.0002*cos(1)) + 1.6e-5*sin(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0)))/(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^3 + + r_4 = -pi^2*mu_const*((0.003*sin(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 1) + 0.001*sin(1))*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1) + 0.1*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2])) - 8.0e-5*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))*sin(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1)^2)/(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^3 + +# r_5 = 0.2*pi*(2.77777777777778e-7*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) -2.58888888888889e-5*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 + 1.38888888888889e-7*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^3 - 1.78888888888889e-5*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 + 0.000547222222222222*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) + 2.77777777777778e-7*pi*mu_const*sin(-sqrt(5)*pi*t +pi*x[1] + 2*pi*x[2] + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) - 2.58888888888889e-5*pi*mu_const*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2 + 1.38888888888889e-7*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^3 - 1.78888888888889e-5*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 + 0.000547222222222223*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] -2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) + 1.0e-7*sqrt(5)*(-sin(-4*sqrt(5)*pi*t + 4*pi*x[1] + 8*pi*x[2] + 2) + 2*sin(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 2) - sin(2)) + 1.0e-7*sqrt(5)*(sin(-4*sqrt(5)*pi*t + 4*pi*x[1] + 8*pi*x[2] + 2) + 2*sin(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 2) + sin(2)) - 8.0e-9*sqrt(5)*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 - 2.0e-5*sqrt(5)*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) - 8.0e-9*sqrt(5)*sin(-sqrt(5)*pi*t+ pi*x[1] + 2*pi*x[2] + 1)*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 - 2.0e-5*sqrt(5)*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2)/(-0.04*sin(pi*x[2])*sin(-sqrt(5)*pi*t + pi*x[1] + pi*x[2] + 1) + 0.02*cos(-sqrt(5)*pi*t + pi*x[1] + 1) - 1.0)^4 + + r_5 = 0.2*pi*(1.7205356741103e-21*pi*mu_const*(cos(-4*sqrt(5)*pi*t + 4*pi*x[1] + 8*pi*x[2] + 2) - cos(2)) - 1.01643953670516e-19*pi*mu_const*(cos(-3*sqrt(5)*pi*t + 3*pi*x[1] + 6*pi*x[2] + 1) - cos(sqrt(5)*pi*t - pi*x[1] - 2*pi*x[2] + 1)) + 2.77777777777778e-7*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) - 2.58888888888889e-5*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 + 1.38888888888889e-7*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^3 - 1.78888888888889e-5*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 + 0.000547222222222222*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) + 4.33680868994202e-18*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2 - 1.17643464896431e-22*pi*mu_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 + 2.77777777777778e-7*pi*mu_const*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) - 2.58888888888889e-5*pi*mu_const*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2 - 2.16840434497101e-25*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^4 + 1.38888888888889e-7*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^3 - 1.78888888888889e-5*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 + 0.000547222222222223*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) - 6.64073830647371e-18*pi*mu_const*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2 + 1.0842021724855e-21*sqrt(5)*(-sin(-3*sqrt(5)*pi*t + 3*pi*x[1] + 6*pi*x[2] + 1) + sin(sqrt(5)*pi*t - pi*x[1] - 2*pi*x[2] + 1)) + 1.0842021724855e-27*sqrt(5)*(cos(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 2) + 1)^2*sin(pi*(-2*sqrt(5)*t + 2*x[1] + 4*x[2])) - 1.62630325872826e-23*sqrt(5)*(-2*sin(pi*(-2*sqrt(5)*t + 2*x[1] + 4*x[2])) - sin(-4*sqrt(5)*pi*t + 4*pi*x[1] + 8*pi*x[2] + 2) + sin(2)) + 1.0e-7*sqrt(5)*(-sin(-4*sqrt(5)*pi*t + 4*pi*x[1] + 8*pi*x[2] + 2) + 2*sin(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 2) - sin(2)) + 1.0e-7*sqrt(5)*(sin(-4*sqrt(5)*pi*t + 4*pi*x[1] + 8*pi*x[2] + 2) + 2*sin(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 2) + sin(2)) + 2.71050543121376e-20*sqrt(5)*sin(pi*(-2*sqrt(5)*t + 2*x[1] + 4*x[2])) - 8.0e-9*sqrt(5)*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 - 2.0e-5*sqrt(5)*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))^2*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1) + 1.73472347597681e-24*sqrt(5)*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2]))*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^3 - 8.0e-9*sqrt(5)*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2*cos(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)^2 - 2.0e-5*sqrt(5)*sin(-sqrt(5)*pi*t + pi*x[1] + 2*pi*x[2] + 1)*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0))^2)/(-0.04*sin(pi*x[2])*sin(-sqrt(5)*pi*t + pi*x[1] + pi*x[2] + 1) + 0.02*cos(-sqrt(5)*pi*t + pi*x[1] + 1) - 1.0)^4 + + r_6 = pi*(-0.04*(sqrt(5)*pi*eta_const*cos(pi*(-sqrt(5)*t + x[1] + 2*x[2])) + sin(pi*(-sqrt(5)*t + x[1] + 2*x[2])))*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2 + 0.04*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2])) + 0.0004*sin(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 1) + 0.0004*sin(1))/(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2 + + r_7 = pi*(0.02*(sqrt(5)*pi*eta_const*cos(pi*(-sqrt(5)*t + x[1] + 2*x[2])) + sin(pi*(-sqrt(5)*t + x[1] + 2*x[2])))*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2 - 0.02*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2])) - 0.0002*sin(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 1) - 0.0002*sin(1))/(0.02*cos(-sqrt(5)*pi*t + pi*(x[1]+ 2*x[2] - 3.0) + 1) + 1)^2 + + r_8 = pi*((0.1*pi*eta_const*sin(pi*(-sqrt(5)*t + x[1] + 2*x[2])) + 0.02*sqrt(5)*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0)))*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2 - 0.02*sqrt(5)*(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) +1)*cos(pi*(sqrt(5)*t - x[1] - 2*x[2] + 3.0)) + 0.0002*sqrt(5)*(cos(-2*sqrt(5)*pi*t + 2*pi*x[1] + 4*pi*x[2] + 1) - cos(1)))/(0.02*cos(-sqrt(5)*pi*t + pi*(x[1] + 2*x[2] - 3.0) + 1) + 1)^2 + + r_9 = 0.0 + + return SVector(r_1, r_2, r_3, r_4, r_5, r_6, r_7, r_8, r_9) +end + +initial_condition = initial_condition_constant_alfven_3d +source_terms = source_terms_mhd_convergence_test_3d + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), + initial_condition, solver; + source_terms = source_terms) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.5) +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 200, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +cfl = 0.5 +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + save_solution, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1e-5, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary. +summary_callback() diff --git a/src/Trixi.jl b/src/Trixi.jl index bf0986084af..a839374687c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -146,6 +146,7 @@ export AcousticPerturbationEquations2D, CompressibleEulerMulticomponentEquations2D, CompressibleEulerEquationsQuasi1D, IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D, + ViscoResistiveMhd3D, IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D, diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl index a063e9f2758..d4a0cf9266c 100644 --- a/src/equations/equations_parabolic.jl +++ b/src/equations/equations_parabolic.jl @@ -21,3 +21,8 @@ include("compressible_navier_stokes.jl") include("compressible_navier_stokes_1d.jl") include("compressible_navier_stokes_2d.jl") include("compressible_navier_stokes_3d.jl") + +# Resistive and visuous MHD +abstract type AbstractViscoResistiveMhd{NDIMS, NVARS} <: + AbstractEquationsParabolic{NDIMS, NVARS, GradientVariablesConservative} end +include("visco_resistive_mhd_3d.jl") diff --git a/src/equations/visco_resistive_mhd_3d.jl b/src/equations/visco_resistive_mhd_3d.jl new file mode 100644 index 00000000000..9bed6a4e79d --- /dev/null +++ b/src/equations/visco_resistive_mhd_3d.jl @@ -0,0 +1,210 @@ +# 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 +#! format: noindent + +@doc raw""" + ViscoResistiveMhd3D(gamma, inv_gamma_minus_one, + μ, Pr, eta, kappa, + equations, gradient_variables) + +These equations contain the viscous Navier-Stokes equations coupled to +the magnetic field together with the magnetic diffusion applied +to mass, momenta, magnetic field and total energy together with the advective terms from +the [`IdealGlmMhdEquations3D`](@ref). +Together they constitute the compressible, viscous and resistive MHD equations with energy. + +- `gamma`: adiabatic constant, +- `mu`: dynamic viscosity, +- `Pr`: Prandtl number, +- `eta`: magnetic diffusion (resistivity) +- `equations`: instance of the [`IdealGlmMhdEquations3D`](@ref) +- `gradient_variables`: which variables the gradients are taken with respect to. +Defaults to `GradientVariablesPrimitive()`. + +Fluid properties such as the dynamic viscosity $\mu$ and magnetic diffusion $\eta$ +can be provided in any consistent unit system, e.g., +[$\mu$] = kg m⁻¹ s⁻¹. + +The equations given here are the visco-resistive part of the MHD equations +in conservation form: +```math +\overleftrightarrow{f}^\mathrm{\mu\eta} = +\begin{pmatrix} +\overrightarrow{0} \\ +\underline{\tau} \\ +\underline{\tau}\overrightarrow{v} - \overrightarrow{\nabla}q - + \eta\left( (\overrightarrow{\nabla}\times\overrightarrow{B})\times\overrightarrow{B} \right) \\ +\eta \left( (\overrightarrow{\nabla}\overrightarrow{B})^\mathrm{T} - \overrightarrow{\nabla}\overrightarrow{B} \right) \\ +\overrightarrow{0} +\end{pmatrix}, +``` +where `\tau` is the viscous stress tensor and `q = \kappa \overrightarrow{\nabla} T`. +For the induction term we have the usual Laplace operator on the magnetic field +but we also include terms with `div(B)`. +Divergence cleaning is done using the `\psi` field. + +For more details see e.g. arXiv:2012.12040. +""" +struct ViscoResistiveMhd3D{GradientVariables, RealT <: Real, + E <: AbstractIdealGlmMhdEquations{3}} <: + AbstractViscoResistiveMhd{3, 9} + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + mu::RealT # viscosity + Pr::RealT # Prandtl number + eta::RealT # magnetic diffusion + kappa::RealT # thermal diffusivity for Fick's law + equations_hyperbolic::E # IdealGlmMhdEquations3D + gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy +end + +# default to primitive gradient variables +function ViscoResistiveMhd3D(equations::IdealGlmMhdEquations3D; + mu, Prandtl, eta, + gradient_variables = GradientVariablesPrimitive()) + gamma = equations.gamma + inv_gamma_minus_one = equations.inv_gamma_minus_one + μ, Pr, eta = promote(mu, Prandtl, eta) + + # Under the assumption of constant Prandtl number the thermal conductivity + # constant is kappa = gamma μ / ((gamma-1) Pr). + # Important note! Factor of μ is accounted for later in `flux`. + kappa = gamma * inv_gamma_minus_one / Pr + + ViscoResistiveMhd3D{typeof(gradient_variables), typeof(gamma), typeof(equations) + }(gamma, inv_gamma_minus_one, + μ, Pr, eta, kappa, + equations, gradient_variables) +end + +# Explicit formulas for the diffusive MHD fluxes are available, e.g., in Section 2 +# of the paper by Rueda-Ramírez, Hennemann, Hindenlang, Winters, and Gassner +# "An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive +# MHD Equations. Part II: Subcell Finite Volume Shock Capturing" +function flux(u, gradients, orientation::Integer, equations::ViscoResistiveMhd3D) + # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`. + rho, v1, v2, v3, E, B1, B2, B3, psi = convert_transformed_to_primitive(u, equations) + # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, v3, T) + # either computed directly or reverse engineered from the gradient of the entropy variables + # by way of the `convert_gradient_variables` function. + + @unpack eta = equations + + _, dv1dx, dv2dx, dv3dx, dTdx, dB1dx, dB2dx, dB3dx, _ = convert_derivative_to_primitive(u, + gradients[1], + equations) + _, dv1dy, dv2dy, dv3dy, dTdy, dB1dy, dB2dy, dB3dy, _ = convert_derivative_to_primitive(u, + gradients[2], + equations) + _, dv1dz, dv2dz, dv3dz, dTdz, dB1dz, dB2dz, dB3dz, _ = convert_derivative_to_primitive(u, + gradients[3], + equations) + + # Components of viscous stress tensor + + # Diagonal parts + # (4/3 * (v1)_x - 2/3 * ((v2)_y + (v3)_z) + tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * (dv2dy + dv3dz) + # (4/3 * (v2)_y - 2/3 * ((v1)_x + (v3)_z) + tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * (dv1dx + dv3dz) + # (4/3 * (v3)_z - 2/3 * ((v1)_x + (v2)_y) + tau_33 = 4.0 / 3.0 * dv3dz - 2.0 / 3.0 * (dv1dx + dv2dy) + + # Off diagonal parts, exploit that stress tensor is symmetric + # ((v1)_y + (v2)_x) + tau_12 = dv1dy + dv2dx # = tau_21 + # ((v1)_z + (v3)_x) + tau_13 = dv1dz + dv3dx # = tau_31 + # ((v2)_z + (v3)_y) + tau_23 = dv2dz + dv3dy # = tau_32 + + # Fick's law q = -kappa * grad(T) = -kappa * grad(p / (R rho)) + # with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr) + # Note, the gas constant cancels under this formulation, so it is not present + # in the implementation + q1 = equations.kappa * dTdx + q2 = equations.kappa * dTdy + q3 = equations.kappa * dTdz + + # Constant dynamic viscosity is copied to a variable for readability. + # Offers flexibility for dynamic viscosity via Sutherland's law where it depends + # on temperature and reference values, Ts and Tref such that mu(T) + mu = equations.mu + + if orientation == 1 + # viscous flux components in the x-direction + f1 = zero(rho) + f2 = tau_11 * mu + f3 = tau_12 * mu + f4 = tau_13 * mu + f5 = (v1 * tau_11 + v2 * tau_12 + v3 * tau_13 + q1) * mu + + (B2 * (dB2dx - dB1dy) + B3 * (dB3dx - dB1dz)) * eta + f6 = zero(rho) + f7 = eta * (dB2dx - dB1dy) + f8 = eta * (dB3dx - dB1dz) + f9 = zero(rho) + + return SVector(f1, f2, f3, f4, f5, f6, f7, f8, f9) + elseif orientation == 2 + # viscous flux components in the y-direction + # Note, symmetry is exploited for tau_12 = tau_21 + g1 = zero(rho) + g2 = tau_12 * mu # tau_21 * mu + g3 = tau_22 * mu + g4 = tau_23 * mu + g5 = (v1 * tau_12 + v2 * tau_22 + v3 * tau_23 + q2) * mu + + (B1 * (dB1dy - dB2dx) + B3 * (dB3dy - dB2dz)) * eta + g6 = eta * (dB1dy - dB2dx) + g7 = zero(rho) + g8 = eta * (dB3dy - dB2dz) + g9 = zero(rho) + + return SVector(g1, g2, g3, g4, g5, g6, g7, g8, g9) + else # if orientation == 3 + # viscous flux components in the z-direction + # Note, symmetry is exploited for tau_13 = tau_31, tau_23 = tau_32 + h1 = zero(rho) + h2 = tau_13 * mu # tau_31 * mu + h3 = tau_23 * mu # tau_32 * mu + h4 = tau_33 * mu + h5 = (v1 * tau_13 + v2 * tau_23 + v3 * tau_33 + q3) * mu + + (B1 * (dB1dz - dB3dx) + B2 * (dB2dz - dB3dy)) * eta + h6 = eta * (dB1dz - dB3dx) + h7 = eta * (dB2dz - dB3dy) + h8 = zero(rho) + h9 = zero(rho) + + return SVector(h1, h2, h3, h4, h5, h6, h7, h8, h9) + end +end + +# the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables. +# For CNS, it is simplest to formulate the viscous terms in primitive variables, so we transform the transformed +# variables into primitive variables. +@inline function convert_transformed_to_primitive(u_transformed, + equations::ViscoResistiveMhd3D{ + GradientVariablesPrimitive + }) + return u_transformed +end + +# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4, w_5) and +# reverse engineers the gradients to be terms of the primitive variables (v1, v2, v3, T). +# Helpful because then the diffusive fluxes have the same form as on paper. +# Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused. +# TODO: parabolic; entropy stable viscous terms +@inline function convert_derivative_to_primitive(u, gradient, + ::ViscoResistiveMhd3D{ + GradientVariablesPrimitive + }) + return gradient +end + +# Calculate the magnetic energy for a conservative state `cons'. +@inline function energy_magnetic_mhd(cons, ::ViscoResistiveMhd3D) + return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2) +end +end # @muladd diff --git a/test/Project.toml b/test/Project.toml index a376c2805ea..4c2e7b577cf 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,15 +1,22 @@ [deps] +AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" +Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95" [compat] Aqua = "0.8" diff --git a/test/test_tree_3d_mhd.jl b/test/test_tree_3d_mhd.jl index e75685f0b43..58b9c4db10a 100644 --- a/test/test_tree_3d_mhd.jl +++ b/test/test_tree_3d_mhd.jl @@ -291,6 +291,54 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_mhd_diffusive_alfven_wave.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_diffusive_alfven_wave.jl"), + l2=[ + 0.009723896964276463, + 0.0013111262214223665, + 0.000639632696039887, + 0.0014689482873199698, + 0.00036585344043155426, + 0.0013294252402214933, + 0.0006389697014878424, + 0.0014688601323265813, + 0.00022638993249319084, + ], + linf=[ + 0.014222841419814558, + 0.002288228859293565, + 0.001110821298443975, + 0.0026433825065787994, + 0.0008390577089585349, + 0.0024549005206054852, + 0.0012253189393660602, + 0.0027498258638988405, + 0.0007904248596911075, + ]) + + @testset "analysis_callback(sol) for AnalysisCallbackViscoResistiveMHD" begin + errors = analysis_callback(sol) + @test errors.l2≈[ + 0.009723896964276463, 0.0013111262214223665, 0.000639632696039887, + 0.0014689482873199698, 0.00036585344043155426, 0.0013294252402214933, + 0.0006389697014878424, 0.0014688601323265813, 0.00022638993249319084, + ] rtol=1.0e-4 + @test errors.linf≈[ + 0.014222841419814558, 0.002288228859293565, 0.001110821298443975, + 0.0026433825065787994, 0.0008390577089585349, 0.0024549005206054852, + 0.0012253189393660602, 0.0027498258638988405, 0.0007904248596911075, + ] rtol=1.0e-4 + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end end end # module