Merge pull request #83 from sintefmath/dev

Hysteretic parameters and MPI memory usage improvements

Author: moyner

Project.toml

@@ -1,7 +1,7 @@
 name = "Jutul"
 uuid = "2b460a1a-8a2b-45b2-b125-b5c536396eb9"
 authors = ["Olav Møyner <[email protected]>"]
-version = "0.2.34"
+version = "0.2.35"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -78,9 +78,9 @@ LayeredLayouts = "0.2.5"
 LinearOperators = "2.3.2"
 LoopVectorization = "0.12.115"
 MAT = "0.10.3"
-Makie = "0.20.0"
+Makie = "0.21.0"
 MappedArrays = "0.4.1"
-Meshes = "0.28 - 0.41, 0.46"
+Meshes = "0.28 - 0.41"
 Metis = "1.1.0"
 NetworkLayout = "0.4.5"
 OrderedCollections = "1.4.1"

ext/JutulMakieExt/interactive_3d.jl

@@ -392,6 +392,7 @@ function plot_interactive_impl(grid, states;
         "batlowK",
         "berlin",
         "brg",
+        "commercial",
         "gnuplot",
         "gray1",
         "hawaii",
@@ -553,7 +554,6 @@ function plot_interactive_impl(grid, states;
         tri = primitives.triangulation
         scat = Makie.mesh!(ax, pts, tri; color = ys,
                                         colorrange = lims,
-                                        size = 60,
                                         backlight = 1,
                                         colormap = cmap,
                                         transparency = transparency,
@@ -675,7 +675,11 @@ function unpack(buffer, x::AbstractVector, ix)
 end
 
 function generate_colormap(colormap_name, alphamap_name, base_alpha, low, high)
-    cmap = to_colormap(colormap_name)
+    if colormap_name == :commercial
+        cmap = commercial_colormap()
+    else
+        cmap = to_colormap(colormap_name)
+    end
     n = length(cmap)
     if alphamap_name != :no_alpha_map
         if alphamap_name == :linear
@@ -878,3 +882,27 @@ function Jutul.plotting_check_interactive(; warn = true)
     end
     return true
 end
+
+function commercial_colormap()
+    blue =   (0, 0, 1)
+    cyan =   (0, 1, 1)
+    green =  (0, 1, 0)
+    yellow = (1, 1, 0)
+    red =    (1, 0, 0)
+
+    function simple_interp(F_0, F_1, x)
+        v = F_0 .+ (F_1 .- F_0).*x
+        return Makie.RGB(v...)
+    end
+    cmap = Vector{typeof(Makie.RGB(0, 0, 0))}()
+    nsteps = [30, 25, 25, 20]
+    colors = (blue, cyan, green, yellow, red)
+    for (i, nstep) in enumerate(nsteps)
+        c1 = colors[i]
+        c2 = colors[i+1]
+        for dx in range(0.0, 1.0, nstep)
+            push!(cmap, simple_interp(c1, c2, dx))
+        end
+    end
+    return cmap
+end

ext/JutulMeshesExt/JutulMeshesExt.jl

@@ -8,7 +8,7 @@ module JutulMeshesExt
         # Global neighborship
         N = zeros(Int64, 2, nf)
         # Normals and centroids
-        normals = Vector{Vec3}(undef, nf)
+        normals = Vector{Meshes.Vec3}(undef, nf)
         cell_centroids = centroid.(grid)
         face_centroids = Vector{Point3}(undef, nf)
         # Measures

ext/JutulPartitionedArraysExt/io.jl

@@ -1,14 +1,21 @@
-function consolidate_distributed_results_on_disk!(pth, np, steps; cleanup = true)
+function Jutul.consolidate_distributed_results_on_disk!(pth, np, steps; cleanup = true, verbose = false)
     @assert isdir(pth)
     partitions = []
+    function print_msg(x)
+        if verbose
+            jutul_message("IO", x)
+        end
+    end
     for rank in 1:np
+        print_msg("Reading partition for process $rank/$np")
         proc_pth = rank_folder(pth, rank)
         part_path = joinpath(proc_pth, "partition.jld2")
         p = load(part_path)
         push!(partitions, p)
     end
     allpaths = []
     for step in steps
+        print_msg("Reading result for state $step/$(length(steps))")
         step::Int
         states, reports, paths = read_step(pth, step, np)
         state, report = consolidate_distributed_results(states, reports, partitions)
@@ -17,6 +24,7 @@ function consolidate_distributed_results_on_disk!(pth, np, steps; cleanup = true
         GC.gc()
     end
     if cleanup
+        print_msg("Cleaning up files...")
         # Delete the parallel results after consolidation.
         for paths in allpaths
             for path in paths

ext/JutulPartitionedArraysExt/overloads.jl

@@ -283,11 +283,12 @@ function Jutul.retrieve_output!(sim::PArraySimulator, states, reports, config, n
     if config[:consolidate_results]
         MPI.Barrier(comm)
         if np > 1 && is_main && pth isa String
-            consolidate_distributed_results_on_disk!(pth, np, 1:n, cleanup = config[:delete_on_consolidate])
+            Jutul.consolidate_distributed_results_on_disk!(pth, np, 1:n, cleanup = config[:delete_on_consolidate], verbose = config[:info_level] > 0)
         end
         MPI.Barrier(comm)
     end
-    Jutul.retrieve_output!(states, reports, config, n)
+    states, reports = Jutul.retrieve_output!(states, reports, config, n, read_states = config[:output_states] && is_main)
+    return (states, reports)
 end
 
 function Jutul.simulator_reports_per_step(psim::PArraySimulator)

src/conservation/conservation.jl

@@ -366,7 +366,6 @@ function update_lsys_sources_theaded!(nz, r, acc, src, rv_src, nz_src, cp, conte
                 pos = get_jacobian_pos(acc, cell, e, d, cp)
                 @inbounds nz[pos] += v.partials[d]
             end
-            # @info "Entry $e, $cell ($nc x $np): $(v.value)"
         end
     end
 end

src/core_types/contexts/interface.jl

@@ -1,6 +1,19 @@
 export minbatch
-minbatch(::Any) = 1000
-nthreads(::Any) = Threads.nthreads()
+
+function minbatch(x::Any)
+    nt = nthreads(x)
+    if nt == 1
+        nb = typemax(Int)
+    else
+        nb = 1000
+    end
+    return nb
+end
+
+function nthreads(::Any)
+    Threads.nthreads()
+end
+
 minbatch(x, n) = max(n ÷ nthreads(x), minbatch(x))
 
 function jacobian_eltype(context, layout, block_size)

src/ext/partitionedarrays_ext.jl

@@ -28,6 +28,10 @@ function mpi_scalar_allreduce
 
 end
 
+function consolidate_distributed_results_on_disk!
+
+end
+
 abstract type PArrayBackend <: JutulBackend end
 
 struct DebugPArrayBackend <: PArrayBackend end

src/interpolation.jl

@@ -284,7 +284,7 @@ function get_dependencies(var::UnaryTabulatedVariable, model)
     return [var.x_symbol]
 end
 
-function update_secondary_variable!(V, var::UnaryTabulatedVariable, model, state)
+function update_secondary_variable!(V, var::UnaryTabulatedVariable, model, state, ix = entity_eachindex(V))
     update_unary_tabulated!(V, var, model, state[var.x_symbol], entity_eachindex(V))
 end
 

src/models.jl

@@ -930,7 +930,10 @@ function update_before_step!(storage, model, dt, forces; kwarg...)
 end
 
 function update_before_step!(storage, ::Any, model, dt, forces; time = NaN, recorder = ProgressRecorder(), update_explicit = true)
-    # Do nothing
+    state = storage.state
+    for (k, prm) in pairs(storage.variable_definitions.parameters)
+        update_parameter_before_step!(state[k], prm, storage, model, dt, forces)
+    end
 end
 
 function update_after_step!(storage, model, dt, forces; kwarg...)
@@ -963,10 +966,20 @@ function update_after_step!(storage, model, dt, forces; kwarg...)
     return report
 end
 
+"""
+    update_parameter_before_step!(prm_val, prm, storage, model, dt, forces)
+
+Update parameters before time-step. Used for hysteretic parameters.
+"""
+function update_parameter_before_step!(prm_val, prm, storage, model, dt, forces)
+    # Do nothing by default.
+    return prm_val
+end
+
 function variable_change_report(X::AbstractArray, X0::AbstractArray{T}, pvar) where T<:Real
     max_dv = max_v = sum_dv = sum_v = zero(T)
     @inbounds @simd for i in eachindex(X)
-        x = value(X[i])
+        x = value(X[i])::T
         dx = x - value(X0[i])
 
         dx_abs = abs(dx)