Nonlinear function per physical region

General
1

Hello,

I’d like to define a nonlinear function which is function of the trial function, per physical region.

For instance, I have c_a(u) = 1+u^2 and c_b = 0.1+u^4 for respectively physical regions \Omega_a and \Omega_b s.t. \Omega = \Omega_a \cup \Omega_b.
I tried to use fem.Function interpolating fem.Expression’s over their respective domain, however it does not yield the expected behavior

DG0 = fem.functionspace(domain, ("DG", 0))
c = fem.Function(DG0) 
q = lambda u: 1+u*u
c.interpolate(
  fem.Expression( q(uh), DG0.element.interpolation_points ),
  mesh_data.cell_tags.find(one) #Omega_a, one is a 2d physical tag
)
c.interpolate(
  fem.Expression( q(uh), DG0.element.interpolation_points ),
  mesh_data.cell_tags.find(two) #Omega_b, two is a 2d physical tag
)
# F = c * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - f * v * ufl.dx

I give below an adapted example from Jørgen Dokken’s Nonlinear Poisson tutorial.

Running the below code with --single yields

image
image761×674 115 KB

while without

image
image761×674 100 KB

Since here the below code interpolates a single nonlinear function over the whole domain, both runs shall give the same result.
I guess the interpolation projects once the field, removing the nonlinearity of the weak formulation.

How can I define a nonlinear weak formulation whose nonlinear law (function) is defined by physical region?

import gmsh
import numpy as np
import ufl
from dolfinx import fem, io, default_scalar_type, mesh, plot
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from mpi4py import MPI
import pyvista

import argparse

def main(single: bool):

    gmsh.initialize()

    rectangle = gmsh.model.occ.addRectangle(-1, -1, 0, 2, 2)
    disk = gmsh.model.occ.addDisk(0, 0, 0, .5, .5)

    out, _ = gmsh.model.occ.fragment([(2, rectangle)], [(2, disk)])
    dim2 = [ o[1] for o in out if o[0]==2 ]

    gmsh.model.occ.synchronize()

    one = gmsh.model.add_physical_group(2, [dim2[0]], name="one")
    two = gmsh.model.add_physical_group(2, [dim2[1]], name="two")

    gmsh.model.mesh.generate(2)

    # What following is from
    # Author: Jørgen S. Dokken, "Nonlinear poisson tutorial"
    def q(u):
        return 1 + u**2


    mesh_data = io.gmsh.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0)
    domain = mesh_data.mesh
    x = ufl.SpatialCoordinate(domain)
    u_ufl = 1 + x[0] + 2 * x[1]
    f = -ufl.div(q(u_ufl) * ufl.grad(u_ufl))

    V = fem.functionspace(domain, ("Lagrange", 1))


    def u_exact(x):
        return eval(str(u_ufl))


    u_D = fem.Function(V)
    u_D.interpolate(u_exact)
    fdim = domain.topology.dim - 1
    boundary_facets = mesh.locate_entities_boundary(
        domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool)
    )
    bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))

    uh = fem.Function(V)
    v = ufl.TestFunction(V)
    DG0 = fem.functionspace(domain, ("DG", 0))
    c = fem.Function(DG0)
    if single:
        c = q(uh)
    else: # the comments shows what is intended, the run code shows what is wrong
        qbis = lambda u: u**4+0.1
        c.interpolate(
            fem.Expression( q(uh), DG0.element.interpolation_points ),
            None#mesh_data.cell_tags.find(one)
        )
        #c.interpolate(
        #    fem.Expression( qbis(uh), DG0.element.interpolation_points ),
        #    mesh_data.cell_tags.find(two)
        #)
    #_end if-else
    F = c * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - f * v * ufl.dx


    petsc_options = {
        "snes_type": "newtonls",
        "snes_linesearch_type": "none",
        "snes_atol": 1e-6,
        "snes_rtol": 1e-6,
        "snes_monitor": None,
        "ksp_error_if_not_converged": True,
        "ksp_type": "gmres",
        "ksp_rtol": 1e-8,
        "ksp_monitor": None,
        "pc_type": "hypre",
        "pc_hypre_type": "boomeramg",
        "pc_hypre_boomeramg_max_iter": 1,
        "pc_hypre_boomeramg_cycle_type": "v",
    }

    problem = NonlinearProblem(
        F,
        uh,
        bcs=[bc],
        petsc_options=petsc_options,
        petsc_options_prefix="nonlinpoisson",
    )

    problem.solve()
    converged = problem.solver.getConvergedReason()
    num_iter = problem.solver.getIterationNumber()
    assert converged > 0, f"Solver did not converge, got {converged}."
    print(
        f"Solver converged after {num_iter} iterations with converged reason {converged}."
    )


    u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)

    u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
    u_grid.point_data["u"] = uh.x.array.real
    u_grid.set_active_scalars("u")
    u_plotter = pyvista.Plotter()
    u_plotter.add_mesh(u_grid, show_edges=True)
    u_plotter.view_xy()
    u_plotter.show()

    gmsh.finalize()
    exit()
