Unclear errorcodes

[Thorinwasher]

What does this mean?

AssertionError                            Traceback (most recent call last)
Cell In[1], line 59
     56 a_cpp = dolfinx.fem.form(a)
     57 L_cpp = dolfinx.fem.form(L)
---> 59 A = multiphenicsx.fem.petsc.assemble_matrix_block(a_cpp, bcs=[],restriction=(restriction, restriction))
     60 F = multiphenicsx.fem.petsc.assemble_vector_block(L_cpp, a_cpp, bcs=[], restriction=restriction)
     62 questionMark = multiphenicsx.fem.petsc.create_vector_block(f_cpp, restriction=restriction)

File /usr/lib/python3.10/functools.py:889, in singledispatch.<locals>.wrapper(*args, **kw)
    885 if not args:
    886     raise TypeError(f'{funcname} requires at least '
    887                     '1 positional argument')
--> 889 return dispatch(args[0].__class__)(*args, **kw)

File /usr/local/lib/python3.10/dist-packages/multiphenicsx/fem/petsc.py:1417, in assemble_matrix_block(a, bcs, mat_type, diagonal, constants, coeffs, restriction)
   1382 @functools.singledispatch
   1383 def assemble_matrix_block(  # type: ignore[no-any-unimported]
   1384     a: typing.List[typing.List[dolfinx.fem.forms.form_types]],
   (...)
   1390         typing.Tuple[typing.List[mcpp.fem.DofMapRestriction], typing.List[mcpp.fem.DofMapRestriction]]] = None
   1391 ) -> petsc4py.PETSc.Mat:
   1392     """
   1393     Assemble bilinear forms into a new block PETSc matrix.
   1394 
   (...)
   1415         The assembled block PETSc matrix.
   1416     """
-> 1417     A = create_matrix_block(a, restriction, mat_type)
   1418     return assemble_matrix_block(A, a, bcs, diagonal, constants, coeffs, restriction)

File /usr/local/lib/python3.10/dist-packages/multiphenicsx/fem/petsc.py:317, in create_matrix_block(a, restriction, mat_type)
    294 def create_matrix_block(  # type: ignore[no-any-unimported]
    295     a: typing.List[typing.List[dolfinx.fem.forms.form_types]],
    296     restriction: typing.Optional[
    297         typing.Tuple[typing.List[mcpp.fem.DofMapRestriction], typing.List[mcpp.fem.DofMapRestriction]]],
    298     mat_type: typing.Optional[str] = None
    299 ) -> petsc4py.PETSc.Mat:
    300     """
    301     Create a block PETSc matrix which can be used to assemble the bilinear forms `a` with restriction `restriction`.
    302 
   (...)
    315         A PETSc matrix with a blocked layout that is compatible with `a` and restriction `restriction`.
    316     """
--> 317     return _create_matrix_block_or_nest(a, restriction, mat_type, mcpp.fem.petsc.create_matrix_block)

