How to define solution dependable material properties using data points in DOLFINx?

evdkmv · September 15, 2021, 1:10pm

Dear Community,

I am trying to set nonlinear material properties, which depend on the solution in DOLFINx (version: 0.1.1.0). The task is pretty simple if polynomials are used. It is demonstrated in the following code (2D rectangle is cooling down during some time through heat exchange with the ambient):

import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [50, 50])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))

# Define initial conditions
u_n = dolfinx.Function(V)
with u_n.vector.localForm() as loc:
    loc.set(1000)

#Define solution variable
uh = dolfinx.Function(V)
with uh.vector.localForm() as loc:
    loc.set(1000)

v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)

## Material properties
# Definition using polynom
def heat_cond_polynom(u):
    return 0.00004*u + 0.02 

t = 0
dt=0.1
#Time loop
for i in range(1,6):
    t += dt
    print("time: " + str(np.round(t, 2)))

    q_ = heat_cond_polynom(uh)

    F = 0.54 * 0.0078 * uh * v * ufl.dx + q_*dt*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n + dt*f) * v * ufl.dx

    # Boundary conditions (Robin on all sides)
    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    marker = 1 
    facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 2) | np.isclose(x[1], 0) | np.isclose(x[1], 2))
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
    facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
    facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
    F += dt * 0.03 * ufl.inner(uh-20, v) * ds(marker)

    #Problem definition 
    problem = dolfinx.fem.NonlinearProblem(F, uh)
    solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
    n, converged = solver.solve(uh)
    assert(converged)
    print(f"Number of interations: {n:d}")

    #Print max temperature
    print(np.max(np.round(np.array(uh.vector.array),2)))

    # Update solution at previous time step (u_n)
    with uh.vector.localForm() as loc, u_n.vector.localForm() as loc_n:
        loc.copy(result=loc_n)

output:

time: 0.1
Number of interations: 3
temperature: 556.22
time: 0.2
Number of interations: 3
temperature: 298.18
time: 0.3
Number of interations: 3
temperature: 164.14
time: 0.4
Number of interations: 2
temperature: 94.73
time: 0.5
Number of interations: 2
temperature: 58.7

However, i need to set material properties using data points rather than polynomials. Just defining e.g. the interp1d function instead of the polynom results in “ValueError: object arrays are not supported”. There is a number of discussions regarding this issue in dolfin such as RemDelaporteMathurin and freya. Combining approaches suggested in these threads the following code is obtained:

import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [50, 50])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))

# Define initial conditions
u_n = dolfinx.Function(V)
with u_n.vector.localForm() as loc:
    loc.set(1000)

#Define solution variable
uh = dolfinx.Function(V)
with uh.vector.localForm() as loc:
    loc.set(1000)

v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)

## Material properties
# Definition using interpolation 
def heat_cond_interpolation(u):
    heat_conduction_data = sp.interp1d([0,1000],[0.02,0.06])
    return heat_conduction_data(u) 

t = 0
dt=0.1
#Time loop
for i in range(1,6):
    t += dt
    print("time: " + str(np.round(t, 2)))

    q_ = dolfinx.Function(V)
    q_.vector.array[:] = heat_cond_interpolation(np.array(uh.vector.array))

    F = 0.54 * 0.0078 * uh * v * ufl.dx + q_*dt*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n + dt*f) * v * ufl.dx

    # Boundary conditions (Robin on all sides)
    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    marker = 1 
    facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 2) | np.isclose(x[1], 0) | np.isclose(x[1], 2))
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
    facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
    facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
    F += dt * 0.03 * ufl.inner(uh-20, v) * ds(marker)

    #Problem definition 
    problem = dolfinx.fem.NonlinearProblem(F, uh)
    solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
    n, converged = solver.solve(uh)
    assert(converged)
    print(f"Number of interations: {n:d}")

    #Print max temperature
    print(np.max(np.round(np.array(uh.vector.array),2)))

    # Update solution at previous time step (u_n)
    with uh.vector.localForm() as loc, u_n.vector.localForm() as loc_n:
        loc.copy(result=loc_n)

which gives the following output and it is obvious that the temperatures as well as the number of iterations are not the same as they were while using the polynom:

time: 0.1
Number of interations: 1
temperature: 516.69
time: 0.2
Number of interations: 1
temperature: 269.04
time: 0.3
Number of interations: 1
temperature: 147.52
time: 0.4
Number of interations: 1
temperature: 85.95
time: 0.5
Number of interations: 1
temperature: 54.17

