Skip to content

Commit

Permalink
Corrected residual for 3d Alfven wave for the visco-resistive MHD con…
Browse files Browse the repository at this point in the history 
…vergence

test.
  • Loading branch information
@SimonCan
SimonCan committed Feb 16, 2024
1 parent 4549ce8 commit 3421b6e
Showing 1 changed file with 16 additions and 40 deletions.
56 changes: 16 additions & 40 deletions examples/tree_3d_dgsem/elixir_mhd_diffusive_convergence.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@ using Trixi
# semidiscretization of the visco-resistive compressible MHD equations

prandtl_number() = 0.72
mu_const = 2e-2
eta_const = 2e-2
mu_const = 1e-2
eta_const = 1e-2
prandtl_const = prandtl_number()

equations = IdealGlmMhdEquations3D(5 / 3)
Expand All @@ -17,19 +17,18 @@ equations_parabolic = ViscoResistiveMhd3D(equations, mu = mu_const,

volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
solver = DGSEM(polydeg = 3,
# surface_flux = (flux_hindenlang_gassner, flux_nonconservative_powell),
surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
surface_flux = (flux_hindenlang_gassner, 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 = 50_000) # set maximum capacity of tree data structure
initial_refinement_level = 1,
n_cells_max = 150_000) # set maximum capacity of tree data structure

function initial_condition_constant_alfven(x, t, equations)
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.
Expand All @@ -44,24 +43,12 @@ function initial_condition_constant_alfven(x, t, equations)
Va = omega / (ny * sqr)
phi_alv = omega / ny * (nx * (x[1] - 0.5 * r) + ny * (x[2] - 0.5 * r)) - Va * t

# # 1d Alfven wave.
# k = 2*pi
# rho = 1.0
# v1 = 0
# v2 = -e*sin(k*x[1])*sqrt(rho)
# v3 = 0
# B1 = 1
# B2 = e*sin(k*x[1])
# B3 = 0
# p = 1
# psi = 0

# 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
p = 1.0# + e*cos(phi_alv + 1.0)
B1 = nx - rho * v1 * sqr
B2 = ny - rho * v2 * sqr
B3 = -rho * v3 * sqr
Expand All @@ -70,22 +57,10 @@ function initial_condition_constant_alfven(x, t, equations)
return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
end

@inline function source_terms_mhd_convergence_test(u, x, t, equations)
# # 1d Alfven wave
# r_1 = 0.0
# r_2 = 0.0004*pi*sin(4*pi*x[1])
# r_3 = -pi*(0.08*pi*mu_const*sin(2*pi*x[1]) + 0.04*cos(2*pi*x[1]))
# r_4 = 0.0
# r_5 = pi*(0.0016*pi*eta_const*sin(2*pi*x[1])^2 - 0.0016*pi*eta_const*cos(2*pi*x[1])^2 + 0.0127111111111111*pi*mu_const*sin(2*pi*x[1])^2 - 0.0127111111111111*pi*mu_const*cos(2*pi*x[1])^2 + 0.0008*sin(4*pi*x[1]))
# r_6 = 0.0
# r_7 = pi*(0.08*pi*eta_const*sin(2*pi*x[1]) + 0.04*cos(2*pi*x[1]))
# r_8 = 0.0
# r_9 = 0.0

# 3d Alfven wave
r_1 = -0.016*sqrt(5)*pi*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + 0.02*sqrt(5)*pi*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)
@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 =-mu_const*(0.04*sqrt(5)*pi^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1) - 0.0016*sqrt(5)*pi^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 - 0.0008*sqrt(5)*pi^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 - 3.2e-5*sqrt(5)*pi^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^3) - 0.0016*sqrt(5)*pi*(-0.01*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + sqrt(5))*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.008*sqrt(5)*pi*(-0.008*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - sqrt(5)/5)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.016*sqrt(5)*pi*(-0.004*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + 2*sqrt(5)/5)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.0016*sqrt(5)*pi*(0.04*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + sqrt(5))*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.0004*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 + 8.0e-6*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^3 + 8.0e-6*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 + 0.04*pi*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))
r_2 = -mu_const*(0.04*sqrt(5)*pi^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1) - 0.0016*sqrt(5)*pi^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 - 0.0008*sqrt(5)*pi^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 - 3.2e-5*sqrt(5)*pi^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^3) - 0.0016*sqrt(5)*pi*(-0.01*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + sqrt(5))*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.008*sqrt(5)*pi*(-0.008*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - sqrt(5)/5)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.016*sqrt(5)*pi*(-0.004*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + 2*sqrt(5)/5)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.0016*sqrt(5)*pi*(0.04*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + sqrt(5))*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.0004*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 + 8.0e-6*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^3 + 8.0e-6*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 + 0.04*pi*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))

r_3 = -mu_const*(-0.02*sqrt(5)*pi^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1) + 0.0008*sqrt(5)*pi^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 + 0.0004*sqrt(5)*pi^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 + 1.6e-5*sqrt(5)*pi^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)^2*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^3) + 0.0032*sqrt(5)*pi*(-0.01*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + sqrt(5))*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.004*sqrt(5)*pi*(-0.008*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - sqrt(5)/5)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.008*sqrt(5)*pi*(-0.004*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + 2*sqrt(5)/5)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + 0.0032*sqrt(5)*pi*(0.04*sqrt(5)*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) + sqrt(5))*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5)) - 0.0008*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))*cos(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 + 1.6e-5*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)^2*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^3 + 1.6e-5*pi*(-0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) - 1)*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))^2*sin(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1)/(0.02*cos(-sqrt(5)*pi*t + sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5) + 1) + 1)^2 - 0.02*pi*sin(sqrt(5)*pi*t - sqrt(5)*pi*(sqrt(5)*(x[1] - 1.0)/5 + 2*sqrt(5)*(x[2] - 1.0)/5))

Expand All @@ -104,17 +79,18 @@ end
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
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_mhd_convergence_test)
source_terms = source_terms)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span `tspan`
tspan = (0.0, 0.01)
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
Expand All @@ -125,11 +101,11 @@ summary_callback = SummaryCallback()
analysis_callback = AnalysisCallback(semi, interval = 100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval = 100,
save_solution = SaveSolutionCallback(interval = 200,
solution_variables = cons2prim)

# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
cfl = 0.1
cfl = 0.5
stepsize_callback = StepsizeCallback(cfl = cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
Expand Down

0 comments on commit 3421b6e

Please sign in to comment.