#_end def main()

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument("--single", help="If true, run code with a single nonlinear function, otherwise tries to run the code with two nonlinear function", action="store_true")
    args = parser.parse_args()
    main(args.single)
2

It is not entirely clear to me what problem you’re trying to solve. Could you confim you’re trying to solve the weak form:

\int_\Omega c(u) \nabla u \cdot \nabla v \, d\Omega = \int_\Omega f \, v \, d\Omega

with c(u) defined per your statement?

Then I’d just split up that c function into two, and split up the integration domain into two:

\int_{\Omega_a} c_a(u) \nabla u \cdot \nabla v \, d\Omega + \int_{\Omega_b} c_b(u) \nabla u \cdot \nabla v \, d\Omega = \int_\Omega f \, v \, d\Omega

For each of the two integrals, you can use exactly the approach from A nonlinear Poisson equation — FEniCSx tutorial and integrating over subdomains is as simple as *dx(one) or *dx(two), assuming you specified the subdomain information in the dx measure (see, e.g. Electromagnetic scattering from a sphere (axisymmetric) — DOLFINx 0.11.0.dev0 documentation )

3

Thanks for your reply and taking the time to understand my issue.

The weak form is only an example to embody what I need, i.e. a way to define many nonlinear laws that can be arbitrarily distinct on each physical subdomain.

I agree that splitting the integral is a way.
However, I expected something less intrusive, where the weak form is a single integral over the whole domain, and under the hood dolfinx respectively uses the nonlinear functions.

Do you confirm that splitting the integral domain is the only way?

4

I guess you could add ufl conditionals on your c definition, but honestly I’d say that splitting the integral is entirely natural. Anything else feels like a hack. You could even do it procedurally. Something like:

c_relations={domain_id_1:cfunc_1, domain_id_2:cfunc_2, ...}
for domain_id, c_func in c_relations.items():
   form += c(u)*ufl.dot(ufl.grad(u),ufl.grad(v))*dx(domain_id)

The problem with your own implementation is that once the c(u) is interpolated onto some other Function, the dependency on u is lost. You really ought to work within the ufl language to keep the knowledge of the u Function in the form.

5

Ok, what you propose corresponds to my fallback solution.

Thanks for your time.

6

You can use ufl.conditional to link different non-linear behaviors to different regions, based on physical quantities, or coordinates. of course, if the behavior is based on cell markers, the solution of @Stein is preferred.

1 Like
7

Late to the party but another option that you could consider, that doesn’t require splitting the variational form or filling it with conditionals, is to assign the dof values of your coefficient directly. Given a subdomain with tag tag, you can do

dofs = domain_tags.find(tdim, tag)
c.x.array[dofs] = 1 + u.x.array[dofs]**2

This works if the functions c and u belong to the same function space; otherwise, you first need to interpolate but that is automagic.