From a941fc0d187d245aad19489071820ac493ab6915 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 29 Nov 2024 12:05:23 +0100 Subject: [PATCH 01/15] Add definition to map_ij_to_active --- src/dd/finite_volume_map.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/dd/finite_volume_map.jl b/src/dd/finite_volume_map.jl index cca0daf2..2edd8f38 100644 --- a/src/dd/finite_volume_map.jl +++ b/src/dd/finite_volume_map.jl @@ -76,6 +76,11 @@ function map_to_active(V, domain, m::FiniteVolumeGlobalMap, ::Cells) # return filter(i -> m.cell_is_boundary[i], V) end +function map_ij_to_active(I, J, domain, m::FiniteVolumeGlobalMap, e) + # At the moment only cells can have inactive status + return (I, J) +end + function map_ij_to_active(I, J, domain, m::FiniteVolumeGlobalMap, ::Cells) n = length(I) @assert n == length(J) From e63a49d559102a599ebe8c0c7c7b029ef789734a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 8 Dec 2024 15:42:45 +0100 Subject: [PATCH 02/15] Fix bug in update direction --- src/variables/utils.jl | 60 ++++++++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 26 deletions(-) diff --git a/src/variables/utils.jl b/src/variables/utils.jl index 53a11967..b1e6cf3b 100644 --- a/src/variables/utils.jl +++ b/src/variables/utils.jl @@ -1,4 +1,4 @@ -const MINIMUM_SAT_RELAX = 1e-6 +const MINIMUM_SAT_RELAX = 1e-3 """ Number of entities (e.g. Cells, Faces) a variable is defined on. By default, each primary variable exists on all cells of a discretized domain @@ -370,6 +370,8 @@ function update_primary_variable!(state, p::FractionVariables, state_symbol, mod unit_sum_update!(s, p, model, dx, w) end +unit_update_preserve_direction(::FractionVariables) = true + function unit_sum_update!(s, p, model, dx, w, entity = Cells()) nf, nu = value_dim(model, p) abs_max = absolute_increment_limit(p) @@ -380,7 +382,7 @@ function unit_sum_update!(s, p, model, dx, w, entity = Cells()) if nf == 2 unit_update_pairs!(s, dx, active_cells, minval, maxval, abs_max, w) else - if true + if unit_update_preserve_direction(p) # Preserve direction unit_update_direction!(s, dx, nf, nu, active_cells, minval, maxval, abs_max, w) else @@ -418,28 +420,32 @@ function unit_update_direction_local!(s, active_ix, full_cell, dx, nf, nu, minva bad_update = w <= MINIMUM_SAT_RELAX if bad_update - w = w0 - end - @inbounds for i in 1:(nf-1) - s[i, full_cell] += w*dx[i, active_ix] - end - @inbounds s[nf, full_cell] += w*dlast0 - if bad_update - # Dampening is tiny, update and renormalize instead - tot = 0.0 - @inbounds for i in 1:nf - s_i = s[i, full_cell] - sat = clamp(value(s_i), minval, maxval) - tot += sat - s_i = replace_value(s_i, sat) - s[i, full_cell] = s_i + # It was not possible to find a reasonable dampening factor. + # We instead go to the magnitude-preserving version. + unit_update_magnitude_local!(s, active_ix, full_cell, dx, nf, nu, minval, maxval, abs_max) + else + @inbounds for i in 1:(nf-1) + s[i, full_cell] += w*dx[i, active_ix] end - @inbounds for i = 1:nf - s_i = s[i, full_cell] - s_i = replace_value(s_i, value(s_i)/tot) - s[i, full_cell] = s_i + @inbounds s[nf, full_cell] += w*dlast0 + if bad_update + # Dampening is tiny, update and renormalize instead + tot = 0.0 + @inbounds for i in 1:nf + s_i = s[i, full_cell] + sat = clamp(value(s_i), minval, maxval) + tot += sat + s_i = replace_value(s_i, sat) + s[i, full_cell] = s_i + end + @inbounds for i = 1:nf + s_i = s[i, full_cell] + s_i = replace_value(s_i, value(s_i)/tot) + s[i, full_cell] = s_i + end end end + return s end function unit_update_pairs!(s, dx, active_cells, minval, maxval, abs_max, w) @@ -456,18 +462,19 @@ function unit_update_pairs!(s, dx, active_cells, minval, maxval, abs_max, w) end function unit_update_magnitude!(s, dx, nf, nu, active_cells, minval, maxval, abs_max) - for cell in active_cells - unit_update_magnitude_local!(s, cell, dx, nf, nu, minval, maxval, abs_max) + for active_ix in eachindex(active_cells) + cell = active_cells[active_ix] + unit_update_magnitude_local!(s, active_ix, cell, dx, nf, nu, minval, maxval, abs_max) end end -function unit_update_magnitude_local!(s, cell, dx, nf, nu, minval, maxval, abs_max) +function unit_update_magnitude_local!(s, ix, cell, dx, nf, nu, minval, maxval, abs_max) # First pass: Find the relaxation factors that keep all fractions in [0, 1] # and obeying the maximum change targets - dlast0 = 0 + dlast0 = zero(eltype(dx)) @inbounds for i = 1:(nf-1) v = value(s[i, cell]) - dv = dx[cell + (i-1)*nu] + dv = dx[i, ix] dv = choose_increment(v, dv, abs_max, nothing, minval, maxval) s[i, cell] += dv dlast0 -= dv @@ -476,6 +483,7 @@ function unit_update_magnitude_local!(s, cell, dx, nf, nu, minval, maxval, abs_m dlast = choose_increment(value(s[nf, cell]), dlast0, abs_max, nothing, minval, maxval) s[nf, cell] += dlast if dlast != dlast0 + # Need to renormalize since the last value was not within bounds. t = 0.0 for i = 1:nf @inbounds t += s[i, cell] From 042d822e049bacdef59431f5eaf543b52cafa94b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 8 Dec 2024 15:50:09 +0100 Subject: [PATCH 03/15] Avoid messing up AD status in update_magnitude --- src/variables/utils.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/variables/utils.jl b/src/variables/utils.jl index b1e6cf3b..0735bbcc 100644 --- a/src/variables/utils.jl +++ b/src/variables/utils.jl @@ -486,7 +486,8 @@ function unit_update_magnitude_local!(s, ix, cell, dx, nf, nu, minval, maxval, a # Need to renormalize since the last value was not within bounds. t = 0.0 for i = 1:nf - @inbounds t += s[i, cell] + # Note: Careful to handle AD values correctly here. + @inbounds t += value(s[i, cell]) end for i = 1:nf @inbounds s[i, cell] /= t From 1f4f3c5887fe04b672dbb7d778c42fa534c82683 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sun, 8 Dec 2024 15:56:01 +0100 Subject: [PATCH 04/15] Fix another AD bug in old update code --- src/variables/utils.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/variables/utils.jl b/src/variables/utils.jl index 0735bbcc..217a16b3 100644 --- a/src/variables/utils.jl +++ b/src/variables/utils.jl @@ -489,8 +489,9 @@ function unit_update_magnitude_local!(s, ix, cell, dx, nf, nu, minval, maxval, a # Note: Careful to handle AD values correctly here. @inbounds t += value(s[i, cell]) end - for i = 1:nf - @inbounds s[i, cell] /= t + @inbounds for i = 1:nf + s_i = s[i, cell] + s[i, cell] = replace_value(s_i, clamp(value(s_i), minval, maxval)/t) end end end From 96a152c22a8ba7539c53384a23c281bb2a45f85e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 11 Dec 2024 20:42:54 +0100 Subject: [PATCH 05/15] Add option to successful_reports --- src/timesteps.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/timesteps.jl b/src/timesteps.jl index cd2b2b3c..1dc43fde 100644 --- a/src/timesteps.jl +++ b/src/timesteps.jl @@ -152,9 +152,10 @@ end successful_reports(old_reports, current_reports, step_index, n = 1) Get the `n` last successful solve reports from all previous reports -(old_reports) and the current ministep set. +(old_reports) and the current ministep set. You can optionally provide a +function that replaces the default definition of `success=r->r[:success]`. """ -function successful_reports(old_reports, current_reports, step_index, n = 1) +function successful_reports(old_reports, current_reports, step_index, n = 1; success = r -> r[:success]) out = similar(current_reports, 0) if isfinite(n) sizehint!(out, n) @@ -169,7 +170,7 @@ function successful_reports(old_reports, current_reports, step_index, n = 1) end for r in Iterators.reverse(reports) - if !ismissing(r) && r[:success] + if !ismissing(r) && success(r) push!(out, r) if length(out) >= n return out @@ -187,12 +188,12 @@ Get last n successful reports starting at the end of `step` and reversing backwards until `n` values have been found. `n` can be set to `Inf` to produce all successful reports. """ -function successful_reports(reports, current_reports = missing; step = length(reports)+1, n = 1) +function successful_reports(reports, current_reports = missing; step = length(reports)+1, n = 1, kwarg...) if ismissing(current_reports) step = clamp(step, 1, length(reports)) current_reports = reports[step][:ministeps] end - return successful_reports(reports, current_reports, step, n) + return successful_reports(reports, current_reports, step, n; kwarg...) end """ From 2f2645f0335635b722f4768bea79a44fef8635a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Fri, 13 Dec 2024 11:43:53 +0100 Subject: [PATCH 06/15] Adjust default WENO epsilon --- src/WENO/WENO.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/WENO/WENO.jl b/src/WENO/WENO.jl index 883b73be..26a13ab3 100644 --- a/src/WENO/WENO.jl +++ b/src/WENO/WENO.jl @@ -24,7 +24,7 @@ module WENO r::WENOHalfFaceDiscretization{D, N, R}; do_clamp = false, threshold = 0.0, - epsilon = 1e-10 + epsilon = 1e-7 ) where {D, N, R} @assert D == N-1 return new{D, N, R}(l, r, do_clamp, threshold, epsilon) From 1159b898f5766043a61fcce3304d810300f83c67 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Sat, 14 Dec 2024 14:52:21 +0100 Subject: [PATCH 07/15] Fix offset bug in tol_factor_final_iteration --- src/simulator/simulator.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulator/simulator.jl b/src/simulator/simulator.jl index 91bd2e4a..a21fc5a2 100644 --- a/src/simulator/simulator.jl +++ b/src/simulator/simulator.jl @@ -406,7 +406,7 @@ function perform_step_check_convergence_impl!(report, prev_report, storage, mode converged = false e = NaN t_conv = @elapsed begin - if iteration == config[:max_nonlinear_iterations] + if iteration == config[:max_nonlinear_iterations]+1 tf = config[:tol_factor_final_iteration] else tf = 1 From d91c9f393286b49612e60d2237735c443cf2e705 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 16 Dec 2024 12:43:59 +0100 Subject: [PATCH 08/15] Speed up subset plot_mesh --- ext/JutulMakieExt/mesh_plots.jl | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/ext/JutulMakieExt/mesh_plots.jl b/ext/JutulMakieExt/mesh_plots.jl index 7ac51df9..25ece817 100644 --- a/ext/JutulMakieExt/mesh_plots.jl +++ b/ext/JutulMakieExt/mesh_plots.jl @@ -27,14 +27,28 @@ function Jutul.plot_mesh_impl!(ax, m; has_face_filter = !isnothing(faces) has_bface_filter = !isnothing(boundaryfaces) if has_cell_filter || has_face_filter || has_bface_filter + keep_cells = Dict{Int, Bool}() + keep_faces = Dict{Int, Bool}() + keep_bf = Dict{Int, Bool}() + if eltype(cells) == Bool @assert length(cells) == number_of_cells(m) cells = findall(cells) end + if has_cell_filter + for c in cells + keep_cells[c] = true + end + end if eltype(faces) == Bool @assert length(faces) == number_of_faces(m) faces = findall(faces) end + if has_face_filter + for f in faces + keep_faces[f] = true + end + end if eltype(boundaryfaces) == Bool @assert length(boundaryfaces) == number_of_boundary_faces(m) boundaryfaces = findall(boundaryfaces) @@ -43,6 +57,9 @@ function Jutul.plot_mesh_impl!(ax, m; nf = number_of_faces(m) boundaryfaces = deepcopy(boundaryfaces) boundaryfaces .+= nf + for f in boundaryfaces + keep_bf[f] = true + end end ntri = size(tri, 1) keep = fill(false, ntri) @@ -53,13 +70,13 @@ function Jutul.plot_mesh_impl!(ax, m; tri_tmp = tri[i, 1] keep_this = true if has_cell_filter - keep_this = keep_this && cell_ix[tri_tmp] in cells + keep_this = keep_this && haskey(keep_cells, cell_ix[tri_tmp]) end if has_face_filter - keep_this = keep_this && face_ix[tri_tmp] in faces + keep_this = keep_this && haskey(keep_faces, face_ix[tri_tmp]) end if has_bface_filter - keep_this = keep_this && face_ix[tri_tmp] in boundaryfaces + keep_this = keep_this && haskey(keep_bf, face_ix[tri_tmp]) end keep[i] = keep_this end From bbbf5e8fa1d977459f156ac0606df7bfbd8e581f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 17 Dec 2024 22:20:50 +0100 Subject: [PATCH 09/15] Add a new timestepper --- src/timesteps.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/timesteps.jl b/src/timesteps.jl index 1dc43fde..33991d2f 100644 --- a/src/timesteps.jl +++ b/src/timesteps.jl @@ -148,6 +148,30 @@ function pick_next_timestep(sel::VariableChangeTimestepSelector, sim, config, dt return linear_timestep_selection(x, x0, x1, dt0, dt1) end +""" + sel = LimitByFailedTimestepSelector(num = 10, factor = 0.9) + +Limit the timestep by the shortest of `num` failed timesteps, reducing the +timestep by `factor` multiplied by the shortest failed timestep. If no +time-steps failed during the last `num` steps, the timestep is not changed. +""" +Base.@kwdef struct LimitByFailedTimestepSelector <: AbstractTimestepSelector + num::Int = 10 + factor::Float64 = 0.9 +end + +function pick_next_timestep(sel::LimitByFailedTimestepSelector, sim, config, dt_prev, dT, forces, reports, current_reports, step_index, new_step) + R = successful_reports(reports, current_reports, step_index, sel.num) + dt = dT + for rep in R + if !rep[:success] + @info "Limiting step." + dt = min(dt, rep[:dt]*sel.factor) + end + end + return dt +end + """ successful_reports(old_reports, current_reports, step_index, n = 1) From dbe3efe8cf59282fc48bb669ce89a8f32fa7eeb2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Tue, 17 Dec 2024 22:27:35 +0100 Subject: [PATCH 10/15] Make new timestepper actually work --- src/timesteps.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/timesteps.jl b/src/timesteps.jl index 33991d2f..e3f75713 100644 --- a/src/timesteps.jl +++ b/src/timesteps.jl @@ -161,11 +161,10 @@ Base.@kwdef struct LimitByFailedTimestepSelector <: AbstractTimestepSelector end function pick_next_timestep(sel::LimitByFailedTimestepSelector, sim, config, dt_prev, dT, forces, reports, current_reports, step_index, new_step) - R = successful_reports(reports, current_reports, step_index, sel.num) + R = successful_reports(reports, current_reports, step_index, sel.num, success = r -> true) dt = dT for rep in R if !rep[:success] - @info "Limiting step." dt = min(dt, rep[:dt]*sel.factor) end end From c9668c12016a1a4e403496747dab350bcd315fd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Wed, 18 Dec 2024 11:30:31 +0100 Subject: [PATCH 11/15] Improvements for new timestepper --- src/simulator/config.jl | 2 +- src/timesteps.jl | 6 +++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/simulator/config.jl b/src/simulator/config.jl index e2c9b603..98addda9 100644 --- a/src/simulator/config.jl +++ b/src/simulator/config.jl @@ -46,7 +46,7 @@ Negative values disable output. The interpretation of this number is subject to # IO options add_option!(cfg, :output_path, nothing, "Path to write output. If nothing, output is not written to disk.", types = Union{String, Nothing}) - add_option!(cfg, :in_memory_reports, 5, "Limit for number of reports kept in memory if output_path is provided.", types = Int) + add_option!(cfg, :in_memory_reports, 10, "Limit for number of reports kept in memory if output_path is provided.", types = Int) add_option!(cfg, :report_level, 0, "Level of information stored in reports when written to disk.", types = Int) add_option!(cfg, :output_substates, false, "Store substates (between report steps) as field on each state.", types = Bool) diff --git a/src/timesteps.jl b/src/timesteps.jl index e3f75713..6d1c751d 100644 --- a/src/timesteps.jl +++ b/src/timesteps.jl @@ -161,7 +161,7 @@ Base.@kwdef struct LimitByFailedTimestepSelector <: AbstractTimestepSelector end function pick_next_timestep(sel::LimitByFailedTimestepSelector, sim, config, dt_prev, dT, forces, reports, current_reports, step_index, new_step) - R = successful_reports(reports, current_reports, step_index, sel.num, success = r -> true) + R = successful_reports(reports, current_reports, n = sel.num, success = r -> true) dt = dT for rep in R if !rep[:success] @@ -189,6 +189,10 @@ function successful_reports(old_reports, current_reports, step_index, n = 1; suc if step == step_index reports = current_reports else + report = old_reports[step] + if ismissing(report) + continue + end reports = old_reports[step][:ministeps] end From 39d58a46331ada02e534eefbc5bc5f03cd97d845 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Thu, 19 Dec 2024 20:56:20 +0100 Subject: [PATCH 12/15] Update the performance plots --- ext/JutulMakieExt/performance.jl | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/ext/JutulMakieExt/performance.jl b/ext/JutulMakieExt/performance.jl index 07efef45..5f851f2a 100644 --- a/ext/JutulMakieExt/performance.jl +++ b/ext/JutulMakieExt/performance.jl @@ -131,14 +131,22 @@ function Jutul.plot_cumulative_solve!(f, allreports, dt = nothing, names = nothi c = colors[mod(i-1, n_rep)+1] data_i = get_data(r_rep[i]) push!(alldata, data_i) - lines!(ax, t[i], data_i, - label = names[i], - linewidth = linewidth, - color = c, - linestyle = get_linestyle(i); - kwarg...) + lstyle = get_linestyle(i) + skip_line = ismissing(lstyle) + if !skip_line + lines!(ax, t[i], data_i, + label = names[i], + linewidth = linewidth, + color = c, + linestyle = lstyle; + kwarg...) + end if scatter_points - scatter!(ax, t[i], data_i, color = c) + if skip_line + scatter!(ax, t[i], data_i, color = c, label = names[i]) + else + scatter!(ax, t[i], data_i, color = c) + end end end if length(r_rep) > 1 && !names_missing && legend From 337a19375cd381f8b42b2852adfdbb00824e5f10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 6 Jan 2025 15:07:06 +0100 Subject: [PATCH 13/15] Add convenience constructor for IndirectionMap --- src/core_types/core_types.jl | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/core_types/core_types.jl b/src/core_types/core_types.jl index 6173147b..cf4c3177 100644 --- a/src/core_types/core_types.jl +++ b/src/core_types/core_types.jl @@ -1087,6 +1087,24 @@ struct IndirectionMap{V} end end +""" + imap = IndirectionMap(vec_of_vec::Vector{V}) where V + +Create indirection map for a variable length dense vector that is represented as +a Vector of Vectors. +""" +function IndirectionMap(vec_of_vec::Vector{Vector{T}}) where T + vals = T[] + pos = Int[1] + for subvec in vec_of_vec + for v in subvec + push!(vals, v) + end + push!(pos, pos[end] + length(subvec)) + end + return IndirectionMap(vals, pos) +end + struct IndexRenumerator{T} indices::Dict{T, Int} end From bca3f0056d06e75dcfcdd1db906d0aa4cf06ce5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 6 Jan 2025 18:27:54 +0100 Subject: [PATCH 14/15] Gmsh extension (#111) * Gmsh extension * Import routine before testing --- Project.toml | 3 + ext/JutulGmshExt/JutulGmshExt.jl | 14 ++ ext/JutulGmshExt/interface.jl | 66 +++++++++ ext/JutulGmshExt/utils.jl | 225 +++++++++++++++++++++++++++++++ src/ext/extensions.jl | 1 + src/ext/gmsh_ext.jl | 12 ++ src/utils.jl | 21 +++ test/utils.jl | 17 ++- 8 files changed, 358 insertions(+), 1 deletion(-) create mode 100644 ext/JutulGmshExt/JutulGmshExt.jl create mode 100644 ext/JutulGmshExt/interface.jl create mode 100644 ext/JutulGmshExt/utils.jl create mode 100644 src/ext/gmsh_ext.jl diff --git a/Project.toml b/Project.toml index a2970e48..0a1e219a 100644 --- a/Project.toml +++ b/Project.toml @@ -41,6 +41,7 @@ Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" +Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" KaHyPar = "2a6221f6-aa48-11e9-3542-2d9e0ef01880" LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e" @@ -54,6 +55,7 @@ PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" JutulAMGCLWrapExt = "AMGCLWrap" JutulGLMakieExt = "GLMakie" JutulGraphMakieExt = ["GraphMakie", "NetworkLayout", "LayeredLayouts", "Makie"] +JutulGmshExt = "Gmsh" JutulHYPREExt = "HYPRE" JutulKaHyParExt = "KaHyPar" JutulMakieExt = "Makie" @@ -101,6 +103,7 @@ julia = "1.7" CPUSummary = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2" +Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771" KaHyPar = "2a6221f6-aa48-11e9-3542-2d9e0ef01880" LayeredLayouts = "f4a74d36-062a-4d48-97cd-1356bad1de4e" diff --git a/ext/JutulGmshExt/JutulGmshExt.jl b/ext/JutulGmshExt/JutulGmshExt.jl new file mode 100644 index 00000000..e37ce8b9 --- /dev/null +++ b/ext/JutulGmshExt/JutulGmshExt.jl @@ -0,0 +1,14 @@ +module JutulGmshExt + using Jutul + import Jutul: UnstructuredMesh, IndirectionMap, check_equal_perm + using Gmsh: Gmsh, gmsh + + using StaticArrays + using OrderedCollections + + const QUAD_T = SVector{4, Int} + const TRI_T = SVector{3, Int} + + include("utils.jl") + include("interface.jl") +end diff --git a/ext/JutulGmshExt/interface.jl b/ext/JutulGmshExt/interface.jl new file mode 100644 index 00000000..c6a90fb6 --- /dev/null +++ b/ext/JutulGmshExt/interface.jl @@ -0,0 +1,66 @@ + +function Jutul.mesh_from_gmsh(pth; verbose = false, kwarg...) + Gmsh.initialize() + ext = pth |> splitext |> last + gmsh.open(pth) + if lowercase(ext) == ".geo" + gmsh.model.mesh.generate() + end + + dim = gmsh.model.getDimension() + dim == 3 || error("Only 3D models are supported") + + gmsh.model.mesh.removeDuplicateNodes() + node_tags, pts, = gmsh.model.mesh.getNodes() + node_remap = Dict{UInt64, Int}() + for (i, tag) in enumerate(node_tags) + tag::UInt64 + node_remap[tag] = i + end + remaps = ( + nodes = node_remap, + faces = Dict{UInt64, Int}(), + cells = Dict{UInt64, Int}() + ) + pts = reshape(pts, Int(dim), :) + pts_s = collect(vec(reinterpret(SVector{3, Float64}, pts))) + + @assert size(pts, 2) == length(node_tags) + faces_to_nodes = parse_faces(remaps, verbose = verbose) + face_lookup = generate_face_lookup(faces_to_nodes) + + cells_to_faces = parse_cells(remaps, faces_to_nodes, face_lookup, verbose = verbose) + # We are done with Gmsh.jl at this point, finalize it. + Gmsh.finalize() + + neighbors = build_neighbors(cells_to_faces, faces_to_nodes, face_lookup) + + # Make both of these in case we have rogue faces that are not connected to any cell. + bnd_faces = Int[] + int_faces = Int[] + n_bad = 0 + for i in axes(neighbors, 2) + l, r = neighbors[:, i] + l_bnd = l == 0 + r_bnd = r == 0 + + if l_bnd || r_bnd + if l_bnd && r_bnd + n_bad += 1 + else + push!(bnd_faces, i) + end + else + push!(int_faces, i) + end + end + bnd_neighbors, bnd_faces_to_nodes, bnd_cells_to_faces = split_boundary(neighbors, faces_to_nodes, cells_to_faces, bnd_faces, boundary = true) + int_neighbors, int_faces_to_nodes, int_cells_to_faces = split_boundary(neighbors, faces_to_nodes, cells_to_faces, int_faces, boundary = false) + + c2f = IndirectionMap(int_cells_to_faces) + c2b = IndirectionMap(bnd_cells_to_faces) + f2n = IndirectionMap(int_faces_to_nodes) + b2n = IndirectionMap(bnd_faces_to_nodes) + print_message("Mesh parsed successfully:\n $(length(c2f)) cells\n $(length(f2n)) internal faces\n $(length(b2n)) boundary faces\n $(length(pts_s)) nodes", verbose) + return UnstructuredMesh(c2f, c2b, f2n, b2n, pts_s, int_neighbors, bnd_neighbors; kwarg...) +end diff --git a/ext/JutulGmshExt/utils.jl b/ext/JutulGmshExt/utils.jl new file mode 100644 index 00000000..70ef07ac --- /dev/null +++ b/ext/JutulGmshExt/utils.jl @@ -0,0 +1,225 @@ + +function add_next!(faces, remap, tags, numpts, offset) + vals = Int[] + for j in 1:numpts + push!(vals, remap[tags[offset + j]]) + end + push!(faces, vals) +end + +function parse_faces(remaps; verbose = false) + node_remap = remaps.nodes + face_remap = remaps.faces + faces = Vector{Int}[] + for (dim, tag) in gmsh.model.getEntities() + if dim != 2 + continue + end + type = gmsh.model.getType(dim, tag) + name = gmsh.model.getEntityName(dim, tag) + elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag) + type = gmsh.model.getType(dim, tag) + ename = gmsh.model.getEntityName(dim, tag) + for (etypes, etags, enodetags) in zip(elemTypes, elemTags, elemNodeTags) + name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(etypes) + if name == "Quadrilateral 4" + numpts = 4 + elseif name == "Triangle 3" + numpts = 3 + else + error("Unsupported element type $name for faces.") + end + @assert length(enodetags) == numpts*length(etags) + print_message("Faces: Processing $(length(etags)) tags of type $name", verbose) + for (i, etag) in enumerate(etags) + offset = (i-1)*numpts + add_next!(faces, node_remap, enodetags, numpts, offset) + face_remap[etag] = length(faces) + end + print_message("Added $(length(etags)) faces of type $name with $(length(unique(enodetags))) unique nodes", verbose) + end + end + return faces +end + +function get_cell_decomposition(name) + if name == "Hexahedron 8" + tris = Tuple{}() + quads = ( + QUAD_T(0, 4, 7, 3), + QUAD_T(1, 2, 6, 5), + QUAD_T(0, 1, 5, 4), + QUAD_T(2, 3, 7, 6), + QUAD_T(0, 3, 2, 1), + QUAD_T(4, 5, 6, 7) + ) + numpts = 8 + elseif name == "Tetrahedron 4" + tris = ( + TRI_T(0, 1, 3), + TRI_T(0, 2, 1), + TRI_T(0, 3, 2), + TRI_T(1, 2, 3) + ) + quads = Tuple{}() + numpts = 4 + elseif name == "Pyramid 5" + # TODO: Not really tested. + tris = ( + TRI_T(0, 1, 4), + TRI_T(0, 4, 3), + TRI_T(3, 4, 2), + TRI_T(1, 2, 4) + ) + quads = (QUAD_T(0, 3, 2, 1),) + numpts = 4 + elseif name == "Prism 6" + # TODO: Not really tested. + tris = ( + TRI_T(0, 2, 1), + TRI_T(3, 4, 5) + ) + quads = ( + QUAD_T(0, 1, 4, 3), + QUAD_T(0, 3, 5, 2), + QUAD_T(1, 2, 5, 4) + ) + numpts = 6 + else + error("Unsupported element type $name for cells.") + end + return (tris, quads, numpts) +end + +function print_message(msg, verbose) + if verbose + println(msg) + end +end + +function parse_cells(remaps, faces, face_lookup; verbose = false) + node_remap = remaps.nodes + face_remap = remaps.faces + cell_remap = remaps.cells + cells = Vector{Tuple{Int, Int}}[] + for (dim, tag) in gmsh.model.getEntities() + if dim != 3 + continue + end + type = gmsh.model.getType(dim, tag) + name = gmsh.model.getEntityName(dim, tag) + # Get the mesh elements for the entity (dim, tag): + elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(dim, tag) + # * Type and name of the entity: + type = gmsh.model.getType(dim, tag) + ename = gmsh.model.getEntityName(dim, tag) + for (etypes, etags, enodetags) in zip(elemTypes, elemTags, elemNodeTags) + name, dim, order, numv, parv, _ = gmsh.model.mesh.getElementProperties(etypes) + tris, quads, numpts = get_cell_decomposition(name) + print_message("Cells: Processing $(length(etags)) tags of type $name", verbose) + @assert length(enodetags) == numpts*length(etags) + nadded = 0 + for (i, etag) in enumerate(etags) + offset = (i-1)*numpts + pt_range = (offset+1):(offset+numpts) + @assert length(pt_range) == numpts + pts = map(i -> node_remap[enodetags[i]], pt_range) + cell = Tuple{Int, Int}[] + cellno = length(cells) + for face_t in (tris, quads) + for (fno, face) in enumerate(face_t) + face_pts = map(i -> pts[i+1], face) + face_pts_sorted = sort(face_pts) + faceno = get(face_lookup, face_pts_sorted, 0) + if faceno == 0 + nadded += 1 + push!(faces, face_pts) + faceno = length(faces) + face_lookup[face_pts_sorted] = faceno + sgn = 1 + else + sgn = check_equal_perm(face_pts, faces[faceno]) ? 1 : 2 + end + push!(cell, (faceno, sgn)) + end + end + push!(cells, cell) + end + print_message("Added $(length(etags)) new cells of type $name and $nadded new faces.", verbose) + end + end + return cells +end + +function build_neighbors(cells, faces, face_lookup) + neighbors = zeros(Int, 2, length(faces)) + for (cellno, cell_to_faces) in enumerate(cells) + for (face_index, lr) in cell_to_faces + face_pts = faces[face_index] + face_pts_sorted = sort(face_pts) + faceno = face_lookup[face_pts_sorted] + face_ref = faces[faceno] + oldn = neighbors[lr, faceno] + oldn == 0 || error("Cannot overwrite face neighbor for cell $cellno - was already defined as $oldn for index $lr: $(neighbors[:, faceno])") + neighbors[lr, faceno] = cellno + end + end + return neighbors +end + +function generate_face_lookup(faces) + face_lookup = Dict{Union{QUAD_T, TRI_T}, Int}() + + for (i, face) in enumerate(faces) + n = length(face) + if n == 3 + ft = sort(TRI_T(face[1], face[2], face[3])) + elseif n == 4 + ft = sort(QUAD_T(face[1], face[2], face[3], face[4])) + else + error("Unsupported face type with $n nodes, only 3 (for tri) and 4 (for quad) are known.") + end + @assert issorted(ft) + face_lookup[ft] = i + end + return face_lookup +end + +function split_boundary(neighbors, faces_to_nodes, cells_to_faces, active_ix::Vector{Int}; boundary::Bool) + remap = OrderedDict{Int, Int}() + for (i, ix) in enumerate(active_ix) + remap[ix] = i + end + # is_active = [false for _ in eachindex(faces_to_nodes)] + # is_active[active_ix] .= true + # Make renumbering here. + if boundary + new_neighbors = Int[] + for ix in active_ix + l, r = neighbors[:, ix] + @assert l == 0 || r == 0 + push!(new_neighbors, max(l, r)) + end + else + new_neighbors = Tuple{Int, Int}[] + for ix in active_ix + l, r = neighbors[:, ix] + @assert l != 0 && r != 0 + push!(new_neighbors, (l, r)) + end + end + new_faces_to_nodes = map(copy, faces_to_nodes[active_ix]) + # Handle cells -> current type of faces + new_cells_to_faces = Vector{Int}[] + for cell_to_faces in cells_to_faces + new_cell = Int[] + for (face, sgn) in cell_to_faces + if haskey(remap, face) + push!(new_cell, remap[face]) + end + end + push!(new_cells_to_faces, new_cell) + end + + return (new_neighbors, new_faces_to_nodes, new_cells_to_faces) +end diff --git a/src/ext/extensions.jl b/src/ext/extensions.jl index 0e38a0ce..e2583949 100644 --- a/src/ext/extensions.jl +++ b/src/ext/extensions.jl @@ -5,3 +5,4 @@ include("graphmakie_ext.jl") include("hypre_ext.jl") include("partitionedarrays_ext.jl") include("kahypar_ext.jl") +include("gmsh_ext.jl") diff --git a/src/ext/gmsh_ext.jl b/src/ext/gmsh_ext.jl new file mode 100644 index 00000000..17147d86 --- /dev/null +++ b/src/ext/gmsh_ext.jl @@ -0,0 +1,12 @@ +export mesh_from_gmsh + +""" + G = mesh_from_gmsh(pth; verbose = false) + +Parse a Gmsh file and return a Jutul `UnstructuredMesh` (in 3D only). Requires +the Gmsh.jl package to be loaded. Please note that Gmsh is GPL licensed unless +you have obtained another type of license from the authors. +""" +function mesh_from_gmsh + +end diff --git a/src/utils.jl b/src/utils.jl index 056137d7..856761a5 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1188,3 +1188,24 @@ function get_mat_writable_file_from_report(report; config = missing) end return out end + +function check_equal_perm(a, b) + N = length(a) + M = length(b) + N == M || throw(ArgumentError("Lengths of a and b do not match ($N != $M)")) + to_cyclic = i -> mod(i - 1, N) + 1 + start = a[1] + offset = findfirst(isequal(start), b)-1 + if isnothing(offset) + is_equal = false + else + is_equal = true + for i in 2:N + if a[i] != b[to_cyclic(i + offset)] + is_equal = false + break + end + end + end + return is_equal +end diff --git a/test/utils.jl b/test/utils.jl index 00435812..4dc6cce4 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -185,7 +185,6 @@ end xfine = range(-1.0, 5.0, nf) yfine = range(-1.0, 6.0, nf) using ForwardDiff, Test - F(X) = I(first(X), last(X)) function F_num(X) ϵ = 1e-6 @@ -351,3 +350,19 @@ end @test ismissing(get_mesh_entity_tag(g, Cells(), :group2, :tag3, throw = false)) end end + +import Jutul: check_equal_perm +@testset "check_equal_perm" begin + @test check_equal_perm(SVector(1, 2, 3), SVector(1, 2, 3)) + @test check_equal_perm(SVector(1, 2, 3), SVector(2, 3, 1)) + @test check_equal_perm(SVector(1, 2, 3), SVector(3, 1, 2)) + @test !check_equal_perm(SVector(1, 2, 3), SVector(3, 2, 1)) + @test !check_equal_perm(SVector(1, 2, 3), SVector(3, 2, 1)) + @test !check_equal_perm(SVector(1, 2, 3), SVector(2, 1, 3)) + @test !check_equal_perm(SVector(1, 2, 3), SVector(1, 3, 2)) + @test !check_equal_perm(SVector(1, 2, 3), SVector(1, 2, 5)) + @test check_equal_perm(SVector(1, 2, 3, 4), SVector(1, 2, 3, 4)) + @test check_equal_perm(SVector(1, 2, 3, 4), SVector(2, 3, 4, 1)) + @test check_equal_perm(SVector(1, 2, 3, 4), SVector(3, 4, 1, 2)) + @test check_equal_perm(SVector(1, 2, 3, 4), SVector(4, 1, 2, 3)) +end From 0b8bef6cd89821de6de2fc27b72da6bde58f01b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olav=20M=C3=B8yner?= Date: Mon, 6 Jan 2025 19:50:40 +0100 Subject: [PATCH 15/15] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0a1e219a..eec3de0c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Jutul" uuid = "2b460a1a-8a2b-45b2-b125-b5c536396eb9" authors = ["Olav Møyner "] -version = "0.3.2" +version = "0.3.3" [deps] AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"