A problem with NonlinearVariationalProblem

Ian · February 18, 2025, 6:13am

Hi，I’m a newhand in fenics. Now I want to do some dynamic simulation using fenics, but I met some problems in nonelinear problem solver. Here I’ve several PDE to describe a physical model. Some of them are Nonlinear problems and some of them are time-dependent. Anyway, I simplify the code and list it behind.

from fenics import *
mesh=BoxMesh(Point(0, 0, 0), Point(1e-6, 1e-6, 5*1e-7), 5,5,5)
T_element=FiniteElement("CG", tetrahedron, 1)
v_element=VectorElement("CG", tetrahedron, 1)
Mix=MixedElement([T_element,v_element])
W = FunctionSpace(mesh, Mix)
(T, v) = TrialFunctions(W)
(v_T, v_v) = TestFunctions(W)
w = Function(W)
def boundary(x, on_boundary):
    return on_boundary
bc_T = DirichletBC(W.sub(0), Constant(300), boundary)
bc_v = DirichletBC(W.sub(1), Constant((0,0,0)), boundary)
bcs = [bc_T, bc_v]
w_n = interpolate(Expression((
    "300",
    "0", "0", "0",
), degree=1),W)
T_n, v_n = w_n.split()
t=0
dt=1e-15
T_final = 1e-12
while t < T_final:
    t += float(dt)
    F_T=((T - T_n) / dt * v_T * dx
              + inner(grad(T), grad(v_T)) * dx +inner(v, v)*dx)
    F_v=(inner((v - v_n) / dt , v_v) * dx+ (1 / 2) * inner(grad(2 / 3  * T), v_v) * dx)
    F=F_T + F_v
    J = derivative(F, w)
    problem = NonlinearVariationalProblem(F, w, bcs=bcs, J=J)
    solver = NonlinearVariationalSolver(problem)
    newton_solver_params = solver.parameters["newton_solver"]
    newton_solver_params["linear_solver"] = "mumps"
    newton_solver_params["relative_tolerance"] = 1e-6 
    newton_solver_params["absolute_tolerance"] = 1e-8
    newton_solver_params["maximum_iterations"] = 50 
    solver.solve()
    w_n.assign(w)

The feedback from python is:

