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
Issue with PETScSENS solver in dolfinx for two field problem dolfinx	2	510	August 11, 2022
Issue with PETScSENS solver in dolfinx dolfinx	2	318	August 10, 2022
Error in assigning function values in Matrix-free solution in dolfinx dolfinx	2	504	March 10, 2023
SNES solver fails when using the line search in FEniCSx dolfinx	9	2034	January 25, 2022
Error with complex numbers in dolfinx dolfinx	15	1227	August 19, 2022

Issue with dolfinx.fem.create_matrix

Related Topics