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.

Sign up for GitHub

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

Factorize matrices for BBMBBMVariableEquations1D #112

Merged
merged 4 commits into from
May 23, 2024
Merged

Factorize matrices for BBMBBMVariableEquations1D #112

Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ee5d61f
factorize matrices for BBMBBMVariableEquations1D
JoshuaLampert May 22, 2024
5bf3378
Update src/equations/bbm_bbm_variable_bathymetry_1d.jl
JoshuaLampert May 23, 2024
b999abd
adjust reference values
JoshuaLampert May 23, 2024
3a07113
Merge branch 'factorize-matrices-variable-bbmbbm' of https://github.c…
JoshuaLampert May 23, 2024
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
49 changes: 35 additions & 14 deletions src/equations/bbm_bbm_variable_bathymetry_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,14 +202,15 @@ function create_cache(mesh, equations::BBMBBMVariableEquations1D,
K = Diagonal(D .^ 2)
if solver.D1 isa PeriodicDerivativeOperator ||
solver.D1 isa UniformPeriodicCoupledOperator
invImDKD = inv(I - 1 / 6 * Matrix(solver.D1) * K * Matrix(solver.D1))
invImDKD = lu(I - 1 / 6 * sparse(solver.D1) * K * sparse(solver.D1))
JoshuaLampert marked this conversation as resolved.
Show resolved Hide resolved
elseif solver.D1 isa PeriodicUpwindOperators
invImDKD = inv(I - 1 / 6 * Matrix(solver.D1.minus) * K * Matrix(solver.D1.plus))
invImDKD = lu(I - 1 / 6 * sparse(solver.D1.minus) * K * sparse(solver.D1.plus))
else
@error "unknown type of first-derivative operator: $(typeof(solver.D1))"
end
invImD2K = inv(I - 1 / 6 * Matrix(solver.D2) * K)
return (invImDKD = invImDKD, invImD2K = invImD2K, D = D)
invImD2K = lu(I - 1 / 6 * sparse(solver.D2) * K)
tmp2 = Array{RealT}(undef, nnodes(mesh))
return (invImDKD = invImDKD, invImD2K = invImD2K, D = D, tmp2 = tmp2)
end

function create_cache(mesh, equations::BBMBBMVariableEquations1D,
Expand All @@ -229,7 +230,7 @@ function create_cache(mesh, equations::BBMBBMVariableEquations1D,
Pd = BandedMatrix((-1 => fill(one(real(mesh)), N - 2),), (N, N - 2))
D2d = (sparse(solver.D2) * Pd)[2:(end - 1), :]
# homogeneous Dirichlet boundary conditions
invImD2Kd = inv(I - 1 / 6 * D2d * K_i)
invImD2Kd = lu(I - 1 / 6 * D2d * K_i)
m = diag(M)
m[1] = 0
m[end] = 0
Expand All @@ -239,14 +240,16 @@ function create_cache(mesh, equations::BBMBBMVariableEquations1D,
if solver.D1 isa DerivativeOperator ||
solver.D1 isa UniformCoupledOperator
D1_b = BandedMatrix(solver.D1)
invImDKDn = inv(I + 1 / 6 * inv(M) * D1_b' * PdM * K * D1_b)
invImDKDn = lu(I + 1 / 6 * inv(M) * D1_b' * PdM * K * D1_b)
elseif solver.D1 isa UpwindOperators
D1plus_b = BandedMatrix(solver.D1.plus)
invImDKDn = inv(I + 1 / 6 * inv(M) * D1plus_b' * PdM * K * D1plus_b)
invImDKDn = lu(I + 1 / 6 * inv(M) * D1plus_b' * PdM * K * D1plus_b)
else
@error "unknown type of first-derivative operator: $(typeof(solver.D1))"
end
return (invImD2Kd = invImD2Kd, invImDKDn = invImDKDn, D = D)
tmp2 = Array{RealT}(undef, nnodes(mesh))
tmp3 = Array{RealT}(undef, nnodes(mesh) - 2)
return (invImD2Kd = invImD2Kd, invImDKDn = invImDKDn, D = D, tmp2 = tmp2, tmp3 = tmp3)
end

# Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
Expand All @@ -256,7 +259,7 @@ end
function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
initial_condition, ::BoundaryConditionPeriodic, source_terms,
solver, cache)
@unpack invImDKD, invImD2K, D = cache
@unpack invImDKD, invImD2K, D, tmp1, tmp2 = cache

q = wrap_array(u_ode, mesh, equations, solver)
dq = wrap_array(du_ode, mesh, equations, solver)
Expand All @@ -283,16 +286,26 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,

@timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver)

@timeit timer() "deta elliptic" deta[:]=invImDKD * deta
@timeit timer() "dv elliptic" dv[:]=invImD2K * dv
# To use the in-place version `ldiv!` instead of `\`, we need temporary arrays
# since `deta` and `dv` are not stored contiguously
@timeit timer() "deta elliptic" begin
tmp1[:] = deta
ldiv!(tmp2, invImDKD, tmp1)
deta[:] = tmp2
end
@timeit timer() "dv elliptic" begin
tmp2[:] = dv
ldiv!(tmp1, invImD2K, tmp2)
dv[:] = tmp1
end

return nothing
end

function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
initial_condition, ::BoundaryConditionReflecting, source_terms,
solver, cache)
@unpack invImDKDn, invImD2Kd, D = cache
@unpack invImDKDn, invImD2Kd, D, tmp1, tmp2, tmp3 = cache

q = wrap_array(u_ode, mesh, equations, solver)
dq = wrap_array(du_ode, mesh, equations, solver)
Expand Down Expand Up @@ -320,10 +333,18 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,

@timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver)

@timeit timer() "deta elliptic" deta[:]=invImDKDn * deta
# To use the in-place version `ldiv!` instead of `\`, we need temporary arrays
# since `deta` and `dv` are not stored contiguously
@timeit timer() "deta elliptic" begin
tmp1[:] = deta
ldiv!(tmp2, invImDKDn, tmp1)
deta[:] = tmp2
end
@timeit timer() "dv elliptic" begin
dv[2:(end - 1)] = invImD2Kd * dv[2:(end - 1)]
tmp2[:] = dv
ldiv!(tmp3, invImD2Kd, tmp2[2:(end - 1)])
dv[1] = dv[end] = zero(eltype(dv))
dv[2:(end - 1)] = tmp3
end

return nothing
Expand Down
Loading