Merge branch 'main' into compathelper/new_version/2023-08-20-01-13-20…

…-980-01032934808

Project.toml

@@ -11,12 +11,12 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
-IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MKL = "33e6dc65-8f57-5167-99aa-e5a354878fb2"
+Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -25,6 +25,9 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 ArgParse = "1"
+StaticArrays = "1"
+TimerOutputs = "0.5"
+
 julia = "1"
 
 [extras]

src/Erebus.jl

@@ -27,7 +27,6 @@ else
         using AppleAccelerate
     else
         using MKL
-        BLAS.set_num_threads(4)
     end
 end
 
@@ -3037,19 +3036,20 @@ $(SIGNATURES)
 # Details
 
     - ps: Instance of pardiso solver
-    - iparms: dictionary of iparm parameters
+    - iparms_dict: dictionary of iparm parameters
 
 # Returns
 
     - nothing
 """
-function initialize_pardiso!(pardiso_solver, iparms)
+function initialize_pardiso!(pardiso_solver, iparms_dict)
     set_msglvl!(pardiso_solver, Pardiso.MESSAGE_LEVEL_OFF)
     set_matrixtype!(pardiso_solver, Pardiso.REAL_NONSYM)
-    set_nprocs!(pardiso_solver, parse(Int32, ENV["OMP_NUM_THREADS"]))
-    for (i, v) in iparms
+    set_nprocs!(pardiso_solver, cache_kwargs.nprocs)
+    for (i, v) in iparms_dict
         set_iparm!(pardiso_solver, i+1, v)
     end
+    set_phase!(pardiso_solver, Pardiso.ANALYSIS)
 end
 
 """
