Skip to content

Commit

Permalink
Merge remote-tracking branch 'jlchan/jc/fix_parabolic_jacobian' into …
Browse files Browse the repository at this point in the history
…jc/fix_parabolic_jacobian
  • Loading branch information
jlchan committed Oct 10, 2023
2 parents 854badf + 02661cb commit 55b2bd1
Show file tree
Hide file tree
Showing 9 changed files with 257 additions and 335 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Trixi"
uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
version = "0.5.45-pre"
version = "0.5.46-pre"

[deps]
CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
Expand Down
215 changes: 215 additions & 0 deletions examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,215 @@
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal compressible Navier-Stokes equations

prandtl_number() = 0.72
mu() = 0.01

equations = CompressibleEulerEquations2D(1.4)
equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu=mu(), Prandtl=prandtl_number(),
gradient_variables=GradientVariablesPrimitive())

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
volume_integral=VolumeIntegralWeakForm())

coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max = ( 1.0, 1.0) # maximum coordinates (max(x), max(y))

trees_per_dimension = (4, 4)
mesh = P4estMesh(trees_per_dimension,
polydeg=3, initial_refinement_level=2,
coordinates_min=coordinates_min, coordinates_max=coordinates_max,
periodicity=(false, false))

# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion2D`
# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion2D`)
# and by the initial condition (which passes in `CompressibleEulerEquations2D`).
# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000)
function initial_condition_navier_stokes_convergence_test(x, t, equations)
# Amplitude and shift
A = 0.5
c = 2.0

# convenience values for trig. functions
pi_x = pi * x[1]
pi_y = pi * x[2]
pi_t = pi * t

rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
v1 = sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A * (x[2] - 1.0)) ) * cos(pi_t)
v2 = v1
p = rho^2

return prim2cons(SVector(rho, v1, v2, p), equations)
end

@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations)
y = x[2]

# TODO: parabolic
# we currently need to hardcode these parameters until we fix the "combined equation" issue
# see also https://github.com/trixi-framework/Trixi.jl/pull/1160
inv_gamma_minus_one = inv(equations.gamma - 1)
Pr = prandtl_number()
mu_ = mu()

# Same settings as in `initial_condition`
# Amplitude and shift
A = 0.5
c = 2.0

# convenience values for trig. functions
pi_x = pi * x[1]
pi_y = pi * x[2]
pi_t = pi * t

# compute the manufactured solution and all necessary derivatives
rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t)
rho_t = -pi * A * sin(pi_x) * cos(pi_y) * sin(pi_t)
rho_x = pi * A * cos(pi_x) * cos(pi_y) * cos(pi_t)
rho_y = -pi * A * sin(pi_x) * sin(pi_y) * cos(pi_t)
rho_xx = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)
rho_yy = -pi * pi * A * sin(pi_x) * cos(pi_y) * cos(pi_t)

v1 = sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
v1_t = -pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * sin(pi_t)
v1_x = pi * cos(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
v1_y = sin(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
v1_xx = -pi * pi * sin(pi_x) * log(y + 2.0) * (1.0 - exp(-A * (y - 1.0))) * cos(pi_t)
v1_xy = pi * cos(pi_x) * (A * log(y + 2.0) * exp(-A * (y - 1.0)) + (1.0 - exp(-A * (y - 1.0))) / (y + 2.0)) * cos(pi_t)
v1_yy = (sin(pi_x) * ( 2.0 * A * exp(-A * (y - 1.0)) / (y + 2.0)
- A * A * log(y + 2.0) * exp(-A * (y - 1.0))
- (1.0 - exp(-A * (y - 1.0))) / ((y + 2.0) * (y + 2.0))) * cos(pi_t))
v2 = v1
v2_t = v1_t
v2_x = v1_x
v2_y = v1_y
v2_xx = v1_xx
v2_xy = v1_xy
v2_yy = v1_yy

p = rho * rho
p_t = 2.0 * rho * rho_t
p_x = 2.0 * rho * rho_x
p_y = 2.0 * rho * rho_y
p_xx = 2.0 * rho * rho_xx + 2.0 * rho_x * rho_x
p_yy = 2.0 * rho * rho_yy + 2.0 * rho_y * rho_y

# Note this simplifies slightly because the ansatz assumes that v1 = v2
E = p * inv_gamma_minus_one + 0.5 * rho * (v1^2 + v2^2)
E_t = p_t * inv_gamma_minus_one + rho_t * v1^2 + 2.0 * rho * v1 * v1_t
E_x = p_x * inv_gamma_minus_one + rho_x * v1^2 + 2.0 * rho * v1 * v1_x
E_y = p_y * inv_gamma_minus_one + rho_y * v1^2 + 2.0 * rho * v1 * v1_y

# Some convenience constants
T_const = equations.gamma * inv_gamma_minus_one / Pr
inv_rho_cubed = 1.0 / (rho^3)

# compute the source terms
# density equation
du1 = rho_t + rho_x * v1 + rho * v1_x + rho_y * v2 + rho * v2_y

# x-momentum equation
du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2
+ 2.0 * rho * v1 * v1_x
+ rho_y * v1 * v2
+ rho * v1_y * v2
+ rho * v1 * v2_y
# stress tensor from x-direction
- 4.0 / 3.0 * v1_xx * mu_
+ 2.0 / 3.0 * v2_xy * mu_
- v1_yy * mu_
- v2_xy * mu_ )
# y-momentum equation
du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2
+ rho * v1_x * v2
+ rho * v1 * v2_x
+ rho_y * v2^2
+ 2.0 * rho * v2 * v2_y
# stress tensor from y-direction
- v1_xy * mu_
- v2_xx * mu_
- 4.0 / 3.0 * v2_yy * mu_
+ 2.0 / 3.0 * v1_xy * mu_ )
# total energy equation
du4 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x)
+ v2_y * (E + p) + v2 * (E_y + p_y)
# stress tensor and temperature gradient terms from x-direction
- 4.0 / 3.0 * v1_xx * v1 * mu_
+ 2.0 / 3.0 * v2_xy * v1 * mu_
- 4.0 / 3.0 * v1_x * v1_x * mu_
+ 2.0 / 3.0 * v2_y * v1_x * mu_
- v1_xy * v2 * mu_
- v2_xx * v2 * mu_
- v1_y * v2_x * mu_
- v2_x * v2_x * mu_
- T_const * inv_rho_cubed * ( p_xx * rho * rho
- 2.0 * p_x * rho * rho_x
+ 2.0 * p * rho_x * rho_x
- p * rho * rho_xx ) * mu_
# stress tensor and temperature gradient terms from y-direction
- v1_yy * v1 * mu_
- v2_xy * v1 * mu_
- v1_y * v1_y * mu_
- v2_x * v1_y * mu_
- 4.0 / 3.0 * v2_yy * v2 * mu_
+ 2.0 / 3.0 * v1_xy * v2 * mu_
- 4.0 / 3.0 * v2_y * v2_y * mu_
+ 2.0 / 3.0 * v1_x * v2_y * mu_
- T_const * inv_rho_cubed * ( p_yy * rho * rho
- 2.0 * p_y * rho * rho_y
+ 2.0 * p * rho_y * rho_y
- p * rho * rho_yy ) * mu_ )