Basically, what i need is to reproduce the output of the first code using the material definition (Definition using interpolation) from the second code. There are other related topics such as RemDelaporteMathurin. However it does not work anymore since the “UserExpression” class has been deprecated in DOLFINx. I tried to come up with another solution using recommendations from Dokken, but have not succeeded. It has been more than two years since RemDelaporteMathurin posted his question. I have been hoping that during this time a simple solution to that issue appeared and is now available in DOLFINx.

Thank you

Sincerely
Anton

conpierce8 · September 20, 2021, 8:11pm

Hi Anton,

One possible solution would be to write your own Newton solver, in which you interpolate the material properties at each Newton iteration based on the current value of uh. See below, where I have adapted this solver to DOLFINx. On my machine, the same max temperature is achieved by both the “polynomial” and the “interpolation” method, although the number of Newton steps is different.

Note that this only works in serial, however.

import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
import scipy.interpolate as sp

mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([2, 2, 0])], [50, 50])
V = dolfinx.FunctionSpace(mesh, ("CG", 1))

# Define initial conditions
u_n_p = dolfinx.Function(V)
with u_n_p.vector.localForm() as loc:
    loc.set(1000)
u_n_i = dolfinx.Function(V)
with u_n_i.vector.localForm() as loc:
    loc.set(1000)

#Define solution variable
uh_p = dolfinx.Function(V)
uh_i = dolfinx.Function(V)
with uh_p.vector.localForm() as loc:
    loc.set(1000)
with uh_i.vector.localForm() as loc:
    loc.set(1000)

v = ufl.TestFunction(V)
f = dolfinx.Constant(mesh, 0)

## Material properties
# Definition using polynom
def heat_cond_polynom(u):
    return 0.00004*u + 0.02 

# Definition using interpolation 
def heat_cond_interpolation(u):
    heat_conduction_data = sp.interp1d([0,1000],[0.02,0.06])
    return heat_conduction_data(u) 

### NewtonSolver implementation adapted from javjiv4
#   https://fenicsproject.discourse.group/t/custom-newton-solver-problem-with-dirichlet-conditions/2329/5
def newtonSolve(G, dG, u, bcs, q, q_interp):
    """
    Solves the nonlinear system of equations G(u; v) == 0 using Newton's method:
        G'(u; du, v) = -G(u; v)
    where the nonlinearity occurs due to the coefficient `q` that depends on
    `u`.
    
    Parameters
    ==========
    `G` :: semilinear form
    `dG` :: Jacobian of G
    `u` :: `Function` to be solved for
    `bcs` :: `DirichletBC`s for the problem
    `q` :: `Function` or list of functions that occurs as a coefficient(s) in 
            `G` and `dG`, which depends on `u` but has to be evaluated
            numerically by interpolation
    `q_interp` :: Python function or list of functions of a vector that
             interpolates for q
    """
    
    if isinstance(q, list):
        if isinstance(q_interp, list):
            if not len(q) == len(q_interp):
                raise ValueError("q and q_interp must have same len()")
        else:
            raise ValueError("Mismatched types for q and q_interp")
    elif isinstance(q, dolfinx.Function):
        if not callable(q_interp):
            raise ValueError("If q is a Function, q_interp must be callable")
    else:
        raise ValueError("q must be a Function or list of Functions")
    
    tol = 1e-3
    E = 1.
    it = 0
    max_it = 10
    
    for bc in bcs:
        bc.apply(u.vector())
        bc.homogenize()
    
    while E > tol and it < max_it:
        # Interpolate material properties using current value of `u`
        if isinstance(q, list):
            for qi, q_int_i in zip(q, q_interp):
                qi.vector.array[:] = q_int_i(np.array(u.vector.array))
        else:
            q.vector.array[:] = q_interp(np.array(u.vector.array))
        
        x = dolfinx.Function(u.function_space)
        problem = dolfinx.fem.LinearProblem(dG, -G, bcs, x)
        problem.solve()
        
        error = (x)**2 * ufl.dx
        E = np.sqrt(dolfinx.fem.assemble_scalar(error))
        # print('Iteration:', str(it), '; Error: %3.3e' % E)
        if E < tol:
            break
        
        # Update
        u.vector.array[:] = u.vector.array[:] + x.vector.array[:]
        it += 1
    return it, it < max_it

