A coupled system with multiple variables

variational formulation
1

Hello everyone, if it is a coupled system of convection diffusion equations with multiple variables, how should we define a variational equation? Is it simply adding several equations? How can we obtain its linear and bilinear formats?

I discovered the variational equation form of a coupled system in Fenics’s book tutorial and wrote a similar code, hoping to obtain its a and L by calling lhs() and rhs().

But this seems to be impractical.

F1 = ((u_1 - u_n1) / dt)*v_1*dx + dot(ve1, grad(u_1))*v_1*dx + D1*dot(grad(u_1), grad(v_1))*dx
    F2 = ((u_2 - u_n2) / dt)*v_2*dx + dot(ve2, grad(u_2))*v_2*dx + D2*dot(grad(u_2), grad(v_2))*dx
    F3 = ((u_3 - u_n3) / dt)*v_3*dx + dot(ve3, grad(u_3))*v_3*dx + D3*dot(grad(u_3), grad(v_3))*dx
    F = F1 + F2 + F3
    a = lhs(F1)
    L = rhs(F)

The running result is:
ValueError: Found Argument in , this is an invalid expression.

2

Without a minimal reproducible example it is hard to give any sort of guidance.
There are in general many ways to form such PDEs, depending on how you define your function-spaces.
For instances, in your snippet above, I don’t understand how aa=lhs(F1), while L=rhs(F).

Also note that you could simply make a python function:

def my_reaction(u, u_n, ve, v, D, dt):
    return ((u - u_n) / dt)*v*dx + dot(ve, grad(u))*v*dx + D*dot(grad(u), grad(v))*dx

F = my_reaction(u_1, u_n1, ve1, v_1, D1, dt) +  my_reaction(u_2, u_n2, ve2, v_2, D2, dt) +  my_reaction(u_3, u_n3, ve3, v_3, D3, dt)
3

I apologise for my unclear expression. Please allow me additional time to attempt to resolve the issue independently or modify the MWE.

4

Thank you for the reminder.

5

Dear dokken,
I have attempted to rebuild the code, and it seems to be running correctly now.
However, I am not certain if the usage of many APIs is correct, especially in the mixed function space.
Therefore, I have tried my best to extract MWE, as shown below.
I really hope to receive your confirmation so that I can make further settings and carry out the next simulation work.

from dolfinx import mesh, fem, plot
from dolfinx.fem.petsc import LinearProblem
import ufl
from ufl import (dot, dx, lhs, grad, rhs, sqrt, exp)
import numpy
from petsc4py import PETSc
from mpi4py import MPI
import pyvista

domain = mesh.create_rectangle(MPI.COMM_WORLD, [numpy.array([-1, 0]), numpy.array([1, 5])],
                               [20, 50], mesh.CellType.triangle)

# Defining FunctionSpace
## Voltage
V = fem.FunctionSpace(domain, ("Lagrange", 1))

## Particle density
P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
element = ufl.MixedElement([P1, P1, P1])
N = fem.FunctionSpace(domain, element)

## Particle velocity
v_cg = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
element = ufl.MixedElement([v_cg, v_cg, v_cg])
W = fem.FunctionSpace(domain, element)


# Define trailfunction and testfunction
## Voltage
u_voltage = ufl.TrialFunction(V)
v_voltage = ufl.TestFunction(V)

## Particle density
v_1, v_2, v_3 = ufl.TestFunctions(N)
u_1, u_2, u_3 = ufl.TrialFunctions(N)
u_n = fem.Function(N)
u_n1, u_n2, u_n3 = ufl.split(u_n)

## Particle velocity
velocity = fem.Function(W)
ve1, ve2, ve3 = ufl.split(velocity)


# Create boundary condition
## anode_marker
def anode_marker(x):
    return numpy.isclose(x[1], 5)

anode_boundary_entity = mesh.locate_entities_boundary(domain, 1, anode_marker)
anode_boundary_dof = fem.locate_dofs_topological(V, 1, anode_boundary_entity)
bc1 = fem.dirichletbc(PETSc.ScalarType(9000), anode_boundary_dof, V)

## cathode_marker
def cathode_marker(x):
    return numpy.isclose(x[1], 0)

