Heterogeneous constants in mixed spaces

General
1

Hi all! I am running a code to integrate a set of PDEs consisten of scalar fields on a spherical surface. One of the parameters, namely k10noise in my code, is heterogeneous in space, and a valueError arises when computing its product with my variable u3. An MME is provided below:

import pygmsh
import meshio
import numpy as np

from tqdm import tqdm
from scipy.spatial.distance import cdist
from scipy.linalg import cholesky

from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx import fem, io, default_real_type
from ufl import TestFunction, inner, dx, grad, split
from basix.ufl import element, mixed_element
from mpi4py import MPI
from petsc4py import PETSc

def generate_noise(coordinates, correlation_length):
    """
    Generate a Gaussian correlated noise field on the surface of a sphere.
    :param coordinates: Numpy array of spherical coordinates on the mesh.
    :param correlation_length: Correlation length for the Gaussian noise.
    :return: Interpolated noise field on the sphere.
    """
    # Calculate pairwise distances between points
    distances = cdist(coordinates, coordinates)
    
    # Exponential correlation function
    noise_covariance = np.exp(-distances / correlation_length)
    
    # Cholesky decomposition of the covariance matrix
    try:
        L = cholesky(noise_covariance, lower=True)
    except np.linalg.LinAlgError:
        # In case of numerical instability, fallback to regularization
        L = cholesky(noise_covariance + 1e-10 * np.eye(len(coordinates)), lower=True)
    
    # Generate uncorrelated Gaussian noise
    rng = np.random.default_rng(42)
    uncorrelated_noise = rng.standard_normal(coordinates.shape[0])
    
    # Apply Cholesky decomposition to generate correlated noise
    noise_values = L @ uncorrelated_noise
    
    # Return the noise values
    return noise_values

# Generate the spherical surface mesh with pygmsh
with pygmsh.geo.Geometry() as geom:
    radius = 1.0  # Radius of the sphere
    sphere = geom.add_ball([0, 0, 0], radius, mesh_size=0.035)
    mesh = geom.generate_mesh(dim=2)

# Filter out mixed cells and retain only triangles
triangle_cells = mesh.cells_dict["triangle"]

mesh = meshio.Mesh(points=mesh.points, cells=[("triangle", triangle_cells)],)

meshio.write("SphericalSurfaceMesh.xdmf", mesh)

# Load the filtered mesh into FEniCSx
with io.XDMFFile(MPI.COMM_WORLD, "SphericalSurfaceMesh.xdmf", "r") as infile:
    fenics_mesh = infile.read_mesh(name="Grid")

coordinates = fenics_mesh.geometry.x  # Nodal coordinates
x, y, z = coordinates[:, 0], coordinates[:, 1], coordinates[:, 2]

# Parameters
cell_radius = 75
Du = fem.Constant(fenics_mesh, 0.08/(cell_radius**2))
Dv = fem.Constant(fenics_mesh, 0.4/(cell_radius**2))
Dw = fem.Constant(fenics_mesh, 0.001/(cell_radius**2))

k0 = fem.Constant(fenics_mesh, 0.00625)
k1 = fem.Constant(fenics_mesh, 0.3125)
k2 = fem.Constant(fenics_mesh, 1.)
k3 = fem.Constant(fenics_mesh, 0.0625)
k4 = fem.Constant(fenics_mesh, 0.05625)
k5 = fem.Constant(fenics_mesh, 0.0625)
k6 = fem.Constant(fenics_mesh, 0.02083)
k7 = fem.Constant(fenics_mesh, 0.001875)
k8 = fem.Constant(fenics_mesh, 0.140625)
k9 = fem.Constant(fenics_mesh, 0.25)

alpha = fem.Constant(fenics_mesh, 1.) 
beta = fem.Constant(fenics_mesh, 0.3) 

unit = fem.Constant(fenics_mesh, 1.) 

noise_values = generate_noise(cell_radius*coordinates, 4)
k10noise = fem.Constant(fenics_mesh, 0.025)           #No noise
k10noise = fem.Constant(fenics_mesh, 0.025*noise_values) #Noise

# Define the function space for the scalar fields in the mixed space
P1 = element("Lagrange", fenics_mesh.basix_cell(), 1, dtype=default_real_type)
ME = fem.functionspace(fenics_mesh, mixed_element([P1, P1, P1]))

u = fem.Function(ME)            # Current solution
u_minus = fem.Function(ME)      # Previous solution
v1, v2, v3 = TestFunction(ME)   # Test functions

# Split mixed functions and assign initial conditions
u1, u2, u3 = split(u)
u1_minus, u2_minus, u3_minus = split(u_minus)