t = 0
dt=0.1
#Time loop
for i in range(1,6):
    t += dt
    print("time: " + str(np.round(t, 2)))
    
    q_p = heat_cond_polynom(uh_p)
    q_i = dolfinx.Function(V)
    q_i.vector.array[:] = heat_cond_interpolation(np.array(uh_i.vector.array))

    Fp = 0.54 * 0.0078 * uh_p * v * ufl.dx + q_p*dt*ufl.dot(ufl.grad(uh_p), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n_p + dt*f) * v * ufl.dx
    Fi = 0.54 * 0.0078 * uh_i * v * ufl.dx + q_i*dt*ufl.dot(ufl.grad(uh_i), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n_i + dt*f) * v * ufl.dx

    # Boundary conditions (Robin on all sides)
    facet_indices, facet_markers = [], []
    fdim = mesh.topology.dim - 1
    marker = 1 
    facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 2) | np.isclose(x[1], 0) | np.isclose(x[1], 2))
    facet_indices.append(facets)
    facet_markers.append(np.full(len(facets), marker))
    facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
    facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
    sorted_facets = np.argsort(facet_indices)
    facet_tag = dolfinx.MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
    ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag)
    Fp += dt * 0.03 * ufl.inner(uh_p-20, v) * ds(marker)
    Fi += dt * 0.03 * ufl.inner(uh_i-20, v) * ds(marker)

    #Problem definition - polynomial
    problem = dolfinx.fem.NonlinearProblem(Fp, uh_p)
    solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
    n, converged = solver.solve(uh_p)
    assert(converged)
    print(f"Polynomial number of iterations: {n:d}")
    
    #Problem definition - interpolation
    u = ufl.TrialFunction(V)
    J = ufl.derivative(Fi, uh_i, u)
    
    n, converged = newtonSolve(Fi, J, uh_i, [], q_i, heat_cond_interpolation)
    assert(converged)
    print(f"Interpolation number of iterations: {n:d}")

    #Print max temperature
    print("Max temp, polynomial:",np.max(np.round(np.array(uh_p.vector.array),2)))
    print("Max temp, interpolation:",np.max(np.round(np.array(uh_i.vector.array),2)))

    # Update solution at previous time step (u_n)
    with uh_p.vector.localForm() as loc, u_n_p.vector.localForm() as loc_n:
        loc.copy(result=loc_n)
    with uh_i.vector.localForm() as loc, u_n_i.vector.localForm() as loc_n:
        loc.copy(result=loc_n)

EDIT: updated my code to deal with the case of multiple nonlinear coefficients q.

evdkmv · September 23, 2021, 10:12am

Hi Connor,

many thanks for your answer. The code works as expected. Following your idea I will definitely try to explore more about writing own solvers.

However, since the “polynom” definition of material properties works smoothly with existing solvers, I have a strong feeling that the “interpolation” definition should not require writing a modified solver (perhaps i am not right here). Both functions return the same (or almost the same, see ufl vs numpy). My main idea is to modify the “interpolation function” in a way that the output would be ok for the existing solver (now the solver gives an error: “AttributeError: ‘numpy.ndarray’ object has no attribute ‘subdomain_data’”).

What do you think, would it be possible?

Regards
Anton

conpierce8 · September 23, 2021, 1:45pm

Hi Anton,

I’m not a FEniCS developer and so I can’t make an authoritative statement on this, but I believe that by using interpolate, your function q_ has no symbolic link to the function uh. The NewtonSolver only updates the dof values of uh at each Newton iteration, and since there is no symbolic link to q_, the values of q_ remain constant. The Newton solver thus solves this problem as if it were a linear problem.

The reason the polynom function works is because it is equivalent to the following form, in which the nonlinearity is symbolically defined:

F = 0.54 * 0.0078*uh*v*ufl.dx + (0.00004*uh+0.02)*dt*ufl.dot(ufl.grad(uh), ufl.grad(v))*ufl.dx - (0.54*0.0078*u_n + dt*f) * v * ufl.dx

I’m not aware of a way to interpolate a symbolic argument, which is why I suggested defining your own solver.

evdkmv · October 5, 2021, 2:29pm

Okay, now I see where the problem is. Thank you for the explanation.
Looks like your suggestion is indeed a proper solution to this issue

Topic		Replies	Views
Vector/Tensor Material Properties dolfinx	3	1123	June 30, 2020
How to set material property as a function of split unknown variable? Legacy DOLFIN material-properties	7	180	December 15, 2023
Update Material Coordinates within Cells dolfinx	1	104	December 15, 2023
Dolfinx matrix-valued Expression dolfinx	2	930	March 11, 2022
Initial Conditions and Function update on time in mixed Element dolfinx	2	278	December 5, 2023

How to define solution dependable material properties using data points in DOLFINx?

Related Topics