Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement EllipsisNotation where useful #189

Merged
merged 5 commits into from
Sep 21, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ authors = ["Michael Schlottke-Lakemper <mschlott@math.uni-koeln.de>", "Gregor Ga
version = "0.2.5-pre"

[deps]
EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
Expand Down
1 change: 1 addition & 0 deletions src/Trixi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ using Printf: @printf, @sprintf, println
using Profile: clear_malloc_data
using Random: seed!

using EllipsisNotation
using HDF5: h5open, attrs
using StaticArrays: @MVector, @SVector, MVector, MMatrix, MArray, SVector, SMatrix, SArray
using TimerOutputs: @notimeit, @timeit, TimerOutput, print_timer, reset_timer!
Expand Down
24 changes: 3 additions & 21 deletions src/io/io.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,13 +38,7 @@ function load_restart_file!(dg::AbstractDg, restart_filename)

# Read variable
println("Reading variables_$v ($name)...")
if ndims(dg) == 2
dg.elements.u[v, :, :, :] = read(file["variables_$v"])
elseif ndims(dg) == 3
dg.elements.u[v, :, :, :, :] = read(file["variables_$v"])
else
error("Unsupported number of spatial dimensions: ", ndims(dg))
end
dg.elements.u[v, .., :] = read(file["variables_$v"])
end
end

Expand Down Expand Up @@ -88,13 +82,7 @@ function save_restart_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep)
# Store each variable of the solution
for v in 1:nvariables(dg)
# Convert to 1D array
if ndims(dg) == 2
file["variables_$v"] = vec(data[v, :, :, :])
elseif ndims(dg) == 3
file["variables_$v"] = vec(data[v, :, :, :, :])
else
error("Unsupported number of spatial dimensions: ", ndims(dg))
end
file["variables_$v"] = vec(data[v, .., :])

# Add variable name as attribute
var = file["variables_$v"]
Expand Down Expand Up @@ -168,13 +156,7 @@ function save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep,
# Store each variable of the solution
for v in 1:nvariables(dg)
# Convert to 1D array
if ndims(dg) == 2
file["variables_$v"] = vec(data[v, :, :, :])
elseif ndims(dg) == 3
file["variables_$v"] = vec(data[v, :, :, :, :])
else
error("Unsupported number of spatial dimensions: ", ndims(dg))
end
file["variables_$v"] = vec(data[v, .., :])

# Add variable name as attribute
var = file["variables_$v"]
Expand Down
16 changes: 3 additions & 13 deletions src/run.jl
Original file line number Diff line number Diff line change
Expand Up @@ -309,18 +309,8 @@ function run_simulation(mesh, solver, time_parameters, time_integration_function
end

# Check steady-state integration residual
if solver.equations isa HyperbolicDiffusionEquations2D
if maximum(abs, view(solver.elements.u_t, 1, :, :, :)) <= solver.equations.resid_tol
println()
println("-"^80)
println(" Steady state tolerance of ",solver.equations.resid_tol," reached at time ",time)
println("-"^80)
println()
finalstep = true
end
end
if solver.equations isa HyperbolicDiffusionEquations3D
if maximum(abs, view(solver.elements.u_t, 1, :, :, :, :)) <= solver.equations.resid_tol
if solver.equations isa AbstractHyperbolicDiffusionEquations
if maximum(abs, view(solver.elements.u_t, 1, .., :)) <= solver.equations.resid_tol
println()
println("-"^80)
println(" Steady state tolerance of ",solver.equations.resid_tol," reached at time ",time)
Expand Down Expand Up @@ -592,7 +582,7 @@ function compute_jacobian_dg(parameters_file; verbose=false, parameters...)
res_m = vec(dg.elements.u_t) |> copy
# restore linearisation state
dg.elements.u[idx] = u0[idx]
# central second order finite difference
# central second order finite difference
J[:,idx] = (res_p - res_m) / (2 * epsilon)
end

Expand Down
28 changes: 6 additions & 22 deletions src/timedisc/timedisc_euler_gravity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -101,18 +101,10 @@ function update_gravity!(solver, u_euler, gravity_parameters)
end

# this is an absolute tolerance check
if ndims(solver) == 2
if maximum(abs, @views solver.elements.u_t[1, :, :, :]) <= solver.equations.resid_tol
# println(" Gravity solution tolerance ", solver.equations.resid_tol,
# " reached in iterations ",iteration)
finalstep = true
end
else
if maximum(abs, @views solver.elements.u_t[1, :, :, :, :]) <= solver.equations.resid_tol
# println(" Gravity solution tolerance ", solver.equations.resid_tol,
# " reached in iterations ",iteration)
finalstep = true
end
if maximum(abs, @views solver.elements.u_t[1, .., :]) <= solver.equations.resid_tol
# println(" Gravity solution tolerance ", solver.equations.resid_tol,
# " reached in iterations ",iteration)
finalstep = true
end
end

Expand All @@ -134,11 +126,7 @@ function timestep_gravity_2N!(solver::AbstractSolver, t, dt, u_euler, gravity_pa

# put in gravity source term proportional to Euler density
# OBS! subtract off the background density ρ_0 (spatial mean value)
if ndims(solver) == 2
@views @. solver.elements.u_t[1,:,:,:] += grav_scale*(u_euler[1,:,:,:] - rho0)
else
@views @. solver.elements.u_t[1,:,:,:,:] += grav_scale*(u_euler[1,:,:,:,:] - rho0)
end
@views @. solver.elements.u_t[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0)

# now take the RK step
@timeit timer() "Runge-Kutta step" begin
Expand Down Expand Up @@ -178,11 +166,7 @@ function timestep_gravity_3Sstar!(solver::AbstractSolver, t, dt, u_euler, gravit
# Source term: Jeans instability OR coupling convergence test OR blast wave
# put in gravity source term proportional to Euler density
# OBS! subtract off the background density ρ_0 around which the Jeans instability is perturbed
if ndims(solver) == 2
@views @. solver.elements.u_t[1,:,:,:] += grav_scale*(u_euler[1,:,:,:] - rho0)
else
@views @. solver.elements.u_t[1,:,:,:,:] += grav_scale*(u_euler[1,:,:,:,:] - rho0)
end
@views @. solver.elements.u_t[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0)

delta_stage = delta[stage]
gamma1_stage = gamma1[stage]
Expand Down