Minor

src/Solvers/petsc.jl

@@ -6,6 +6,10 @@ get_block_solver(::Val{:petsc_gmres_schwarz},params) = PETScLinearSolver(petsc_g
 get_block_solver(::Val{:petsc_gmres_amg},params)     = PETScLinearSolver(petsc_gmres_amg_setup)
 get_block_solver(::Val{:petsc_from_options},params)  = PETScLinearSolver()
 
+# MUMPS
+
+MUMPS_MAX_MEM_INCREASE = 40 # TODO: Create preferences?
+
 function petsc_mumps_setup(ksp)
   pc       = Ref{GridapPETSc.PETSC.PC}()
   mumpsmat = Ref{GridapPETSc.PETSC.Mat}()
@@ -17,9 +21,9 @@ function petsc_mumps_setup(ksp)
   @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat)
   # @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  4, 1)
   @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[],  7, 0)
+  @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 14, MUMPS_MAX_MEM_INCREASE)
   @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2)
   @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 1)
-  # @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6)
   @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL)
 end
 

test/dev/bugfix.jl

@@ -11,7 +11,7 @@ using GridapGmsh
 using GridapMHD
 
 np = (2,2)
-ranks = with_mpi() do distribute
+ranks = with_debug() do distribute
   distribute(LinearIndices((prod(np),)))
 end
 
@@ -20,12 +20,18 @@ mesh = Dict{Symbol,Any}(
   :base_mesh => "solid",
 )
 
+mesh = Dict{Symbol,Any}(
+  :mesher => :p4est_SG,
+  :base_mesh => "solid",
+  :num_refs => 0,
+)
+
 params = Dict{Symbol,Any}(
   :debug=>false,
   :solve=>true,
   :res_assemble=>false,
   :jac_assemble=>false,
-  :solver=> GridapMHD.default_solver_params(Val(:julia))
+  :solver=> GridapMHD.default_solver_params(Val(:badia2024))
 )
 params[:fluid] = Dict{Symbol,Any}(
   :domain => "fluid",
@@ -34,20 +40,25 @@ params[:fluid] = Dict{Symbol,Any}(
   :γ => 0.0,
   :f => VectorValue(0.0,0.0,0.0),
   :B => VectorValue(0.0,1.0,0.0),
-  :ζ => 0.0,
+  :ζ => 1.0,
   :convection => true,
 )
+u_in = GridapMHD.u_inlet(:parabolic,1.0,4.0,1.0)
 params[:bcs] = Dict{Symbol,Any}(
   :u => Dict(
     :tags => ["inlet", "wall"],
-    :values => [VectorValue(0.0, 0.0, 0.0), VectorValue(0.0, 0.0, 0.0)]
+    :values => [u_in, VectorValue(0.0, 0.0, 0.0)]
   ),
   :j => Dict(
     :tags => ["inlet", "outlet"],
     :values => [VectorValue(0.0,0.0,0.0), VectorValue(0.0,0.0,0.0)],
   ),
 )
 params[:solid] = Dict(:domain => "wall",:σ => 1.0)
+params[:fespaces] = Dict{Symbol,Any}(
+  :order_u => 2,
+  :order_j => 3,
+)
 
 model = GridapMHD.expansion_mesh(mesh,ranks,params)
 
@@ -63,19 +74,39 @@ mfs = GridapMHD._multi_field_style(params)
 V = MultiFieldFESpace([V_u,V_p,V_j,V_φ];style=mfs)
 U = MultiFieldFESpace([U_u,U_p,U_j,U_φ];style=mfs)
 
-res, jac = GridapMHD.weak_form(params)
+op = GridapMHD._fe_operator(U,V,params)
+xh = GridapMHD._allocate_solution(op,params)
+r = residual(op,xh)
+j = jacobian(op,xh)
 
-xh = zero(U);
 u = get_trial_fe_basis(U);
 v = get_fe_basis(V);
 
 contr = jac(xh,u,v)
 map(local_views(contr)) do a 
   for strian in Gridap.CellData.get_domains(a)
-    scell_mat = Gridap.CellData.get_contribution(a,strian)
-    for mat in scell_mat
-      println(mat)
-    end
+    cell_mat = Gridap.CellData.get_contribution(a,strian)
+    #for mat in cell_mat
+    #  println(mat)
+    #end
+    println(first(cell_mat))
   end
 end
+
+res, jac = GridapMHD.weak_form(params)
+
+# res, jac = GridapMHD.weak_form(params)
+# xh = zero(U);
+# u = get_trial_fe_basis(U);
+# v = get_fe_basis(V);
+# 
+# contr = jac(xh,u,v)
+# map(local_views(contr)) do a 
+#   for strian in Gridap.CellData.get_domains(a)
+#     scell_mat = Gridap.CellData.get_contribution(a,strian)
+#     for mat in scell_mat
+#       println(mat)
+#     end
+#   end
+# end
 #cell_mat = collect_cell_matrix(U, V, contr)