cathode_boundary_entity = mesh.locate_entities_boundary(domain, 1, cathode_marker)
cathode_boundary_dof = fem.locate_dofs_topological(V, 1, cathode_boundary_entity)
bc2 = fem.dirichletbc(PETSc.ScalarType(0), cathode_boundary_dof, V)


# Create initial condition of partical density
# compute electric field intensity at the initial moment
def gaussian_distribution(x, y, x0, y0, sigma_x, sigma_y, amplitude):
    return amplitude * numpy.exp(-((x - x0)**2 / (2 * sigma_x**2) + (y - y0)**2 / (2 * sigma_y**2)))

def electron_initial_condition(x):
    return gaussian_distribution(x[0], x[1], 0, 4, 0.1, 0.1, 10**30)

u_n.sub(0).interpolate(electron_initial_condition)
u_n.sub(1).interpolate(electron_initial_condition)

# Create variational formulation
air_dielectric = fem.Constant(domain, PETSc.ScalarType(1))
elementary_charge = fem.Constant(domain, PETSc.ScalarType(1.6 * 10**(-19)))

a = dot(grad(u_voltage), grad(v_voltage)) * dx
L = elementary_charge / air_dielectric * (u_n2-u_n1-u_n3) * v_voltage * dx
problem = LinearProblem(a, L, bcs=[bc1, bc2], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
u_voltage = problem.solve()


# set simulation time, step size, and some parameters
T = 5e-6 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size
D1 = 0.18  # diffusion coefficient
D2 = 0.028e-4 # diffusion coefficient
D3 = 0.043e-4 # diffusion coefficient

t = 0
import tqdm.autonotebook
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
for i in range(num_steps):
    t += dt
    progress.update(1)
    # compute electric field intensity
    E = -grad(u_voltage)

    # compute particle mobility
    μ_e = 1.9163 * sqrt(dot(E, E)) ** (-0.25)
    μ_p = 2.43 * 1e-4
    μ_n = 2.7 * 1e-4

    # compute particle velocity
    v_e = - E * μ_e
    v_p = E * μ_p
    v_n = - E * μ_n
    
    # compute ionization and attachment coefficient
    alpha = 3.5 * 1e5 * exp(-1.65 * 1e7 / sqrt(dot(E, E)))
    eta = 1.5 * 1e3 * exp(-2.5 * 1e6 / sqrt(dot(E, E)))

    # create variational formulation
    def diffusion_advection(u, u_n, ve, v, D):
        return ((u - u_n) / dt)*v*dx + dot(ve, grad(u))*v*dx - D*dot(grad(u), grad(v))*dx
    
    ## redefine trialfunctions
    u_1, u_2, u_3 = ufl.TrialFunctions(N)
    F1 = diffusion_advection(u_1, u_n1, v_e, v_1, D1) - (alpha - eta) * u_1 * sqrt(dot(μ_e, μ_e)) * v_1 *dx
    F2 = diffusion_advection(u_2, u_n2, v_p, v_2, D2) - alpha * u_1 * sqrt(dot(μ_p, μ_p)) * v_1 *dx
    F3 = diffusion_advection(u_3, u_n3, v_n, v_3, D3) + eta * u_1 * sqrt(dot(v_n, v_n)) * v_1 *dx
    F = F1 + F2 + F3
    a = lhs(F)
    L = rhs(F)
    problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

    u = problem.solve()
    u_n1, u_n2, u_n3 = ufl.split(u)

    u_voltage = ufl.TrialFunction(V)
    a = ufl.dot(ufl.grad(u_voltage), ufl.grad(v_voltage)) * ufl.dx
    L = elementary_charge / air_dielectric * (u_n2-u_n1-u_n3) * v_voltage * ufl.dx
    problem = LinearProblem(a, L, bcs=[bc1, bc2], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    u_voltage = problem.solve()
6

What is the difference between them like this？

Is this interpolation method for components correct？

And is the components called in this place correct?

7

I don’t have time to review the code in detail.

However, note that you should define all variational forms and linear problems outside the for loop as explained in Form compilation — FEniCS Tutorial @ Sorbonne

8

Okay, thank you for your reminder.

So what about these three questions?

9

There is a fundamental different between a trialfunction and a function.
A trialfunction is used to create bilinear forms, it does not have unknowns as in DOFs.
A function contains degrees of freedom either being unknown (to be determined by a solve) or being known a priori, for instance through interpolation and assignment.

Seems fine to me.

Not sure what you mean by this.

Please bear in mind that there are many posts on this forum, and some of your questions are already answered.

10

Thank you for your prompt response. I will do it more carefully.

11

Hello, dokken.
Later on, I seriously pondered over the problem and worked hard to find similar questions and answers in the community. However, I am still very confused about the usage scenarios of several variable named functions.
Because I saw in the tutorial books and guessed that they seem to be interchangeable with each other.

u_n = fem.Function(N)
u_n1, u_n2, u_n3 = ufl.split(u_n)

&&

u_n1, u_n2, u_n3 = ufl.TrialFunctions(N)

What is the essential difference between them？
Looking forward to your answer.

12

From my understanding, the “Function” is used when you need to solve a nonlinear problem. For example

problem = NonlinearProblem(F, u, bcs=[bc1, bc2])

here, u is a function. And inside the NonlinearProblem, derivatives of the coefficients of u will be calculated automatically

V = u.function_space
du = ufl.TrialFunction(V)
J = ufl.derivative(F, u, du)

Trial function is used when you solve a linear problem.

2 Likes
13

Thank you for your answer.
I know my thoughts are a bit redundant. Can I ask you another question？
When I solve a linear problem, can the “Function()” and “LinearProblem()” be used?

14

No, in fact if you use Function in a LinearSolver, error occurs. For example, the poisson equation example in the tutorial. If you make the following change:

import ufl
#u = ufl.TrialFunction(V)
u = dolfinx.fem.Function(V)

you will get some errors like this

RuntimeError: Cannot create sparsity pattern. Form is not a bilinear.
Exception ignored in: <function LinearProblem.__del__ at 0x7fb7d3f71620>
Traceback (most recent call last):
  File "/mnt/d/software_install/fenics/lib/python3.12/site-packages/dolfinx/fem/petsc.py", line 627, in __del__
    self._solver.destroy()
    ^^^^^^^^^^^^
AttributeError: 'LinearProblem' object has no attribute '_solver'

you can try it.

1 Like
15

Thank you for your help. Your statement is completely correct.
In fact, I also tried multivariable coupling functions:

# u_1, u_2, u_3 = ufl.TrialFunctions(V)
u = fem.Function(V)
u_1, u_2, u_3 = ufl.split(V)

It will result in exactly the same error.

    problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 589, in __init__
    self._A = create_matrix(self._a)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 139, in create_matrix
    return _cpp.fem.petsc.create_matrix(a._cpp_object)
RuntimeError: Cannot create sparsity pattern. Form is not a bilinear.
Exception ignored in: <function LinearProblem.__del__ at 0x7f9ef5b46710>
Traceback (most recent call last):
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py", line 627, in __del__
    self._solver.destroy()
AttributeError: 'LinearProblem' object has no attribute '_solver'
16

I’m sorry.
Although I have studied the link tutorial you provided, I still have no way to define all variational forms and linear problems outside the for loop. The difficulty I encountered is that in my coupled equation system, the coefficients and source terms are not expressions defined by time and space, but are coupled with other terms.
Here is the coupled equation system: (Mainly three diffusion-drift equations and one Poisson’s equation)

image
image804×499 37.7 KB

image

At each time step, I need to calculate the voltage based on the density of each particle, and then calculate the mobility and ionization coefficient at different positions based on the negative gradient of the voltage, which is the electric field strength. For example, the formula for calculating electron mobility is as follows:

image

Then, the particle density is calculated through the diffusion-drift equations of three types of particles.

Therefore, my approach is to define the function space of voltage and particle density, while calculating the mobility and ionization coefficient using functions and UFL. For example, calculating the electron mobility is as follows:

# u_voltage is the calculated voltage function
E = -grad(u_voltage)
μ_e = 1.9163 * sqrt(dot(E, E)) ** (-0.25)

Combine them with functions to form variational formulas, repeating this process for each step.

In fact, I want to know if my approach is feasible.
And, is there a way to define all variational forms and linear problems in advance? If possible, how much resources can be saved approximately? I will consider it because I have many other jobs to attend to.