File /usr/local/lib/python3.10/dist-packages/multiphenicsx/fem/petsc.py:244, in _create_matrix_block_or_nest(a, restriction, mat_type, cpp_create_function)
    237 def _create_matrix_block_or_nest(  # type: ignore[no-any-unimported]
    238     a: typing.List[typing.List[dolfinx.fem.forms.form_types]],
    239     restriction: typing.Optional[
   (...)
    242     cpp_create_function: typing.Callable  # type: ignore[type-arg]
    243 ) -> petsc4py.PETSc.Mat:
--> 244     function_spaces = _get_block_function_spaces(a)
    245     rows, cols = len(function_spaces[0]), len(function_spaces[1])
    246     mesh = None

File /usr/local/lib/python3.10/dist-packages/multiphenicsx/fem/petsc.py:32, in _get_block_function_spaces(block_form)
     30 def _get_block_function_spaces(block_form: typing.List[typing.Any]) -> typing.List[typing.Any]:
     31     if isinstance(block_form[0], list):
---> 32         return _get_block_function_spaces_rank_2(block_form)
     33     else:
     34         return _get_block_function_spaces_rank_1(block_form)

File /usr/local/lib/python3.10/dist-packages/multiphenicsx/fem/petsc.py:78, in _get_block_function_spaces_rank_2(block_form)
     75 function_spaces = [function_spaces_0, function_spaces_1]
     76 assert all(a[i][j] is None or a[i][j].function_spaces[0] == function_spaces[0][i]
     77            for i in range(rows) for j in range(cols))
---> 78 assert all(a[i][j] is None or a[i][j].function_spaces[1] == function_spaces[1][j]
     79            for i in range(rows) for j in range(cols))
     80 return function_spaces

AssertionError:

I can send the code as well if you want, but my main issue with this is that there is no text saying the issue.

Maybe it should be something more like this?

assert all(a[i][j] is None or a[i][j].function_spaces[1] == function_spaces[1][j] for i in range(rows) for j in range(cols)), "functionspace missmatch or something"

I'm happy to contribute, but I'm a beginner at using dolfinx (any fem solver really), so I doubt I could help that much with fem related stuff.

There's many other spots where this is the scenario, just highlighting one I ran into.

[francesco-ballarin]

That probably means that you trial functions in the block form are not in the right column. Feel free to send an example to reproduce this (as minimal as you can).

I understand your point about not having error message in asserts, but I would give such a task such a low priority that I'd rather tell you now that I won't have the time to fix that.

[Thorinwasher]

I see where you're coming from, here's my try on a minimally reproducible code. Probably something stupid I missed to do. Yes, I haven't defined any bcs yet, but that should not be an issue when talking about functionspaces (although I bet the solution would not converge or similar things). Seems like I missed to do something with the definition of a, now that I'm looking at it again

What I'm trying to achieve: solve anything over multiple domains.

from mpi4py import MPI
from dolfinx import mesh
import dolfinx.fem
import multiphenicsx.fem
import ufl
import gmsh
from dolfinx.io import gmshio

kA = 1.0
kC = 2.0

#Creating the mesh with gmsh
if not gmsh.isInitialized():
    gmsh.initialize()
gmsh.model.add("battery2d")
model = gmsh.model()
model.occ.addRectangle(0, 0, 0, 1, 0.5, tag = 1)
model.occ.addRectangle(0, 0.5, 0, 1, 0.5, tag = 2)

model.occ.synchronize()

model.addPhysicalGroup(2, [1], 1)
model.addPhysicalGroup(2, [2], 2)
model.mesh.generate(dim=2)

(domain, cellTags, facetTags) = gmshio.model_to_mesh(model=model,comm=MPI.COMM_WORLD,rank=0,gdim=2)
gmsh.finalize()

anodeCells = cellTags.find(1)
cathodeCells = cellTags.find(2)

V = dolfinx.fem.FunctionSpace(domain, ("CG", 1))
VA = V.clone()
VC = V.clone()

dofs_anode = dolfinx.fem.locate_dofs_topological(VA, cellTags.dim, anodeCells)
dofs_cathode = dolfinx.fem.locate_dofs_topological(VC, cellTags.dim, cathodeCells)

restriction_anode = multiphenicsx.fem.DofMapRestriction(VA.dofmap, dofs_anode)
restriction_cathode = multiphenicsx.fem.DofMapRestriction(VC.dofmap, dofs_cathode)
restriction = [restriction_anode, restriction_cathode]

(uA, uC) = (ufl.TrialFunction(VA), ufl.TrialFunction(VC))
(vA, vC) = (ufl.TestFunction(VA), ufl.TestFunction(VC))

def poissonEquation_a(u : ufl.TrialFunction,v : ufl.TestFunction,k : float):
    return ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx

def poissonEquation_L(u : ufl.TrialFunction,v : ufl.TestFunction):
    return ufl.inner(300 * ufl.sin(20), v) * ufl.dx

a = [[poissonEquation_a(u=uA,v=vA,k=kA)],
     [poissonEquation_a(u=uC,v=vC,k=kC)]]

L = [poissonEquation_L(u=uA,v=vA),
     poissonEquation_L(u=uC,v=vC)]

a_cpp = dolfinx.fem.form(a)
L_cpp = dolfinx.fem.form(L)

A = multiphenicsx.fem.petsc.assemble_matrix_block(a_cpp, bcs=[],restriction=(restriction, restriction))
F = multiphenicsx.fem.petsc.assemble_vector_block(L_cpp, a_cpp, bcs=[], restriction=restriction)

[Thorinwasher]

Yep, I solved it: this was just a stupid definition of a, it should be:

a = [[poissonEquation_a(u=uA,v=vA,k=kA),None],
     [None,poissonEquation_a(u=uC,v=vC,k=kC)]]

I probably have some reading to do, as I don't understand why a would need to be defined in a matrix form here

[francesco-ballarin]

You are right that you do not need to use the block form here, since the equations for uA and for uC are fully decoupled, thus you could solve them separately. However, as soon as one of the off-diagonal terms is non-None, that will couple the A and the C problems.

[Thorinwasher]

Okay, I see. I still don't understand how uA and uC would couple if I add off diagonal terms. What equation is being solved when there's a matrix on the left side and a vector on the right side?

a * q = f

what is the vector q?

Sorry for asking too many questions : )

[francesco-ballarin]

For instance, something like

- \Delta u_A - k_{AC} \Delta u_C = f_A
- k_{CA} \Delta u_A - \Delta u_C = f_C

for k_{AC} and k_{CA} non zero.

[Thorinwasher]

In that scenario, would 'a' be defined as:

$$
\left(\begin{array}{cc} 
- \nabla u_A * \nabla v_A dx & -k_{AC} \nabla u_C * \nabla v_C dx\\
- k_{CA} \nabla u_A * \nabla v_A dx & - \nabla u_C * \nabla v_C dx
\end{array}\right)
$$ 

?

[francesco-ballarin]

There are two errors in your formulation:

$$
\left(\begin{array}{cc} 
- \nabla u_A * \nabla v_A dx & -k_{AC} \nabla u_C * \nabla v_{!!!! A !!!!} dx\\
- k_{CA} \nabla u_A * \nabla v_{!!!! C !!!!}  dx & - \nabla u_C * \nabla v_C dx
\end{array}\right)
$$ 

[Thorinwasher]

Aha. Thanks for the help <3