Updating dof of dolfinx.fem.function.Function using Midpoint time stepping

Hello all,

I am constructing a predictor-corrector method, using a calculated solution (C_n) from a linear solve, to then correct using a midpoint time stepping scheme. I am having trouble accessing the dof I need as the resulting data structures are a combination of ufl.algebra datatypes rather than fem.Function types.
I suspect that the biggest issues here are evaluating the nonlinear term g(C_mid) , flux term div(phi(C_mid)), and the source function f such that I can access their values at the dof. I believe that the evaluate and expression functions may be necessary here, but am confused on where to go from here.

Thank you!

Below is a MWE that produces the error

import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")

import ufl
from ufl import *
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import mesh, fem, io, nls, log , plot
from dolfinx.fem import (Expression, Function, FunctionSpace)
from numpy import linalg as LA

## Set up Function Space and initial parameters
u = as_vector([1.0,1.0])
t=0
dt = 1/16
domain = mesh.create_unit_square(MPI.COMM_WORLD, 4, 4 , mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("CG", 1))
fdim = domain.topology.dim - 1
tdim = domain.topology.dim
domain.topology.create_connectivity(fdim, tdim)

## Construct nodal coordinates
x = ufl.SpatialCoordinate(domain)


class exact_solution(): 
    def __init__(self, t):
        self.t = t
    def __call__(self, x):
         return np.exp(-self.t) * (x[0] - x[0]**2) * (x[1] - x[1]**2)

C_exact = exact_solution(t)
## Form source from manufactured solution
f = C_exact(x) * -1 * np.exp(-t) * (x[0] - x[0]**2) * (x[1] - x[1]**2)
print(type(f)) #ufl.algebra.product

def g(c): # Nonlinear Term (Scalar-Valued)
    return 1/((1+c)**2)
    
def g_array(c): # Possible fix for assignment to Ch.x.array?
    return 1/((1+c.x.array)**2)

def phi(c): #Flux Term (Vector-Valued)
        return c * u  - grad(c)

C_mid = fem.Function(V)
C = fem.Function(V)
C_n = fem.Function(V)
Ch = fem.Function(V)

C_n.interpolate(lambda x: x[0]) # known initial value made random for MWE
C.interpolate(lambda x: x[0]) #known value from predictor solve. This is a linear solve so it could be a trial function

C_mid.x.array[:] = (C_n.x.array + C.x.array)/2 # C_mid is still a fem.Function
## Evaluate nonlinear term and flux term at midpoint
g_mid = g(C_mid)
print(type(g_mid)) # ufl.algebra.division

phi_mid = phi(C_mid)
print(type(phi_mid)) # ufl.algebra.sum
## Midpoint TIme step
Ch = C_n + dt/g_mid*(f- div(phi_mid))
#Ch.x.array[:] = C_n.x.array + dt / g_mid.x.array * (f - div(phi_mid.x.array)) # Original Idea
print(Ch.x.array)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In [11], line 58
     56 Ch = C_n + dt/g_mid*(f- div(phi_mid))
     57 # Ch.x.array[:] = C_n.x.array + dt / g_mid.x.array * (f - div(phi_mid.x.array)) # Original Idea
---> 58 print(Ch.x.array)

AttributeError: 'Sum' object has no attribute 'x'

Referencing https://fenicsproject.discourse.group/t/dolfinx-fem-function-function-add-with-dolfinx-fem-function-function/9328

I replaced

Ch = C_n + dt/g_mid*(f- div(phi_mid))

with

expr = dolfinx.fem.Expression(C_n + dt/g_mid*(f- div(phi_mid)),V.element.interpolation_points())
Ch.interpolate(expr)

Is this on the right track?

Please use 3x` encapsulation of all code and error messages, ie.

```python
import dolfinx

# add code here
```
1 Like