rng = np.random.default_rng(42)
u.sub(0).interpolate(lambda x: 1 + 1.*(rng.random(x.shape[1]))-0.5)
u.sub(1).interpolate(lambda x: 1 + 1.*(rng.random(x.shape[1]))-0.5)
u.sub(2).interpolate(lambda x: 1 + 1.*(rng.random(x.shape[1]))-0.5)
u.x.scatter_forward()

u_minus.x.array[:] = u.x.array[:]

# Time-stepping parameters and data storage
iterations = 20000
dt = 0.025
samples = 100

u_values = np.empty([coordinates[:, 0].shape[0], samples])

# Define nonlinear functions
def f(u1, u2, u3, k0, k1, k2, k3, k4, alpha, beta):
    return (k0+alpha*k1*np.pow(u1, 3)/(unit+k2*np.pow(u1, 2)))*u2-(k3+k4*(unit+beta)*u3)*u1

def g(u1, u2, u3, k0, k1, k2, k3, k4, k5, k6, alpha, beta):
    return k5-k6*u2-((k0+alpha*k1*np.pow(u1, 3)/(unit+k2*np.pow(u1, 2)))*u2-(k3+k4*(unit+beta)*u3)*u1)

def h(u1, u2, u3, k7, k8, k9, k10noise):
    return k7+k8*np.pow(u1, 2)/(unit+k9*np.pow(u1, 2))-k10noise*u3

# Construct the weak formulation
F1 = (inner(u1, v1) - inner(u1_minus, v1) + dt*Du*inner(grad(u1), grad(v1)) - dt*inner(f(u1, u2, u3, k0, k1, k2, k3, k4, alpha, beta), v1))*dx
F2 = (inner(u2, v2) - inner(u2_minus, v2) + dt*Dv*inner(grad(u2), grad(v2)) - dt*inner(g(u1, u2, u3, k0, k1, k2, k3, k4, k5, k6, alpha, beta), v2))*dx
F3 = (inner(u3, v3) - inner(u3_minus, v3) + dt*Dw*inner(grad(u3), grad(v3)) - dt*inner(h(u1, u2, u3, k7, k8, k9, k10noise), v3))*dx

F_total = F1 + F2 + F3

# Define the nonlinear problem
problem = NonlinearProblem(F_total, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

# Configure the solver
solver.convergence_criterion = "incremental"
solver.rtol = 1e-4
solver.report = True

# PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()
opts["ksp_type"] = "gmres"
opts["pc_type"] = "hypre"
opts["pc_hypre_type"] = "boomeramg"
opts["ksp_rtol"] = 1e-6

ksp.setFromOptions()

# Time-stepping loop
for it in tqdm(range(iterations)):
    # Solve the nonlinear problem
    n, converged = solver.solve(u)  # Use u1 as the overall solution
    if not converged:
        raise RuntimeError(f"Newton solver did not converge at iteration {it}")

    # Update previous solutions
    u_minus.x.array[:] = u.x.array[:] 
    
    if it%(iterations//samples) == 0:
        u_values[:, int(it/(iterations//samples))] = u.sub(0).collapse().x.array[:]

np.save("u_values.npy", u_values)

Giving the following error:

Traceback (most recent call last):
  File "/home/dlacasa/FEMSphere_CD.py", line 133, in <module>
    F3 = (inner(u3, v3) - inner(u3_minus, v3) + dt*Dw*inner(grad(u3), grad(v3)) - dt*inner(h(u1, u2, u3, k7, k8, k9, k10noise), v3))*dx
                                                                                           ~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/dlacasa/FEMSphere_CD.py", line 128, in h
    return k7+k8*np.pow(u1, 2)/(unit+k9*np.pow(u1, 2))-k10noise*u3
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~
  File "/home/dlacasa/miniconda3/envs/fenicsx-env/lib/python3.13/site-packages/ufl/exproperators.py", line 231, in _sub
    return Sum(self, -o)
  File "/home/dlacasa/miniconda3/envs/fenicsx-env/lib/python3.13/site-packages/ufl/algebra.py", line 45, in __new__
    raise ValueError("Can't add expressions with different shapes.")
ValueError: Can't add expressions with different shapes.

Thanks for your help!

2

Hi danixala5,

As the error indicates, there is a shape mismatch in the expression referenced, probably due to k10noise being an Ndarray of shape equal to your coordinates (from what I can understand of the generate_noise() function), whereas the other constants (k7, k8 etc.) you are summing in the expression are scalar valued expressions, so the expressions can’t be added together.

I think you could add it all together if you access the actual values instead, e.g.:
k7.value + k10noise.value
However, this would make the sum an Ndarray instead of a form expression.

Since k10noise is a spatially varying function and not a scalar constant, I would rather consider using a dolfinx.fem.Function object for it, instead of dolfinx.fem.Constant.

Cheers

3

It does compile and work now thank you! I thought that anything that does not change in time or in space would be better described by a constant. But your approach does indeed makes sense.

Thanks again!

1 Like