@@ -3567,7 +3567,8 @@ function assemble_hydromechanical_lse!(
     end # @inbounds 
     flush!(L) # finalize CSC matrix
 # end # @timeit to "assemble_hydromechanical_lse()"
-    return L
+    # return L
+    return L.cscmatrix
 end # function assemble_hydromechanical_lse!
 
 """
@@ -6338,11 +6339,11 @@ function simulation_loop(output_path)
     # thermal solver
     RT, ST = setup_thermal_lse()
     # gravitational solver
-    RP, SP= setup_gravitational_lse()
+    RP, SP = setup_gravitational_lse()
     # Pardiso MKL solver
     if use_pardiso
-        pardiso_solver = MKLPardisoSolver()
-        initialize_pardiso!(pardiso_solver, iparms)
+        pardiso_solver = Pardiso.MKLPardisoSolver()
+        initialize_pardiso!(pardiso_solver, iparms_dict)
     # else
     #     F = lu(fdrand(Nx1*Ny1*6, 1, 1, matrixtype=ExtendableSparseMatrix))
     end
@@ -6704,34 +6705,48 @@ function simulation_loop(output_path)
                 @info "starting hydro-mechanical solver $titer-$iplast"
 #     @timeit to "solve hydromechanical system" begin
                 if use_pardiso
-                    # S = solve(pardiso_solver, L.cscmatrix, R)
                     set_phase!(
                         pardiso_solver, Pardiso.ANALYSIS_NUM_FACT_SOLVE_REFINE)
                     pardiso(
                         pardiso_solver,
                         S,
-                        get_matrix(pardiso_solver, L.cscmatrix, :N),
+                        get_matrix(pardiso_solver, L, :N),
                         R
                     )
                     set_phase!(pardiso_solver, Pardiso.RELEASE_ALL)
-                    pardiso(pardiso_solver, S, L.cscmatrix, R)
+                    pardiso(pardiso_solver, S, L, R)
                 else
-                    # S = L \ R
                     if timestep == 1 && iplast == 1
-                        s1 = rand(Ny1*Nx1*6)
                         hydromech_prob = LinearProblem(
-                            L, R; u0 = s1, alias_A = true, alias_b = true)
+                            L,
+                            R;
+                            u0 = rand(Ny1*Nx1*6),
+                            alias_A = true,
+                            alias_b = true
+                            )
+                        # if use_pardiso
+                        #     hydromech_cache = init(
+                        #         hydromech_prob,
+                        #         MKLPardisoIterate();
+                        #         # MKLPardisoFactorize();
+                        #         cache_kwargs...
+                        #         )
+                        #     @info hydromech_cache.cacheval.iparm 
+                        # else
                         hydromech_cache = init(
                             hydromech_prob,
-                            UMFPACKFactorization(reuse_symbolic=true);
+                            UMFPACKFactorization();
                             cache_kwargs...
-                        )
+                            )
+                        # end
                     else
-                        hydromech_cache = LinearSolve.set_A(hydromech_sol.cache, L)
-                        hydromech_cache = LinearSolve.set_b(hydromech_sol.cache, R)
+                        hydromech_cache = hydromech_sol.cache
+                        hydromech_cache.b = R
+                        hydromech_cache.A = L
                     end
-                    hydromech_sol = solve(hydromech_cache)
-                    # @info "SOL CHECK" norm(S-hydromech_sol.u)
+                    Pl = ILU0Preconditioner(L)
+                    hydromech_sol = LinearSolve.solve(hydromech_cache, Pl = Pl)
+                    S = hydromech_sol.u
                 end
 #     end # @timeit to "solve hydromechanical system"
                 @info "finished hydro-mechanical solver $titer-$iplast"

src/constants.jl

@@ -363,18 +363,15 @@ const dphimax = 100.01
 # starting timestep
 const start_step = 1
 # maximum number of timesteps to run
-# const n_steps = 10 
-const n_steps = 30_000 
+const n_steps = 10 
+# const n_steps = 30_000 
 # random number generator seed
 const seed = 42
-# iterative solver keyword arguments
-const cache_kwargs = (
-    ; verbose = true, abstol = 1e-8, reltol = 1e-8, maxiter = 30)
 # using MKL Pardiso solver
 const use_pardiso = false
 # MKL Pardiso solver IPARM control parameters -> ∇ATTN: zero-indexed as in docs:
 # https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html
-const iparms = Dict([
+const iparms_dict = Dict([
     (0, 1), # in: do not use default
     (1, 2), # in: nested dissection from METIS
     (2, 0), # in: reserved, set to zero
@@ -384,7 +381,7 @@ const iparms = Dict([
     (6, 0), # out: number of iterative refinement steps performed
     (7, 20), # in: maximum number of iterative refinement steps
     (8, 0), # in: tolerance level for relative residual, only with iparm[23]=1
-    (9, 13), # in: pivoting perturbation
+    (9, 12), # in: pivoting perturbation
     (10, 1), # in: scaling vectors
     (11, 1), # in: solve AX=B (no transpose, conjugate transpose): CSC matrix
     (12, 1), # in: improved accuracy using (non-)symmetric weighted matching
@@ -439,4 +436,16 @@ const iparms = Dict([
     (61, 0), # reserved, set to zero 
     (62, 0), # out: size of the minimum OOC memory for factorization and sol 
     (63, 0), # reserved, set to zero 
-])
+])
+# MKL Pardiso solver IPARM control parameters formatted for LinearSolve.jl
+const iparms = collect(
+    (key + 1, iparms_dict[key]) for key in sort!(collect(keys(iparms_dict))))
+# LinearSolve.jl solver keyword arguments
+const cache_kwargs = (;
+    nprocs = 16,
+    verbose = true,
+    abstol = 1e-8,
+    reltol = 1e-8,
+    maxiter = 30,
+    iparm = iparms,
+    )

src/test_constants.jl

@@ -364,14 +364,11 @@ const n_steps = 10
 # const n_steps = 30_000 
 # random number generator seed
 const seed = 42
-# iterative solver keyword arguments
-const cache_kwargs = (
-    ; verbose = true, abstol = 1e-8, reltol = 1e-8, maxiter = 30)
 # using MKL Pardiso solver
 const use_pardiso = false
 # MKL Pardiso solver IPARM control parameters -> ∇ATTN: zero-indexed as in docs:
 # https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-c/top/sparse-solver-routines/onemkl-pardiso-parallel-direct-sparse-solver-iface/pardiso-iparm-parameter.html
-const iparms = Dict([
+const iparms_dict = Dict([
     (0, 1), # in: do not use default
     (1, 2), # in: nested dissection from METIS
     (2, 0), # in: reserved, set to zero
@@ -436,4 +433,16 @@ const iparms = Dict([
     (61, 0), # reserved, set to zero 
     (62, 0), # out: size of the minimum OOC memory for factorization and sol 
     (63, 0), # reserved, set to zero 
-])
+])
+# MKL Pardiso solver IPARM control parameters formatted for LinearSolve.jl
+const iparms = collect(
+    (key + 1, iparms_dict[key]) for key in sort!(collect(keys(iparms_dict))))
+# LinearSolve.jl solver keyword arguments
+const cache_kwargs = (;
+    nprocs = 4,
+    verbose = true,
+    abstol = 1e-8,
+    reltol = 1e-8,
+    maxiter = 30,
+    iparm = iparms,
+    )