Hi Roberto,
I’m not amazingly surprised that a mixed function space with 23 components takes a long time. The matrix system is a block of 23^2 matrices and there will be some assembly infrastructure for each of them, even the ones which are zero. This is really pushing the boundaries of what we can do. It seems odd to me that there are so many real function space objects. Can you provide a brief maths write-up of what you are trying to do?
Regards,
David
From: <firedrake-bounces@imperial.ac.uk> on behalf of Roberto Federico Ausas <rfausas@icmc.usp.br>
Date: Wednesday, 29 April 2020 at 15:08
To: firedrake <firedrake@imperial.ac.uk>
Subject: [firedrake] very slow code when using a lot of functions belonging to Real spaces
Hello,
I have recently posted some questions related to the use of Real function spaces for some FSI problems I'm trying to solve.
The problem I'm facing now is that the solving phase is consuming absurd amount of time. It does not seem to be the PETSc resolution itself, but the assembly phase, since, once the KSP iterations begin they proceed reasonably fast.
The piece of code I believe is relevant is:
#------------------------------------------------------------------------
.
.
.
EU = VectorElement("Lagrange", 'triangle', orderU)
EP = FiniteElement("Lagrange", 'triangle', orderP)
EL = FiniteElement("Real", 'triangle', 0)
ES = []ndofs = 20
for i in range(ndofs):
ES.append(FiniteElement("Real", 'triangle', 0))
EM = MixedElement([EU,EP,EL,ES])
W_f = FunctionSpace(mesh, EM) # mesh has a few hundred nodes
sol = Function(W_f)
auxs = split(sol)
u = auxs[0]
p = auxs[1]
l = auxs[2]
s = []
for i in range(ndofs):
s.append(auxs[i+3])
auxT = TestFunctions(W_f)
v = auxT[0]
q = auxT[1]
r = auxT[2]
d = []
for i in range(ndofs):
d.append(auxT[i+3])
uT = u
for i in range(ndofs):
uT = uT + s[i]*phifluid[i]
# ( phifluid[:] are functions belonging to W_f that are zero everywhere except on some boundary nodes )
print(' |> - Fluid part')
Res = momentum_residual(uT, p, l, s, v, mesh) # Stokes operator
Res += mass_residual(uT, p, l, s, q, r, mesh) # Incompressibility restriction
print(' |> - Beam part')
aream1 = 1/assemble(Constant(1) * dx(mesh))
for i in range(ndofs):
Res += momentum_residual(uT, p, l, s, d[i]*phifluid[i], mesh)
force_b = Constant(force_beam(vphi, dphibeam[i]))
Res -= d[i]*force_b*aream1*dx
# Until this point the code runs reasonably fast
# Here starts the bottleneck
idsplit = '2'
for i in range(ndofs):
idsplit = idsplit+','+str(i+3) # idsplit ends up being 2,3,4,5,...
sp = {
"mat_type": "matfree",
"snes_type": "newtonls",
"snes_monitor": None,
"snes_converged_reason": None,
"snes_linesearch_type": "basic",
"snes_max_it": 8,
"snes_rtol": 1e-10,
"ksp_type": "fgmres",
"ksp_monitor_true_residual": None,
"ksp_max_it": 10,
"pc_type": "fieldsplit",
"pc_fieldsplit_type": "schur",
"pc_fieldsplit_schur_fact_type": "full",
"pc_fieldsplit_0_fields": "0,1",
"pc_fieldsplit_1_fields": idsplit,
"fieldsplit_0_ksp_type": "preonly",
"fieldsplit_0_pc_type": "python",
"fieldsplit_0_ksp_monitor_true_residual": None,
"fieldsplit_0_pc_python_type": "firedrake.AssembledPC",
"fieldsplit_0_assembled_pc_type": "lu",
"fieldsplit_0_assembled_pc_factor_mat_solver_type": "mumps",
"fieldsplit_0_assembled_mat_mumps_icntl_14": 200,
"mat_mumps_icntl_14": 200,
"fieldsplit_1_ksp_type": "gmres",
"fieldsplit_1_ksp_monitor_true_residual": None,
"fieldsplit_1_ksp_max_it": ndofs+1,
"fieldsplit_1_ksp_convergence_test": "skip",
"fieldsplit_1_pc_type": "none"
}solve(Res == 0, sol, bcs=bcs, J=Jac, solver_parameters=sp)
#------------------------------------------------------------------------
Any tip on what could be going on would be welcome.
Regards
--
Roberto F. Ausas
Instituto de Ciências Matemáticas e de Computação - ICMC/USP
São Carlos, SP - Brasil