Traceback (most recent call last):
  File "/mnt/d/pypro/THz/newton problem rest.py", line 30, in <module>
    problem = NonlinearVariationalProblem(F, w, bcs=bcs, J=J)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py", line 157, in __init__
    F = Form(F, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
    ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 51, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 101, in ffc_jit
    return ffc.jit(ufl_form, parameters=p)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
    module = jit_build(ufl_object, module_name, parameters)
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
    module, signature = dijitso.jit(jitable=ufl_object,
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
    header, source, dependencies = generate(jitable, name, signature, params["generator"])
  File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
    code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
    return compile_ufl_objects(forms, "form", object_names,
  File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 185, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, kind, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in analyze_ufl_objects
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 89, in <genexpr>
    form_datas = tuple(_analyze_form(form, parameters)
  File "/usr/lib/python3/dist-packages/ffc/analysis.py", line 169, in _analyze_form
    form_data = compute_form_data(form,
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/compute_form_data.py", line 407, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/check_arities.py", line 177, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/check_arities.py", line 159, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
  File "/usr/lib/python3/dist-packages/ufl_legacy/corealg/map_dag.py", line 36, in map_expr_dag
    result, = map_expr_dags(function, [expression], compress=compress,
  File "/usr/lib/python3/dist-packages/ufl_legacy/corealg/map_dag.py", line 99, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/check_arities.py", line 48, in sum
    raise ArityMismatch("Adding expressions with non-matching form arguments {0} vs {1}.".format(_afmt(a), _afmt(b)))
ufl_legacy.algorithms.check_arities.ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ().

I checked my code and tried to read something， but I still don’t know where the problem is.

dokken · February 18, 2025, 6:41am

Ian:

while t < T_final:
    t += float(dt)
    F_T=((T - T_n) / dt * v_T * dx
              + inner(grad(T), grad(v_T)) * dx +inner(v, v)*dx)
    F_v=(inner((v - v_n) / dt , v_v) * dx+ (1 / 2) * inner(grad(2 / 3  * T), v_v) * dx)
    F=F_T + F_v
    J = derivative(F, w)

There are multiple issues here.

if the problem is non-linear, there shouldn’t be a TrialFunction (T or v) in the form, as you are forming the residual, not a matrix.
Secondly, this means that the formation of J doesn’t make sense, as w is not part of F
finally, there is no reason to define F, J, problem and solver within the for loop. This just adds extra work, with no benefit.

dokken · February 18, 2025, 6:47am

As a resource, I would consider:
http://jsdokken.com/FEniCS-workshop/src/form_compilation.html
which explains what can be pulled out of loops.

Ian · February 18, 2025, 10:02am

Thank you for your advice! I changed some of my code to fit the requirement.
However, there are still some problems using dolfinx.fem.petsc.LinearProblem. I used the code like the code in the link you offer Efficient usage of the Unified Form Language — FEniCS Workshop.
The feedback is that :

Traceback (most recent call last):
  File "/mnt/d/pypro/THz/newton problem rest.py", line 150, in <module>
    problem = dolfinx.fem.petsc.LinearProblem(a_compiled, L_compiled, u=w, bcs=[], petsc_options=petsc_options)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 788, in __init__
    self._A = create_matrix(self._a)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 176, in create_matrix
    return _cpp.fem.petsc.create_matrix(a._cpp_object)
AttributeError: 'Form' object has no attribute '_cpp_object'
Exception ignored in: <function LinearProblem.__del__ at 0x7f88a262e950>
Traceback (most recent call last):
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/petsc.py", line 831, in __del__
    self._solver.destroy()
AttributeError: 'LinearProblem' object has no attribute '_solver'

And I will simplify and list my code behind:

from fenics import *
import numpy as np
import dolfinx
import dolfinx.fem.petsc
import matplotlib.pyplot as plt
mesh=BoxMesh(Point(0, 0, 0), Point(1e-6, 1e-6, 5*1e-7), 5,5,5)
T_element=FiniteElement("CG", tetrahedron, 1)
v_element=VectorElement("CG", tetrahedron, 1)
Mix=MixedElement([T_element,v_element])
W = FunctionSpace(mesh, Mix)
(T, v) = TrialFunctions(W)
(v_T, v_v) = TestFunctions(W)
w = Function(W)
def boundary(x, on_boundary):
    return on_boundary
bc_T = DirichletBC(W.sub(0), Constant(300), boundary)
bc_v = DirichletBC(W.sub(1), Constant((0,0,0)), boundary)
bcs = [bc_T, bc_v]
w_n = interpolate(Expression((
    "300",
    "0", "0", "0",
), degree=1),W)
T_n, v_n = w_n.split()
t = dolfinx.fem.Constant(mesh, 0.0)
dt=1e-15
T_final = 1e-12
F_T=((T - T_n) / dt * v_T * dx
              + inner(grad(T), grad(v_T)) * dx +inner(v, v)*dx)
F_v=(inner((v - v_n) / dt , v_v) * dx+ (1 / 2) * inner(grad(2 / 3  * T), v_v) * dx)
F=F_T + F_v
a, L = system(F)
a_compiled = dolfinx.fem.form(a)
L_compiled = dolfinx.fem.form(L)
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
}
problem = dolfinx.fem.petsc.LinearProblem(a_compiled, L_compiled, u=w, bcs=[], petsc_options=petsc_options)
while t < T_final:
    t += dt
    problem.solve()
    w_n.assign(w)

dokken · February 18, 2025, 10:10am

This has now become a mixture of DOLFINx and DOLFIN (sorry, I missed that you was using legacy fenics when I read through it the first time).

I am not even sure how these can co-exist in the same environment.

If we take some steps back, and use legacy FEniCS, ignore the link that I previously gave you, and simply move the definitions out of the loop, and then change the form to use the Function rather than test functions, before differentiating to get the Jacobian

Ian · February 20, 2025, 7:48am

Thank you Prof. dokken, I solve the problem.

Topic		Replies	Views
What changed between using the linear solver and the nonlinear solver? General	2	443	November 7, 2022
Cannot solve nonlinear PDE General	3	219	March 19, 2024
Doubt with a non linear system problem General	13	294	December 10, 2023
Unable to define nonlinear problem Errors	3	435	March 27, 2023
Can anyone lend this noob a hand with his nonlinear Navier-Stokes solver?	1	620	May 15, 2019

A problem with NonlinearVariationalProblem

Related topics