Issue with dolfinx.fem.create_matrix

arianaqyp · January 26, 2022, 6:04am

I am trying to use dolfinx.fem.create_matrix, but its giving me an error. This code is based on the fenicsx tutorial for phase field analysis.

E_alpha = ufl.derivative(total_energy,alpha,ufl.TestFunction(V_alpha))
E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))

damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)

b = dolfinx.cpp.la.create_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
J = dolfinx.fem.create_matrix(damage_problem.a)

# Create Newton solver and solve
solver_alpha_snes = PETSc.SNES().create()
solver_alpha_snes.setType("vinewtonrsls")
solver_alpha_snes.setFunction(damage_problem.F, b)
solver_alpha_snes.setJacobian(damage_problem.J, J)
solver_alpha_snes.setTolerances(rtol=1.0e-9, max_it=50)
solver_alpha_snes.getKSP().setType("preonly")
solver_alpha_snes.getKSP().setTolerances(rtol=1.0e-9)
solver_alpha_snes.getKSP().getPC().setType("lu")

# We set the bound (Note: they are passed as reference and not as values)
solver_alpha_snes.setVariableBounds(alpha_lb.vector,alpha_ub.vector)

This is the error I am getting

TypeError                                 Traceback (most recent call last)
Input In [11], in <module>
      1 E_alpha = ufl.derivative(total_energy,alpha,ufl.TestFunction(V_alpha))
      2 E_alpha_alpha = ufl.derivative(E_alpha,alpha,ufl.TrialFunction(V_alpha))
----> 4 damage_problem = SNESProblem(E_alpha, E_alpha_alpha, alpha, bcs_alpha)
      6 b = dolfinx.cpp.la.create_vector(V_alpha.dofmap.index_map, V_alpha.dofmap.index_map_bs)
      7 J = dolfinx.fem.create_matrix(damage_problem.a)

File /workspaces/newfrac-core-numerics/snes_problem.py:30, in SNESProblem.__init__(self, F, J, u, bcs)
     26 self.u = u
     28 # Create matrix and vector to be used for assembly
     29 # of the non-linear problem
---> 30 self.A = dolfinx.fem.create_matrix(self.a)
     31 self.b = dolfinx.fem.create_vector(self.L)

File /usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/assemble.py:87, in create_matrix(a, mat_type)
     85 def create_matrix(a: FormMetaClass, mat_type=None) -> PETSc.Mat:
     86     if mat_type is None:
---> 87         return _cpp.fem.petsc.create_matrix(a)
     88     else:
     89         return _cpp.fem.petsc.create_matrix(a, mat_type)

TypeError: create_matrix(): incompatible function arguments. The following argument types are supported:
    1. (a: dolfinx::fem::Form<double>, type: str = '') -> mat

Invoked with: Form([Integral(CoefficientDerivative(CoefficientDerivative(Product(FloatValue(0.5), Inner(ComponentTensor(Product(Indexed(Sum(ComponentTensor(Product(Indexed(Sym(Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)), 3))), MultiIndex((Index(8), Index(9)))), Product(FloatValue(2.0), Division(Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 0), Product(FloatValue(2.0), Sum(FloatValue(1.0), Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 1)))))), MultiIndex((Index(8), Index(9)))), ComponentTensor(Product(Indexed(Identity(2), MultiIndex((Index(10), Index(11)))), Product(Division(Product(Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 0), Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 1)), Sum(FloatValue(1.0), Product(IntValue(-1), Power(Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 1), IntValue(2))))), Trace(Sym(Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)), 3)))))), MultiIndex((Index(10), Index(11))))), MultiIndex((Index(12), Index(13)))), Sum(FloatValue(1e-06), Power(Sum(IntValue(1), Product(IntValue(-1), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4))), IntValue(2)))), MultiIndex((Index(12), Index(13)))), Sym(Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)), 3))))), ExprList(*(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4),)), ExprList(*(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None),)), ExprMapping(*())), ExprList(*(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4),)), ExprList(*(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None),)), ExprMapping(*())), 'cell', Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), 'everywhere', {}, None), Integral(CoefficientDerivative(CoefficientDerivative(Product(Sum(Product(Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 3), Dot(Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4)), Grad(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4)))), Division(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4), Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 3))), Division(Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (), 2), FloatValue(2.6666666666666665))), ExprList(*(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4),)), ExprList(*(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None),)), ExprMapping(*())), ExprList(*(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4),)), ExprList(*(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None),)), ExprMapping(*())), 'cell', Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), 'everywhere', {}, None), Integral(CoefficientDerivative(CoefficientDerivative(Product(IntValue(-1), Dot(Constant(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), (2,), 4), Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2)), 3))), ExprList(*(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4),)), ExprList(*(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None),)), ExprMapping(*())), ExprList(*(Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 4),)), ExprList(*(Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None),)), ExprMapping(*())), 'cell', Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), 'everywhere', {}, None)])

dokken · January 26, 2022, 7:34am

As shown in: dolfinx/test_nonlinear_assembler.py at 611361eb24a6b8f6a9480b2034f53dc67feecf7e · FEniCS/dolfinx · GitHub
you need to call fem.form(E_alpha) and fem.form(E_alpha_alpha).

The reason for this is explained in:

github.com/FEniCS/dolfinx

Dynamically create Form classes to inherit from appropriate typed cpp base

FEniCS:main ← FEniCS:garth/dynamic-form

opened 04:50PM - 27 Dec 21 UTC

garth-wells

+910 -889

This change creates `fem.Form` instances with the correct `fem::Form<T>` base. T…his removes the need to access '`_cpp_object`'. Forms can be passed directly into pybind11 wrapped functions that expect a `fem::Form<T>`. This PR also requires that forms passed into lower-level functions, e.g. assembly, be explicitly instantiated. This fixes #1209 and is consistent with the DOLFINx design principle of favouring explicit over implicit (especially in lower level interfaces). Being explicit highlighted some cases in the test where the logic was not transparent. Previously one could do: ```cpp ... from dolfinx.fem import assemble_matrix a = inner(u, v) * dx A = assemble_matrix(a) ``` this change requires: ```cpp ... from dolfinx.fem import assemble_matrix, form a = form(inner(u, v) * dx) A = assemble_matrix(a) ``` Subsequent calls to `assemble_matrix(a)` will no longer involve the overhead of UFL and FFCx processing of the form, Even if the form is in a cache, the cost of pre-processing can in cases be noticeable. Moreover, this change is consistent with forthcoming changes for `Expression` that will deliberately not hide the behind-the-scenes UFL/FFCx processing. Higher level interfaces, e.g. `fem.LinearProblem`, are not affected.

and suggested input parameters to fem.form is highlighted in
https://jorgensd.github.io/dolfinx-tutorial/chapter4/compiler_parameters.html

Topic		Replies	Views
Error and I can't follow up with new fixes Errors	5	31	October 22, 2024
Error occuring in fem.form dolfinx	9	960	October 2, 2023
Issue with PETScSENS solver in dolfinx for two field problem dolfinx	2	584	August 11, 2022
How do I solve "Expecting scalar." error when defining NonlinearProblem? Errors dolfinx	1	183	March 12, 2024
Issue with PETScSENS solver in dolfinx dolfinx	2	351	August 10, 2022

Issue with dolfinx.fem.create_matrix

Related Topics