From 0e8905246e2d42186450c93edc0f8a4cf6ff9692 Mon Sep 17 00:00:00 2001 From: hnil Date: Mon, 2 Sep 2024 14:21:19 +0200 Subject: [PATCH 01/15] - made it possible to wrap krylov solvers as preconditioneres example in BattMo (solver_as_preconditioner.jl) battmo_dev_hmn --- src/discretization/finite-volume.jl | 1 + src/linsolve/krylov.jl | 40 ++++++++++++++++------------- src/linsolve/precond/amg.jl | 5 ++++ src/linsolve/precond/ilu.jl | 33 +++++++++++++++++++----- src/linsolve/precond/utils.jl | 11 ++++---- src/linsolve/precond/various.jl | 16 +++++++++--- src/linsolve/scalar_cpu.jl | 6 ++--- src/models.jl | 3 ++- src/multimodel/model.jl | 5 ++++ 9 files changed, 82 insertions(+), 38 deletions(-) diff --git a/src/discretization/finite-volume.jl b/src/discretization/finite-volume.jl index f76d15d6..bb697433 100644 --- a/src/discretization/finite-volume.jl +++ b/src/discretization/finite-volume.jl @@ -116,6 +116,7 @@ function compute_half_face_trans(cell_centroids, face_centroids, face_normals, f T_hf[fpos] = T end end + #@assert minimum(T_hf) >= 0 return T_hf end diff --git a/src/linsolve/krylov.jl b/src/linsolve/krylov.jl index cf14d01b..2c27fec2 100644 --- a/src/linsolve/krylov.jl +++ b/src/linsolve/krylov.jl @@ -54,12 +54,14 @@ function Base.show(io::IO, krylov::GenericKrylov) print(io, "Generic Krylov using $(krylov.solver) (ϵₐ=$atol, ϵ=$rtol) with prec = $(typeof(krylov.preconditioner))") end -function preconditioner(krylov::GenericKrylov, sys, model, storage, recorder, float_t = Float64) +function preconditioner(krylov::GenericKrylov, sys, context, model, storage, recorder, side, arg...) M = krylov.preconditioner + Ft = float_type(context) if isnothing(M) op = I else - op = PrecondWrapper(linear_operator(M, float_t)) + linop = linear_operator(M, side, Ft, sys, context, model, storage, recorder, arg...) + op = PrecondWrapper(linop) end return op end @@ -70,27 +72,29 @@ end # end function linear_solve!(sys::LSystem, - krylov::GenericKrylov, - model, - storage = nothing, - dt = nothing, - recorder = ProgressRecorder(), - executor = default_executor(); - dx = sys.dx_buffer, - r = vector_residual(sys), - atol = linear_solver_tolerance(krylov, :absolute), - rtol = linear_solver_tolerance(krylov, :relative), - rtol_nl = linear_solver_tolerance(krylov, :nonlinear_relative), - rtol_relaxed = linear_solver_tolerance(krylov, :relaxed_relative) - ) + krylov::GenericKrylov, + context, + model, + storage = nothing, + dt = nothing, + recorder = ProgressRecorder(), + executor = default_executor(); + dx = sys.dx_buffer, + r = vector_residual(sys), + atol = linear_solver_tolerance(krylov, :absolute), + rtol = linear_solver_tolerance(krylov, :relative), + rtol_nl = linear_solver_tolerance(krylov, :nonlinear_relative), + rtol_relaxed = linear_solver_tolerance(krylov, :relaxed_relative) + ) cfg = krylov.config prec = krylov.preconditioner - Ft = float_type(model.context) + Ft = float_type(context) sys = krylov_scale_system!(sys, krylov, dt) t_prep = @elapsed @tic "prepare" prepare_linear_solve!(sys) op = linear_operator(sys) - t_prec = @elapsed @tic "precond" update_preconditioner!(prec, sys, model, storage, recorder, executor) - prec_op = preconditioner(krylov, sys, model, storage, recorder, Ft) + t_prec = @elapsed @tic "precond" update_preconditioner!(prec, sys, context, model, storage, recorder, executor) + prec_op = preconditioner(krylov, sys, context, model, storage, recorder, :left) + # R = preconditioner(krylov, sys, model, storage, recorder, :right, Ft) v = Int64(cfg.verbose) max_it = cfg.max_iterations min_it = cfg.min_iterations diff --git a/src/linsolve/precond/amg.jl b/src/linsolve/precond/amg.jl index 9cf256af..f54b7735 100644 --- a/src/linsolve/precond/amg.jl +++ b/src/linsolve/precond/amg.jl @@ -21,6 +21,11 @@ end matrix_for_amg(A) = A matrix_for_amg(A::StaticSparsityMatrixCSR) = copy(A.At') +#function apply!(x, prec::AMGPreconditioner, y, float_t, sys, args...) +# apply!(x,prec,y) +#end + + function update_preconditioner!(amg::AMGPreconditioner{flavor}, A, b, context, executor) where flavor kw = amg.method_kwarg A_amg = matrix_for_amg(A) diff --git a/src/linsolve/precond/ilu.jl b/src/linsolve/precond/ilu.jl index 6eb2ff0a..7419c8de 100644 --- a/src/linsolve/precond/ilu.jl +++ b/src/linsolve/precond/ilu.jl @@ -1,4 +1,3 @@ - """ ILU(0) preconditioner on CPU """ @@ -56,12 +55,24 @@ function update_preconditioner!(ilu::ILUZeroPreconditioner, A::StaticSparsityMat end end -function apply!(x, ilu::ILUZeroPreconditioner, y) +function apply!(x, ilu::ILUZeroPreconditioner, y, type, arg...) factor = get_factorization(ilu) - ilu_apply!(x, factor, y) + mytype = :both + ilu_apply!(x, factor, y, mytype, arg...) end -function ilu_apply!(x::AbstractArray{F}, f::AbstractILUFactorization, y::AbstractArray{F}) where {F<:Real} +function ilu_f(type::Symbol) + # Why must this be qualified? + if type == :left + f = forward_substitution! + elseif type == :right + f = backward_substitution! + else + f = ldiv! + end +end + +function ilu_apply!(x::AbstractArray{F}, f::AbstractILUFactorization, y::AbstractArray{F}, type::Symbol = :both, args...) where {F<:Real} T = eltype(f) if T == Float64 ldiv!(x, f, y) @@ -79,11 +90,19 @@ function ilu_apply!(x, f::AbstractILUFactorization, y) ldiv!(x, f, y) end -function ilu_apply!(x::AbstractArray{F}, f::ILU0Precon{F}, y::AbstractArray{F}) where {F<:Real} - ldiv!(x, f, y) +function ilu_apply!(x::AbstractArray{F}, f::ILU0Precon{F}, y::AbstractArray{F}, type::Symbol = :both, args...) where {F<:Real} + f! = ilu_f(type) + f!(x, f, y) end -function ilu_apply!(x, ilu::ILU0Precon, y) +# function ilu_apply!(x::AbstractArray{F}, f::CuSparseMatrix{F}, y::AbstractArray{F}, type::Symbol = :both) where {F<:Real} +# x .= y +# ix = 'O' +# sv2!('N', 'L', 'N', 1.0, f, x, ix) +# sv2!('N', 'U', 'U', 1.0, f, x, ix) +# end + +function ilu_apply!(x, ilu::ILU0Precon, y, type::Symbol = :both,args...) T = eltype(ilu.l_nzval) N = size(T, 1) T = eltype(T) diff --git a/src/linsolve/precond/utils.jl b/src/linsolve/precond/utils.jl index 1cf60e22..903ae1ae 100644 --- a/src/linsolve/precond/utils.jl +++ b/src/linsolve/precond/utils.jl @@ -1,11 +1,11 @@ - function update_preconditioner!(preconditioner::Nothing, arg...) # Do nothing. end -function update_preconditioner!(preconditioner, lsys, model, storage, recorder, executor) +function update_preconditioner!(preconditioner, lsys, context, model, storage, recorder, executor) J = jacobian(lsys) r = residual(lsys) - ctx = linear_system_context(model, lsys) + #ctx = linear_system_context(model, lsys) + ctx = context update_preconditioner!(preconditioner, J, r, ctx, executor) end @@ -20,11 +20,11 @@ end is_left_preconditioner(::JutulPreconditioner) = true is_right_preconditioner(::JutulPreconditioner) = false -function linear_operator(precond::JutulPreconditioner, float_t = Float64) +function linear_operator(precond::JutulPreconditioner, side::Symbol = :left, float_t = Float64, args...) n = operator_nrows(precond) function precond_apply!(res, x, α, β::T) where T if β == zero(T) - apply!(res, precond, x) + apply!(res, precond, x, float_t,args...) if α != one(T) lmul!(α, res) end @@ -36,6 +36,7 @@ function linear_operator(precond::JutulPreconditioner, float_t = Float64) return op end +#nead to be spesilized on type not all JutulPreconditioners has get_factor function apply!(x, p::JutulPreconditioner, y, arg...) factor = get_factorization(p) if is_left_preconditioner(p) diff --git a/src/linsolve/precond/various.jl b/src/linsolve/precond/various.jl index 775dfc8d..35215651 100644 --- a/src/linsolve/precond/various.jl +++ b/src/linsolve/precond/various.jl @@ -15,7 +15,7 @@ mutable struct LUPreconditioner <: JutulPreconditioner end end -function update_preconditioner!(lup::LUPreconditioner, A, b, context, executor) +function update_preconditioner!(lup::LUPreconditioner, A, b, context, executor,args...) if isnothing(lup.factor) lup.factor = lu(A) else @@ -29,19 +29,27 @@ function operator_nrows(lup::LUPreconditioner) return size(f.L, 1) end +#function operator_nrows(prec::TrivialPreconditioner) +# return prec.dim[1] +#end # LU factor as precond for wells? """ Trivial / identity preconditioner with size for use in subsystems. """ +function apply!(x,tp::TrivialPreconditioner,r, args...) + x = copy(r) +end + + # Trivial precond -function update_preconditioner!(tp::TrivialPreconditioner, lsys, model, storage, recorder, executor) +function update_preconditioner!(tp::TrivialPreconditioner, lsys, model, storage, recorder, executor, args...) A = jacobian(lsys) b = residual(lsys) tp.dim = size(A).*length(b[1]) end - -function linear_operator(id::TrivialPreconditioner, ::Symbol) +export linear_operator +function linear_operator(id::TrivialPreconditioner, ::Symbol, args...) return opEye(id.dim...) end diff --git a/src/linsolve/scalar_cpu.jl b/src/linsolve/scalar_cpu.jl index b3379275..2bbe864f 100644 --- a/src/linsolve/scalar_cpu.jl +++ b/src/linsolve/scalar_cpu.jl @@ -15,12 +15,12 @@ mutable struct LUSolver end end -function linear_solve!(sys, solver::LUSolver, arg...; kwarg...) +function linear_solve!(sys, solver::LUSolver, arg...;dx = sys.dx, r = sys.r, kwargs...) if length(sys.dx) > solver.max_size error("System too big for LU solver. You can increase max_size at your own peril.") end J = sys.jac - r = sys.r + #r = sys.r if !solver.reuse_memory F = lu(J) else @@ -32,7 +32,7 @@ function linear_solve!(sys, solver::LUSolver, arg...; kwarg...) F = solver.F end - sys.dx .= -(F\r) + dx .= -(F\r) @assert all(isfinite, sys.dx) "Linear solve resulted in non-finite values." return linear_solve_return() end diff --git a/src/models.jl b/src/models.jl index 31e5a164..1931574e 100644 --- a/src/models.jl +++ b/src/models.jl @@ -884,8 +884,9 @@ end function solve_and_update!(storage, model::JutulModel, dt = nothing; linear_solver = nothing, recorder = nothing, executor = default_executor(), kwarg...) lsys = storage.LinearizedSystem + context = model.context t_solve = @elapsed begin - @tic "linear solve" (ok, n, history) = linear_solve!(lsys, linear_solver, model, storage, dt, recorder, executor) + @tic "linear solve" (ok, n, history) = linear_solve!(lsys, linear_solver, context, model, storage, dt, recorder, executor) end t_update = @elapsed @tic "primary variables" update = update_primary_variables!(storage, model; kwarg...) return (t_solve, t_update, n, history, update) diff --git a/src/multimodel/model.jl b/src/multimodel/model.jl index 5261aa56..cb58ae01 100644 --- a/src/multimodel/model.jl +++ b/src/multimodel/model.jl @@ -130,6 +130,10 @@ function setup_storage!(storage, model::MultiModel; return storage end +mutable struct MutableWrapper + maps +end + function setup_multimodel_maps!(storage, model) groups = model.groups if isnothing(groups) @@ -138,6 +142,7 @@ function setup_multimodel_maps!(storage, model) offset_map = map(g -> get_submodel_offsets(model, g), unique(groups)) end storage[:multi_model_maps] = (offset_map = offset_map, ); + storage[:eq_maps] = MutableWrapper(nothing) end function setup_equations_and_primary_variable_views!(storage, model::MultiModel, lsys) From 4a20590c7b36051f81ce21c4ef0e9e8b8847bd28 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 3 Oct 2024 11:22:44 +0200 Subject: [PATCH 02/15] work in scaling --- src/equations.jl | 26 ++++++++++++++++++++++++++ src/models.jl | 10 ++++++++++ src/multimodel/model.jl | 17 +++++++++++++++++ 3 files changed, 53 insertions(+) diff --git a/src/equations.jl b/src/equations.jl index 67e6b47a..46a9b824 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -529,6 +529,32 @@ function update_equation!(eq_s, eq::JutulEquation, storage, model, dt) update_equation_for_entity!(cache, eq, state, state0, model, dt) end end + + apply_scaling_for_entity!(eq_s, eq, model) + +end + +function apply_scaling_for_entity!(eq_s, eq, model) + + scaling = get_scaling(model, eq) + @infiltrate + scaling = 1.0 + + for k in keys(eq_s) + if k == :numeric + continue + end + cache = eq_s[k] + cache.entries .*= 1/scaling + end + +end + +function apply_scaling_for_entity!(eq_s::AbstractMatrix, eq, model) + + scaling = get_scaling(model, eq) + eq_s .*= 1/scaling + end @inline prepare_equation_in_entity!(i, eq, eq_s, state, state0, model, dt) = nothing diff --git a/src/models.jl b/src/models.jl index 31e5a164..6a3e6820 100644 --- a/src/models.jl +++ b/src/models.jl @@ -721,6 +721,16 @@ function update_equations!(storage, equations_storage, equations, model, dt) end end + +function apply_scaling_equations!(storage, equations_storage, equations, model, dt) + for (key, eq) in pairs(equations) + @tic "$key" apply_scaling_equation!(equations_storage[key], eq, storage, model, dt) + end +end + + + + """ update_linearized_system!(storage, model::JutulModel; ) diff --git a/src/multimodel/model.jl b/src/multimodel/model.jl index 5261aa56..586f874f 100644 --- a/src/multimodel/model.jl +++ b/src/multimodel/model.jl @@ -569,6 +569,13 @@ function update_equations!(storage, model::MultiModel, dt; targets = submodels_s end end +function apply_scaling_equations!(storage, model::MultiModel, dt; targets = submodels_symbols(model)) + @tic "scaling equations" for k in targets + apply_scaling_equations!(storage[k], model[k], dt) + end +end + + function update_equations_and_apply_forces!(storage, model::MultiModel, dt, forces; time = NaN, kwarg...) # First update all equations @tic "equations" update_equations!(storage, model, dt; kwarg...) @@ -580,6 +587,10 @@ function update_equations_and_apply_forces!(storage, model::MultiModel, dt, forc @tic "boundary conditions" apply_boundary_conditions!(storage, model; kwarg...) # Apply forces to cross-terms @tic "crossterm forces" apply_forces_to_cross_terms!(storage, model, dt, forces; time = time, kwarg...) + # Apply scaling to equations + # @tic "crossterm forces" apply_scaling_equations!(storage, model, dt; kwarg...) + # Apply scaling to cross-terms + # @tic "crossterm forces" apply_scaling_cross_terms!(storage, model, dt; kwarg...) end function update_cross_terms!(storage, model::MultiModel, dt; targets = submodels_symbols(model), sources = submodels_symbols(model)) @@ -659,6 +670,10 @@ function update_cross_term_inner_target!(cache::GenericAutoDiffCache{<:Any, <:An update_cross_term_for_entity!(cache, ct, eq, state_t_local, state0_t_local, state_s, state0_s, model_t, model_s, dt) end +function get_scaling(model, eq) + return 1.0 +end + function update_cross_term_for_entity!(cache, ct, eq, state_t, state0_t, state_s, state0_s, model_t, model_s, dt) v = cache.entries vars = cache.variables @@ -672,6 +687,8 @@ function update_cross_term_for_entity!(cache, ct, eq, state_t, state0_t, state_s state_s = new_entity_index(state_s, var) state0_s = new_entity_index(state0_s, var) update_cross_term_in_entity!(v_i, i, state_t, state0_t, state_s, state0_s, model_t, model_s, ct, eq, dt, ldisc) + scaling = get_scaling(model_t, eq) + v_i = 1/scaling*v_i end end end From 6ab190458881b08e4038ab2a1fd0f43bb00169cf Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 3 Oct 2024 16:43:29 +0200 Subject: [PATCH 03/15] added scaling for cross terms --- src/equations.jl | 1 - src/multimodel/model.jl | 32 ++++++++++++++++++++++++++++++-- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/src/equations.jl b/src/equations.jl index 46a9b824..1dcae67e 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -537,7 +537,6 @@ end function apply_scaling_for_entity!(eq_s, eq, model) scaling = get_scaling(model, eq) - @infiltrate scaling = 1.0 for k in keys(eq_s) diff --git a/src/multimodel/model.jl b/src/multimodel/model.jl index 586f874f..ce9370f9 100644 --- a/src/multimodel/model.jl +++ b/src/multimodel/model.jl @@ -618,11 +618,41 @@ function update_cross_term!(ct_s, ct::CrossTerm, eq, storage_t, storage_s, model state0_s = storage_s.state0 if ct_s[:helper_mode] update_cross_term_helper_impl!(state_t, state0_t, state_s, state0_s, ct_s.target, ct_s.source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) + apply_scaling_cross_terms_helper!(ct_s.target, ct_s.source, eq, model_t) else update_cross_term_impl!(state_t, state0_t, state_s, state0_s, ct_s.target, ct_s.source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) + apply_scaling_cross_terms!(ct_s.target, ct_s.source, eq, model_t) end end + +function apply_scaling_cross_terms!(ct_s_target, + ct_s_source, + equation, + model) + + scaling = get_scaling(model, equation) + scaling = 1.0 + + for cache in values(ct_s_target) + cache.entries .*= 1/scaling + end + + for cache in values(ct_s_source) + cache.entries .*= 1/scaling + end + +end + +function apply_scaling_cross_terms_helper!(ct_s_target, + ct_s_source, + equation, + model) + + error("not yet implemented") +end + + function update_cross_term_impl!(state_t, state0_t, state_s, state0_s, ct_s_target, ct_s_source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) for i in 1:ct_s.N prepare_cross_term_in_entity!(i, state_t, state0_t, state_s, state0_s, model_t, model_s, ct, eq, dt) @@ -687,8 +717,6 @@ function update_cross_term_for_entity!(cache, ct, eq, state_t, state0_t, state_s state_s = new_entity_index(state_s, var) state0_s = new_entity_index(state0_s, var) update_cross_term_in_entity!(v_i, i, state_t, state0_t, state_s, state0_s, model_t, model_s, ct, eq, dt, ldisc) - scaling = get_scaling(model_t, eq) - v_i = 1/scaling*v_i end end end From 5e0df257eb50982a0eeec2cc1a306c1e9296f003 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 3 Oct 2024 17:00:41 +0200 Subject: [PATCH 04/15] working on scaling --- src/equations.jl | 1 - src/models.jl | 12 +++++++++--- src/multimodel/model.jl | 2 +- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/equations.jl b/src/equations.jl index 1dcae67e..8f57649a 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -530,7 +530,6 @@ function update_equation!(eq_s, eq::JutulEquation, storage, model, dt) end end - apply_scaling_for_entity!(eq_s, eq, model) end diff --git a/src/models.jl b/src/models.jl index 6a3e6820..2e978742 100644 --- a/src/models.jl +++ b/src/models.jl @@ -722,13 +722,19 @@ function update_equations!(storage, equations_storage, equations, model, dt) end -function apply_scaling_equations!(storage, equations_storage, equations, model, dt) - for (key, eq) in pairs(equations) - @tic "$key" apply_scaling_equation!(equations_storage[key], eq, storage, model, dt) +function apply_scaling_equations!(storage, model) + + for (key, equation) in pairs(model.equations) + @tic "$key" apply_scaling_equation!(storage.equations[key], equation, model) end + end +function apply_scaling_equation!(eq_s, equation::JutulEquation, model) + + apply_scaling_for_entity!(eq_s, equation, model) +end """ diff --git a/src/multimodel/model.jl b/src/multimodel/model.jl index ce9370f9..5117ef37 100644 --- a/src/multimodel/model.jl +++ b/src/multimodel/model.jl @@ -571,7 +571,7 @@ end function apply_scaling_equations!(storage, model::MultiModel, dt; targets = submodels_symbols(model)) @tic "scaling equations" for k in targets - apply_scaling_equations!(storage[k], model[k], dt) + apply_scaling_equations!(storage[k], model[k]) end end From 0383d696397dcc3ab631f53ae8f0d85c5a5137fd Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 3 Oct 2024 19:30:55 +0200 Subject: [PATCH 05/15] added scaling for conservation law storage --- src/conservation/conservation.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/conservation/conservation.jl b/src/conservation/conservation.jl index 449727e0..18f2e23d 100644 --- a/src/conservation/conservation.jl +++ b/src/conservation/conservation.jl @@ -496,6 +496,22 @@ function update_equation!(eq_s::ConservationLawTPFAStorage, law::ConservationLaw @tic "fluxes" update_half_face_flux!(eq_s, law, storage, model, dt) end +function apply_scaling_equation!(eq_s::ConservationLawTPFAStorage, equation::ConservationLaw, model) + + scaling = get_scaling(model, equation) + + fields = (:accumulation, :half_face_flux_cells, :half_face_flux_faces) + + for fd in fields + eq = getfield(eq_s, fd) + if !isnothing(eq) + eq.entries .*= 1/scaling + end + end + +end + + function update_half_face_flux!(eq_s::ConservationLawTPFAStorage, law::ConservationLaw, storage, model, dt) fd = law.flow_discretization state = storage.state From d82a5bb73738c983ecf54a3a6fdb6175f9868e79 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 3 Oct 2024 19:31:35 +0200 Subject: [PATCH 06/15] minor clean-up for scaling --- src/equations.jl | 3 --- src/multimodel/model.jl | 13 +++++-------- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/equations.jl b/src/equations.jl index 8f57649a..4a673250 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -529,14 +529,11 @@ function update_equation!(eq_s, eq::JutulEquation, storage, model, dt) update_equation_for_entity!(cache, eq, state, state0, model, dt) end end - - end function apply_scaling_for_entity!(eq_s, eq, model) scaling = get_scaling(model, eq) - scaling = 1.0 for k in keys(eq_s) if k == :numeric diff --git a/src/multimodel/model.jl b/src/multimodel/model.jl index 5117ef37..84ddefa2 100644 --- a/src/multimodel/model.jl +++ b/src/multimodel/model.jl @@ -588,9 +588,7 @@ function update_equations_and_apply_forces!(storage, model::MultiModel, dt, forc # Apply forces to cross-terms @tic "crossterm forces" apply_forces_to_cross_terms!(storage, model, dt, forces; time = time, kwarg...) # Apply scaling to equations - # @tic "crossterm forces" apply_scaling_equations!(storage, model, dt; kwarg...) - # Apply scaling to cross-terms - # @tic "crossterm forces" apply_scaling_cross_terms!(storage, model, dt; kwarg...) + @tic "scaling equations" apply_scaling_equations!(storage, model, dt; kwarg...) end function update_cross_terms!(storage, model::MultiModel, dt; targets = submodels_symbols(model), sources = submodels_symbols(model)) @@ -618,21 +616,20 @@ function update_cross_term!(ct_s, ct::CrossTerm, eq, storage_t, storage_s, model state0_s = storage_s.state0 if ct_s[:helper_mode] update_cross_term_helper_impl!(state_t, state0_t, state_s, state0_s, ct_s.target, ct_s.source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) - apply_scaling_cross_terms_helper!(ct_s.target, ct_s.source, eq, model_t) + apply_scaling_cross_term_helper!(ct_s.target, ct_s.source, eq, model_t) else update_cross_term_impl!(state_t, state0_t, state_s, state0_s, ct_s.target, ct_s.source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) - apply_scaling_cross_terms!(ct_s.target, ct_s.source, eq, model_t) + apply_scaling_cross_term!(ct_s.target, ct_s.source, eq, model_t) end end -function apply_scaling_cross_terms!(ct_s_target, +function apply_scaling_cross_term!(ct_s_target, ct_s_source, equation, model) scaling = get_scaling(model, equation) - scaling = 1.0 for cache in values(ct_s_target) cache.entries .*= 1/scaling @@ -644,7 +641,7 @@ function apply_scaling_cross_terms!(ct_s_target, end -function apply_scaling_cross_terms_helper!(ct_s_target, +function apply_scaling_cross_term_helper!(ct_s_target, ct_s_source, equation, model) From 32b803183d652f6135e9a454e8a3396527a689c0 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 7 Oct 2024 11:20:33 +0200 Subject: [PATCH 07/15] minor clean-up --- src/multimodel/model.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/multimodel/model.jl b/src/multimodel/model.jl index 84ddefa2..a6d529da 100644 --- a/src/multimodel/model.jl +++ b/src/multimodel/model.jl @@ -625,10 +625,10 @@ end function apply_scaling_cross_term!(ct_s_target, - ct_s_source, - equation, - model) - + ct_s_source, + equation, + model) + scaling = get_scaling(model, equation) for cache in values(ct_s_target) @@ -642,9 +642,9 @@ function apply_scaling_cross_term!(ct_s_target, end function apply_scaling_cross_term_helper!(ct_s_target, - ct_s_source, - equation, - model) + ct_s_source, + equation, + model) error("not yet implemented") end From 09d6fb3c9ac647fa3fc271d6b74a98322309efe7 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 13 Nov 2024 16:06:35 +0100 Subject: [PATCH 08/15] removed scaling of equations --- src/conservation/conservation.jl | 16 ------------ src/equations.jl | 21 ---------------- src/models.jl | 16 ------------ src/multimodel/model.jl | 42 -------------------------------- 4 files changed, 95 deletions(-) diff --git a/src/conservation/conservation.jl b/src/conservation/conservation.jl index 18f2e23d..449727e0 100644 --- a/src/conservation/conservation.jl +++ b/src/conservation/conservation.jl @@ -496,22 +496,6 @@ function update_equation!(eq_s::ConservationLawTPFAStorage, law::ConservationLaw @tic "fluxes" update_half_face_flux!(eq_s, law, storage, model, dt) end -function apply_scaling_equation!(eq_s::ConservationLawTPFAStorage, equation::ConservationLaw, model) - - scaling = get_scaling(model, equation) - - fields = (:accumulation, :half_face_flux_cells, :half_face_flux_faces) - - for fd in fields - eq = getfield(eq_s, fd) - if !isnothing(eq) - eq.entries .*= 1/scaling - end - end - -end - - function update_half_face_flux!(eq_s::ConservationLawTPFAStorage, law::ConservationLaw, storage, model, dt) fd = law.flow_discretization state = storage.state diff --git a/src/equations.jl b/src/equations.jl index 4a673250..67e6b47a 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -531,27 +531,6 @@ function update_equation!(eq_s, eq::JutulEquation, storage, model, dt) end end -function apply_scaling_for_entity!(eq_s, eq, model) - - scaling = get_scaling(model, eq) - - for k in keys(eq_s) - if k == :numeric - continue - end - cache = eq_s[k] - cache.entries .*= 1/scaling - end - -end - -function apply_scaling_for_entity!(eq_s::AbstractMatrix, eq, model) - - scaling = get_scaling(model, eq) - eq_s .*= 1/scaling - -end - @inline prepare_equation_in_entity!(i, eq, eq_s, state, state0, model, dt) = nothing function update_equation_for_entity!(cache, eq, state, state0, model, dt) diff --git a/src/models.jl b/src/models.jl index 46d4db58..1931574e 100644 --- a/src/models.jl +++ b/src/models.jl @@ -721,22 +721,6 @@ function update_equations!(storage, equations_storage, equations, model, dt) end end - -function apply_scaling_equations!(storage, model) - - for (key, equation) in pairs(model.equations) - @tic "$key" apply_scaling_equation!(storage.equations[key], equation, model) - end - -end - -function apply_scaling_equation!(eq_s, equation::JutulEquation, model) - - apply_scaling_for_entity!(eq_s, equation, model) - -end - - """ update_linearized_system!(storage, model::JutulModel; ) diff --git a/src/multimodel/model.jl b/src/multimodel/model.jl index 73537800..cb58ae01 100644 --- a/src/multimodel/model.jl +++ b/src/multimodel/model.jl @@ -574,13 +574,6 @@ function update_equations!(storage, model::MultiModel, dt; targets = submodels_s end end -function apply_scaling_equations!(storage, model::MultiModel, dt; targets = submodels_symbols(model)) - @tic "scaling equations" for k in targets - apply_scaling_equations!(storage[k], model[k]) - end -end - - function update_equations_and_apply_forces!(storage, model::MultiModel, dt, forces; time = NaN, kwarg...) # First update all equations @tic "equations" update_equations!(storage, model, dt; kwarg...) @@ -592,8 +585,6 @@ function update_equations_and_apply_forces!(storage, model::MultiModel, dt, forc @tic "boundary conditions" apply_boundary_conditions!(storage, model; kwarg...) # Apply forces to cross-terms @tic "crossterm forces" apply_forces_to_cross_terms!(storage, model, dt, forces; time = time, kwarg...) - # Apply scaling to equations - @tic "scaling equations" apply_scaling_equations!(storage, model, dt; kwarg...) end function update_cross_terms!(storage, model::MultiModel, dt; targets = submodels_symbols(model), sources = submodels_symbols(model)) @@ -621,40 +612,11 @@ function update_cross_term!(ct_s, ct::CrossTerm, eq, storage_t, storage_s, model state0_s = storage_s.state0 if ct_s[:helper_mode] update_cross_term_helper_impl!(state_t, state0_t, state_s, state0_s, ct_s.target, ct_s.source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) - apply_scaling_cross_term_helper!(ct_s.target, ct_s.source, eq, model_t) else update_cross_term_impl!(state_t, state0_t, state_s, state0_s, ct_s.target, ct_s.source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) - apply_scaling_cross_term!(ct_s.target, ct_s.source, eq, model_t) end end - -function apply_scaling_cross_term!(ct_s_target, - ct_s_source, - equation, - model) - - scaling = get_scaling(model, equation) - - for cache in values(ct_s_target) - cache.entries .*= 1/scaling - end - - for cache in values(ct_s_source) - cache.entries .*= 1/scaling - end - -end - -function apply_scaling_cross_term_helper!(ct_s_target, - ct_s_source, - equation, - model) - - error("not yet implemented") -end - - function update_cross_term_impl!(state_t, state0_t, state_s, state0_s, ct_s_target, ct_s_source, ct_s, ct::CrossTerm, eq, storage_t, storage_s, model_t, model_s, dt) for i in 1:ct_s.N prepare_cross_term_in_entity!(i, state_t, state0_t, state_s, state0_s, model_t, model_s, ct, eq, dt) @@ -702,10 +664,6 @@ function update_cross_term_inner_target!(cache::GenericAutoDiffCache{<:Any, <:An update_cross_term_for_entity!(cache, ct, eq, state_t_local, state0_t_local, state_s, state0_s, model_t, model_s, dt) end -function get_scaling(model, eq) - return 1.0 -end - function update_cross_term_for_entity!(cache, ct, eq, state_t, state0_t, state_s, state0_s, model_t, model_s, dt) v = cache.entries vars = cache.variables From 0f5e9aa27aa6abc605b21022a12c2bed763cfd84 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 20 Nov 2024 11:54:12 +0100 Subject: [PATCH 09/15] Clean up merge + remove "side" arg from ILU(0) --- src/discretization/finite-volume.jl | 1 - src/linsolve/krylov.jl | 2 +- src/linsolve/precond/amg.jl | 5 ----- src/linsolve/precond/ilu.jl | 30 +++++------------------------ 4 files changed, 6 insertions(+), 32 deletions(-) diff --git a/src/discretization/finite-volume.jl b/src/discretization/finite-volume.jl index bb697433..f76d15d6 100644 --- a/src/discretization/finite-volume.jl +++ b/src/discretization/finite-volume.jl @@ -116,7 +116,6 @@ function compute_half_face_trans(cell_centroids, face_centroids, face_normals, f T_hf[fpos] = T end end - #@assert minimum(T_hf) >= 0 return T_hf end diff --git a/src/linsolve/krylov.jl b/src/linsolve/krylov.jl index 447961c0..86ee33cd 100644 --- a/src/linsolve/krylov.jl +++ b/src/linsolve/krylov.jl @@ -56,7 +56,7 @@ function Base.show(io::IO, krylov::GenericKrylov) print(io, "Generic Krylov using $(krylov.solver) (ϵₐ=$atol, ϵ=$rtol) with prec = $(typeof(krylov.preconditioner))") end -function preconditioner(krylov::GenericKrylov, sys, context, model, storage, recorder, side, arg...) +function preconditioner(krylov::AbstractKrylov, sys, context, model, storage, recorder, side, arg...) M = krylov.preconditioner Ft = float_type(context) if isnothing(M) diff --git a/src/linsolve/precond/amg.jl b/src/linsolve/precond/amg.jl index f54b7735..9cf256af 100644 --- a/src/linsolve/precond/amg.jl +++ b/src/linsolve/precond/amg.jl @@ -21,11 +21,6 @@ end matrix_for_amg(A) = A matrix_for_amg(A::StaticSparsityMatrixCSR) = copy(A.At') -#function apply!(x, prec::AMGPreconditioner, y, float_t, sys, args...) -# apply!(x,prec,y) -#end - - function update_preconditioner!(amg::AMGPreconditioner{flavor}, A, b, context, executor) where flavor kw = amg.method_kwarg A_amg = matrix_for_amg(A) diff --git a/src/linsolve/precond/ilu.jl b/src/linsolve/precond/ilu.jl index 7419c8de..d96a9326 100644 --- a/src/linsolve/precond/ilu.jl +++ b/src/linsolve/precond/ilu.jl @@ -57,22 +57,10 @@ end function apply!(x, ilu::ILUZeroPreconditioner, y, type, arg...) factor = get_factorization(ilu) - mytype = :both - ilu_apply!(x, factor, y, mytype, arg...) + ilu_apply!(x, factor, y, arg...) end -function ilu_f(type::Symbol) - # Why must this be qualified? - if type == :left - f = forward_substitution! - elseif type == :right - f = backward_substitution! - else - f = ldiv! - end -end - -function ilu_apply!(x::AbstractArray{F}, f::AbstractILUFactorization, y::AbstractArray{F}, type::Symbol = :both, args...) where {F<:Real} +function ilu_apply!(x::AbstractArray{F}, f::AbstractILUFactorization, y::AbstractArray{F}, args...) where {F<:Real} T = eltype(f) if T == Float64 ldiv!(x, f, y) @@ -90,19 +78,11 @@ function ilu_apply!(x, f::AbstractILUFactorization, y) ldiv!(x, f, y) end -function ilu_apply!(x::AbstractArray{F}, f::ILU0Precon{F}, y::AbstractArray{F}, type::Symbol = :both, args...) where {F<:Real} - f! = ilu_f(type) - f!(x, f, y) +function ilu_apply!(x::AbstractArray{F}, f::ILU0Precon{F}, y::AbstractArray{F}, args...) where {F<:Real} + ldiv!(x, f, y) end -# function ilu_apply!(x::AbstractArray{F}, f::CuSparseMatrix{F}, y::AbstractArray{F}, type::Symbol = :both) where {F<:Real} -# x .= y -# ix = 'O' -# sv2!('N', 'L', 'N', 1.0, f, x, ix) -# sv2!('N', 'U', 'U', 1.0, f, x, ix) -# end - -function ilu_apply!(x, ilu::ILU0Precon, y, type::Symbol = :both,args...) +function ilu_apply!(x, ilu::ILU0Precon, y,args...) T = eltype(ilu.l_nzval) N = size(T, 1) T = eltype(T) From 5ed66b0d7b2aeaa971d5e8eaa1bd941d6c6bb1c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 20 Nov 2024 12:28:32 +0100 Subject: [PATCH 10/15] WIP: Remove some args... from linear solver, specify types --- src/linsolve/precond/utils.jl | 10 ++++------ src/linsolve/precond/various.jl | 4 ++-- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/linsolve/precond/utils.jl b/src/linsolve/precond/utils.jl index 903ae1ae..ad3db115 100644 --- a/src/linsolve/precond/utils.jl +++ b/src/linsolve/precond/utils.jl @@ -1,12 +1,10 @@ function update_preconditioner!(preconditioner::Nothing, arg...) # Do nothing. end -function update_preconditioner!(preconditioner, lsys, context, model, storage, recorder, executor) +function update_preconditioner!(preconditioner::JutulPreconditioner, lsys::JutulLinearSystem, context, model, storage, recorder, executor) J = jacobian(lsys) r = residual(lsys) - #ctx = linear_system_context(model, lsys) - ctx = context - update_preconditioner!(preconditioner, J, r, ctx, executor) + update_preconditioner!(preconditioner, J, r, context, executor) end function partial_update_preconditioner!(p, A, b, context, executor) @@ -20,11 +18,11 @@ end is_left_preconditioner(::JutulPreconditioner) = true is_right_preconditioner(::JutulPreconditioner) = false -function linear_operator(precond::JutulPreconditioner, side::Symbol = :left, float_t = Float64, args...) +function linear_operator(precond::JutulPreconditioner, side::Symbol = :left, float_t = Float64) n = operator_nrows(precond) function precond_apply!(res, x, α, β::T) where T if β == zero(T) - apply!(res, precond, x, float_t,args...) + apply!(res, precond, x, float_t) if α != one(T) lmul!(α, res) end diff --git a/src/linsolve/precond/various.jl b/src/linsolve/precond/various.jl index 35215651..d91232df 100644 --- a/src/linsolve/precond/various.jl +++ b/src/linsolve/precond/various.jl @@ -15,7 +15,7 @@ mutable struct LUPreconditioner <: JutulPreconditioner end end -function update_preconditioner!(lup::LUPreconditioner, A, b, context, executor,args...) +function update_preconditioner!(lup::LUPreconditioner, A, b, context, executor) if isnothing(lup.factor) lup.factor = lu(A) else @@ -43,7 +43,7 @@ end # Trivial precond -function update_preconditioner!(tp::TrivialPreconditioner, lsys, model, storage, recorder, executor, args...) +function update_preconditioner!(tp::TrivialPreconditioner, lsys, model, storage, recorder, executor) A = jacobian(lsys) b = residual(lsys) tp.dim = size(A).*length(b[1]) From ef43d1bded7dc95e4d3537f7fba6102231d66d98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 20 Nov 2024 14:05:49 +0100 Subject: [PATCH 11/15] Linear solver cleanup --- src/linsolve/krylov.jl | 7 +++---- src/linsolve/precond/ilu.jl | 10 +++++----- src/linsolve/precond/utils.jl | 16 +++++++--------- src/linsolve/precond/various.jl | 4 +++- 4 files changed, 18 insertions(+), 19 deletions(-) diff --git a/src/linsolve/krylov.jl b/src/linsolve/krylov.jl index 86ee33cd..95a65c60 100644 --- a/src/linsolve/krylov.jl +++ b/src/linsolve/krylov.jl @@ -56,13 +56,13 @@ function Base.show(io::IO, krylov::GenericKrylov) print(io, "Generic Krylov using $(krylov.solver) (ϵₐ=$atol, ϵ=$rtol) with prec = $(typeof(krylov.preconditioner))") end -function preconditioner(krylov::AbstractKrylov, sys, context, model, storage, recorder, side, arg...) +function preconditioner(krylov::AbstractKrylov, sys, context, model, storage, recorder) M = krylov.preconditioner Ft = float_type(context) if isnothing(M) op = I else - linop = linear_operator(M, side, Ft, sys, context, model, storage, recorder, arg...) + linop = linear_operator(M, Ft, sys, context, model, storage, recorder) op = PrecondWrapper(linop) end return op @@ -95,8 +95,7 @@ function linear_solve!(sys::LSystem, t_prep = @elapsed @tic "prepare" prepare_linear_solve!(sys) op = linear_operator(sys) t_prec = @elapsed @tic "precond" update_preconditioner!(prec, sys, context, model, storage, recorder, executor) - prec_op = preconditioner(krylov, sys, context, model, storage, recorder, :left) - # R = preconditioner(krylov, sys, model, storage, recorder, :right, Ft) + prec_op = preconditioner(krylov, sys, context, model, storage, recorder) v = Int64(cfg.verbose) max_it = cfg.max_iterations min_it = cfg.min_iterations diff --git a/src/linsolve/precond/ilu.jl b/src/linsolve/precond/ilu.jl index d96a9326..b244a88b 100644 --- a/src/linsolve/precond/ilu.jl +++ b/src/linsolve/precond/ilu.jl @@ -55,12 +55,12 @@ function update_preconditioner!(ilu::ILUZeroPreconditioner, A::StaticSparsityMat end end -function apply!(x, ilu::ILUZeroPreconditioner, y, type, arg...) +function apply!(x, ilu::ILUZeroPreconditioner, y) factor = get_factorization(ilu) - ilu_apply!(x, factor, y, arg...) + ilu_apply!(x, factor, y) end -function ilu_apply!(x::AbstractArray{F}, f::AbstractILUFactorization, y::AbstractArray{F}, args...) where {F<:Real} +function ilu_apply!(x::AbstractArray{F}, f::AbstractILUFactorization, y::AbstractArray{F}) where {F<:Real} T = eltype(f) if T == Float64 ldiv!(x, f, y) @@ -78,11 +78,11 @@ function ilu_apply!(x, f::AbstractILUFactorization, y) ldiv!(x, f, y) end -function ilu_apply!(x::AbstractArray{F}, f::ILU0Precon{F}, y::AbstractArray{F}, args...) where {F<:Real} +function ilu_apply!(x::AbstractArray{F}, f::ILU0Precon{F}, y::AbstractArray{F}) where {F<:Real} ldiv!(x, f, y) end -function ilu_apply!(x, ilu::ILU0Precon, y,args...) +function ilu_apply!(x, ilu::ILU0Precon, y) T = eltype(ilu.l_nzval) N = size(T, 1) T = eltype(T) diff --git a/src/linsolve/precond/utils.jl b/src/linsolve/precond/utils.jl index ad3db115..4bb569f9 100644 --- a/src/linsolve/precond/utils.jl +++ b/src/linsolve/precond/utils.jl @@ -18,7 +18,11 @@ end is_left_preconditioner(::JutulPreconditioner) = true is_right_preconditioner(::JutulPreconditioner) = false -function linear_operator(precond::JutulPreconditioner, side::Symbol = :left, float_t = Float64) +function linear_operator(precond::JutulPreconditioner) + return linear_operator(precond, Float64, nothing, nothing, nothing, nothing, nothing) +end + +function linear_operator(precond::JutulPreconditioner, float_t, sys, context, model, storage, recorder) n = operator_nrows(precond) function precond_apply!(res, x, α, β::T) where T if β == zero(T) @@ -35,13 +39,7 @@ function linear_operator(precond::JutulPreconditioner, side::Symbol = :left, flo end #nead to be spesilized on type not all JutulPreconditioners has get_factor -function apply!(x, p::JutulPreconditioner, y, arg...) +function apply!(x, p::JutulPreconditioner, y) factor = get_factorization(p) - if is_left_preconditioner(p) - ldiv!(x, factor, y) - elseif is_right_preconditioner(p) - error("Not supported.") - else - error("Neither left or right preconditioner?") - end + ldiv!(x, factor, y) end diff --git a/src/linsolve/precond/various.jl b/src/linsolve/precond/various.jl index d91232df..bb965c59 100644 --- a/src/linsolve/precond/various.jl +++ b/src/linsolve/precond/various.jl @@ -48,8 +48,10 @@ function update_preconditioner!(tp::TrivialPreconditioner, lsys, model, storage, b = residual(lsys) tp.dim = size(A).*length(b[1]) end + export linear_operator -function linear_operator(id::TrivialPreconditioner, ::Symbol, args...) + +function linear_operator(id::TrivialPreconditioner, args...) return opEye(id.dim...) end From 66611f76735de594e4a4bcae4dab420289a99165 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 20 Nov 2024 14:14:39 +0100 Subject: [PATCH 12/15] Remove float_t from apply --- src/linsolve/precond/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linsolve/precond/utils.jl b/src/linsolve/precond/utils.jl index 4bb569f9..c2c999d4 100644 --- a/src/linsolve/precond/utils.jl +++ b/src/linsolve/precond/utils.jl @@ -26,7 +26,7 @@ function linear_operator(precond::JutulPreconditioner, float_t, sys, context, mo n = operator_nrows(precond) function precond_apply!(res, x, α, β::T) where T if β == zero(T) - apply!(res, precond, x, float_t) + apply!(res, precond, x) if α != one(T) lmul!(α, res) end From 89df330af1553ad9ad31ba91b33180c10d822ba3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 20 Nov 2024 14:52:17 +0100 Subject: [PATCH 13/15] Fix PArray to new interface --- ext/JutulPartitionedArraysExt/linalg.jl | 3 ++- src/linsolve/precond/ilu.jl | 10 +++++++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/ext/JutulPartitionedArraysExt/linalg.jl b/ext/JutulPartitionedArraysExt/linalg.jl index 39dccffc..827c4b6a 100644 --- a/ext/JutulPartitionedArraysExt/linalg.jl +++ b/ext/JutulPartitionedArraysExt/linalg.jl @@ -67,7 +67,8 @@ function Jutul.parray_update_preconditioners!(sim, preconditioner_base, precondi model = Jutul.get_simulator_model(sim) storage = Jutul.get_simulator_storage(sim) sys = storage.LinearizedSystem - Jutul.update_preconditioner!(prec, sys, model, storage, recorder, sim.executor) + ctx = model.context + Jutul.update_preconditioner!(prec, sys, ctx, model, storage, recorder, sim.executor) prec end return (preconditioner_base, preconditioners) diff --git a/src/linsolve/precond/ilu.jl b/src/linsolve/precond/ilu.jl index b244a88b..9757f823 100644 --- a/src/linsolve/precond/ilu.jl +++ b/src/linsolve/precond/ilu.jl @@ -34,17 +34,21 @@ function update_preconditioner!(ilu::ILUZeroPreconditioner, A, b, context, execu end end -function update_preconditioner!(ilu::ILUZeroPreconditioner, A::StaticSparsityMatrixCSR, b, context::ParallelCSRContext, executor) +function update_preconditioner!(ilu::ILUZeroPreconditioner, A::StaticSparsityMatrixCSR, b, context, executor) if isnothing(ilu.factor) mb = A.minbatch max_t = max(size(A, 1) ÷ mb, 1) - nt = min(nthreads(context), max_t) + nt = min(A.nthreads, max_t) if nt == 1 @debug "Setting up serial ILU(0)-CSR" F = ilu0_csr(A) else @debug "Setting up parallel ILU(0)-CSR with $(nthreads(td)) threads" - part = context.partitioner + if context isa ParallelCSRContext + part = context.partitioner + else + part = MetisPartitioner() + end lookup = generate_lookup(part, A, nt) F = ilu0_csr(A, lookup) end From 318c1acc6c247ed87718c3f7f26f66907d489bfc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 20 Nov 2024 14:59:31 +0100 Subject: [PATCH 14/15] Krylov: Add backwards compatible interface --- src/linsolve/krylov.jl | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/linsolve/krylov.jl b/src/linsolve/krylov.jl index 95a65c60..4fca8ead 100644 --- a/src/linsolve/krylov.jl +++ b/src/linsolve/krylov.jl @@ -68,15 +68,10 @@ function preconditioner(krylov::AbstractKrylov, sys, context, model, storage, re return op end -# export update! -# function update_preconditioner!(prec, sys, model, storage, recorder) -# update!(prec, sys, model, storage, recorder) -# end - function linear_solve!(sys::LSystem, krylov::GenericKrylov, - context, - model, + context::JutulContext, + model::JutulModel, storage = nothing, dt = nothing, recorder = ProgressRecorder(), @@ -186,6 +181,16 @@ function linear_solve!(sys::LSystem, return linear_solve_return(solved, n, stats, precond = t_prec_op, precond_count = t_prec_count, prepare = t_prec + t_prep) end +function linear_solve!( + sys::LSystem, + krylov::GenericKrylov, + model::JutulModel, + arg...; + kwarg... + ) + return linear_solve!(sys, krylov, model.context, model, arg...; kwarg...) +end + function krylov_scale_system!(sys, krylov::GenericKrylov, dt) sys, krylov.storage_scaling = apply_scaling_to_linearized_system!(sys, krylov.storage_scaling, krylov.scaling, dt) return sys From 209459604e1ec75a27f928de146abc7beb40616d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 20 Nov 2024 15:23:03 +0100 Subject: [PATCH 15/15] Always have context in MultiModel --- src/core_types/core_types.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/core_types/core_types.jl b/src/core_types/core_types.jl index 9e5d5996..6173147b 100644 --- a/src/core_types/core_types.jl +++ b/src/core_types/core_types.jl @@ -965,7 +965,17 @@ struct MultiModel{label, T, CT, G, C, GL} <: AbstractMultiModel{label} group_lookup::GL end -function MultiModel(models, label::Union{Nothing, Symbol} = nothing; cross_terms = Vector{CrossTermPair}(), groups = nothing, context = nothing, reduction = missing, specialize = false, specialize_ad = false) +function MultiModel(models, label::Union{Nothing, Symbol} = nothing; + cross_terms = Vector{CrossTermPair}(), + groups = nothing, + context = nothing, + reduction = missing, + specialize = false, + specialize_ad = false + ) + if isnothing(context) + context = models[first(keys(models))].context + end group_lookup = Dict{Symbol, Int}() if isnothing(groups) num_groups = 1