How to apply the radiation heat transfer?

variational formulation
1

Hello everyone. I am trying to add the radiation heat transfer to Robin bc

image
image213Ă—118 4.17 KB

When i add q_r part it crashes at A = assemble_matrix(bilinear_form, bcs)

        u_next = dolfinx.fem.Function(self.V)
        u_n = dolfinx.fem.Function(self.V)
        u_n.x.array[:] = 26
        Heat = dolfinx.fem.Function(self.V)
        u = ufl.TrialFunction(self.V)
        # u = dolfinx.fem.Function(self.V)
        v = ufl.TestFunction(self.V)
        a = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx 
        L = ( u_n + dt * Heat) * v * ufl.dx
        class BoundaryCondition():
            def __init__(self, type, marker, values):
                self._type = type
                if type == "Dirichlet":
                    facets = facet_tag.find(marker)
                    dofs = dolfinx.fem.locate_dofs_topological(V, fdim, facets)
                    self._bc = dolfinx.fem.dirichletbc(values, dofs, V)
                elif type == "Neumann":
                    self._bc = ufl.inner(values, v) * ds(marker)
                elif type == "Robin":
                    self._bc = [(values[1] * \
                        ufl.inner(u , v) + values[2]*values[3]*ufl.inner(u**4,v)) * ds(marker),(values[1] * \
                        ufl.inner(values[0], v) + values[2]*values[3]*ufl.inner(values[0]**4,v)) * ds(marker)]
                else:
                    raise TypeError(
                        "Unknown boundary condition: {0:s}".format(type))

            @property
            def bc(self):
                return self._bc

            @property
            def type(self):
                return self._type

        # Define the Dirichlet condition
        r = dolfinx.fem.Constant(self.domain, dolfinx.default_scalar_type(15 * 10**-2))
        s = dolfinx.fem.Constant(self.domain, dolfinx.default_scalar_type(20))
        H_c = dolfinx.fem.Constant(self.domain, dolfinx.default_scalar_type(0.8))
        SB_coeff = dolfinx.fem.Constant(self.domain, dolfinx.default_scalar_type(5.670374419*10**-1))
    
        bcs = []
        boundary_conditions = [ BoundaryCondition("Robin"    , 1, (s, r,SB_coeff,H_c)),
                               BoundaryCondition("Robin", 2, (s,r*3,SB_coeff,H_c)) ]
        for condition in boundary_conditions:
            if condition.type == "Dirichlet":
                bcs.append(condition.bc)
            else:
                a += condition.bc[0]
                L += condition.bc[1]
        bilinear_form = fem.form(a)
        linear_form = fem.form(L)
        A = assemble_matrix(bilinear_form, bcs)

Here is the error message

Traceback (most recent call last):
  File "/home/arshidin/Documents/Fenicsx/fenicsx/math_model/matmodel_example.py", line 115, in <module>
    Model.solve_trajectory_load_meshtags(np.array(points),False)
  File "/home/arshidin/Documents/Fenicsx/fenicsx/math_model/matmodel.py", line 550, in solve_trajectory_load_meshtags
    bilinear_form = fem.form(a)
                    ^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 337, in form
    return _create_form(form)
           ^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 331, in _create_form
    return _form(form)
           ^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 254, in _form
    ufcx_form, module, code = jit.ffcx_jit(
                              ^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/dolfinx/jit.py", line 62, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/dolfinx/jit.py", line 212, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 205, in compile_forms
    impl = _compile_objects(
           ^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 330, in _compile_objects
    _, code_body = ffcx.compiler.compile_ufl_objects(
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ffcx/compiler.py", line 108, in compile_ufl_objects
    analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ffcx/analysis.py", line 94, in analyze_ufl_objects
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ffcx/analysis.py", line 94, in <genexpr>
    form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ffcx/analysis.py", line 180, in _analyze_form
    form_data = ufl.algorithms.compute_form_data(
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py", line 427, in compute_form_data
    check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ufl/algorithms/check_arities.py", line 213, in check_form_arity
    check_integrand_arity(itg.integrand(), arguments, complex_mode)
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ufl/algorithms/check_arities.py", line 194, in check_integrand_arity
    arg_tuples = map_expr_dag(rules, expr, compress=False)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ufl/corealg/map_dag.py", line 35, in map_expr_dag
    (result,) = map_expr_dags(
                ^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ufl/corealg/map_dag.py", line 101, in map_expr_dags
    r = handlers[v._ufl_typecode_](v)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/arshidin/miniconda3/envs/dolfinx-env/lib/python3.11/site-packages/ufl/algorithms/check_arities.py", line 48, in nonlinear_operator
    raise ArityMismatch(
ufl.algorithms.check_arities.ArityMismatch: Applying nonlinear operator Power to expression depending on form argument v_1.

i have tried using function u = dolfinx.fem.Function(self.V)

but i get this error

    A = _cpp.fem.petsc.create_matrix(a._cpp_object)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Cannot create sparsity pattern. Form is not a bilinear.

This code and solution similar to diffusion_code.py time dependent problem
Have someone faced similar issues?

2

Since I can’t run your example it’s difficult to give advice. It’s likely your issue lies here. The error you’re getting is related to the fact that a does not correspond to a bilinear formulation. I.e. you’re mixing up linear and bilinear terms via your Robin condition.

There are tools in UFL which help extract the “left”- and “right”-hand side of a variational formulation. But for your own sense of rigour, it may be worth carefully separating these terms in your formulation.

3

With radiation boundary conditions, you end up with a nonlinear problem, you should therefore formulate your heat problem with a nonlinear residual form, not a linear and a bilinear form like you do here. See problems such as Implementation — FEniCSx tutorial
In this demo, the nonlinear residual is

F = q(uh) * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - f * v * ufl.dx

Adapting to your problem in the stationary case would give something such as:

F = q(uh) * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - const * (uh ** 4 - u0 ** 4) * v * ufl.ds
1 Like
4

Thans for pointing that out, i successfully runed my code, however i wasn`t satisfied with the speed so i switched to a linear approximation of radiation boundary condition
h_{rad}=4igma psilon T_{char}^3