return SVector(du1, du2, du3, du4)
end

initial_condition = initial_condition_navier_stokes_convergence_test

# BC types
velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:3])
heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom)

boundary_condition_left_right = BoundaryConditionDirichlet(initial_condition_navier_stokes_convergence_test)

# define inviscid boundary conditions
boundary_conditions = Dict(:x_neg => boundary_condition_left_right,
:x_pos => boundary_condition_left_right,
:y_neg => boundary_condition_slip_wall,
:y_pos => boundary_condition_slip_wall)

# define viscous boundary conditions
boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_left_right,
:x_pos => boundary_condition_left_right,
:y_neg => boundary_condition_top_bottom,
:y_pos => boundary_condition_top_bottom)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver;
boundary_conditions=(boundary_conditions, boundary_conditions_parabolic),
source_terms=source_terms_navier_stokes_convergence_test)

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

# Create ODE problem with time span `tspan`
tspan = (0.0, 0.5)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)

###############################################################################
# run the simulation

time_int_tol = 1e-8
sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5,
ode_default_options()..., callback=callbacks)
summary_callback() # print the timer summary

2 changes: 1 addition & 1 deletion src/callbacks_step/analysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -297,7 +297,7 @@ function (analysis_callback::AnalysisCallback)(u_ode, du_ode, integrator, semi)
mpi_println(" " * " " *
" " *
" PID: " * @sprintf("%10.8e s", performance_index))
mpi_println(" #DOF: " * @sprintf("% 14d", ndofs(semi)) *
mpi_println(" #DOFs per field:" * @sprintf("% 14d", ndofs(semi)) *
" " *
" alloc'd memory: " * @sprintf("%14.3f MiB", memory_use))
mpi_println(" #elements: " *
Expand Down
29 changes: 29 additions & 0 deletions src/equations/compressible_navier_stokes_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -456,4 +456,33 @@ end
})
return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4])
end

# Dirichlet Boundary Condition for P4est mesh

@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
u_inner,
normal::AbstractVector,
x, t,
operator_type::Gradient,
equations::CompressibleNavierStokesDiffusion2D{
GradientVariablesPrimitive
})
# BCs are usually specified as conservative variables so we convert them to primitive variables
# because the gradients are assumed to be with respect to the primitive variables
u_boundary = boundary_condition.boundary_value_function(x, t, equations)

return cons2prim(u_boundary, equations)
end

@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
u_inner,
normal::AbstractVector,
x, t,
operator_type::Divergence,
equations::CompressibleNavierStokesDiffusion2D{
GradientVariablesPrimitive
})
# for Dirichlet boundary conditions, we do not impose any conditions on the viscous fluxes
return flux_inner
end
end # @muladd
2 changes: 1 addition & 1 deletion src/semidiscretization/semidiscretization_coupled.jl
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationCoupled)
summary_line(increment_indent(io), "source terms", semi.semis[i].source_terms)
summary_line(increment_indent(io), "solver", solver |> typeof |> nameof)
end
summary_line(io, "total #DOFs", ndofs(semi))
summary_line(io, "total #DOFs per field", ndofs(semi))
summary_footer(io)
end
end
Expand Down
2 changes: 1 addition & 1 deletion src/semidiscretization/semidiscretization_hyperbolic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli

summary_line(io, "source terms", semi.source_terms)
summary_line(io, "solver", semi.solver |> typeof |> nameof)
summary_line(io, "total #DOFs", ndofs(semi))
summary_line(io, "total #DOFs per field", ndofs(semi))
summary_footer(io)
end
end
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -228,7 +228,7 @@ function Base.show(io::IO, ::MIME"text/plain",
summary_line(io, "source terms", semi.source_terms)
summary_line(io, "solver", semi.solver |> typeof |> nameof)
summary_line(io, "parabolic solver", semi.solver_parabolic |> typeof |> nameof)
summary_line(io, "total #DOFs", ndofs(semi))
summary_line(io, "total #DOFs per field", ndofs(semi))
summary_footer(io)
end
end
Expand Down
Loading

0 comments on commit 55b2bd1

Please sign in to comment.