Fixing instability for a simple sedimentation PDE

General
1

I am trying to use SFEpy to simulate a simple 1 D sedimentation problem which I do have an analytical solution for. I am getting a lot of instability in my solutions. The part of the solution in the lower radial range looks OK but near the bottom of the range there is so much instability. I tried implement a crank-nicolson scheme but it didn’t make a difference so I left it out here to keep the code simple. Any help I could get on this would be much appreciated!

Here is my example solution. The analytical solution would look like a step function for each timepoint.

image
image1164×882 89 KB

from mpi4py import MPI
import numpy as np
from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
import ufl
from petsc4py import PETSc
import matplotlib.pyplot as plt

# Parameters
s_svedberg = 20  # Sedimentation coefficient (s)
s_s = s_svedberg * 10 ** (-13)  # Convert to seconds
rpm = 30000  # Rotations per minute
omega_rad_per_s = 2 * np.pi * rpm / 60  # Angular velocity (rad/s)
t_s = 0.0 # start time
T_s = 1500.0  # Total simulation time (s)
num_steps = 1500  # Number of time steps
dt_s = T_s / num_steps  # Time step size (s)

# Define the radial positions
r_min_cm = 6.2  # Minimum radial position (cm)
r_max_cm = 7.2  # Maximum radial position (cm)

# Create mesh and define function space
domain = mesh.create_interval(MPI.COMM_WORLD, nx=100, points=[r_min_cm, r_max_cm])
V = fem.functionspace(domain, ("Lagrange", 1))

# Concentration at the previous time step
u_n = fem.Function(V) 
u_n.name = "u_n"
u_n.x.array[:] = 1 # initial concentration

# Define solution variable
c = fem.Function(V)  # Concentration at the current time step

# Define variational problem
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
r = ufl.SpatialCoordinate(domain)[0]  # Radial coordinate

# Strong form of the equation: dc/dt = -s * omega^2 * (r * dc/dr + 2c)
a = u *v* ufl.dx - dt_s * s_s * omega_rad_per_s**2 * (2 * u * v * ufl.dx + u * r * ufl.grad(v)[0] * ufl.dx)
L = u_n * v * ufl.dx  

# ## Preparing linear algebra structures for time dependent problems
bilinear_form = fem.form(a)
linear_form = fem.form(L)

# Create PETSc matrix and vector for the system
A = assemble_matrix(bilinear_form)
A.assemble()
b = create_vector(linear_form)

# Create a PETSc solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# Store solutions at 20 evenly spaced time points
time_points = np.linspace(0, T_s, 20)
solutions = []

# Time-stepping loop
for n in range(num_steps):
    t_s += dt_s

    # Assemble the right-hand side vector
    with b.localForm() as b_local:
        b_local.set(0.0)
    assemble_vector(b, linear_form)

    # Solve the linear system
    solver.solve(b, c.x.petsc_vec)
    c.x.scatter_forward()

    # Update previous solution
    u_n.x.array[:] = c.x.array.copy()

    # Save solution at the specified time points
    if n % (num_steps // 20) == 0:
        solutions.append((t_s, c.x.array.copy()))

# Plot solution using Matplotlib
if domain.comm.rank == 0:
    r_coords_cm = domain.geometry.x[:, 0]  # Radial coordinates (x-axis, in cm)

    plt.figure()
    for t_s, c_values in solutions:
        plt.plot(r_coords_cm, c_values, label=f"t = {t_s:.0f} s")
        plt.ylim(-0.01, 10.1)

    plt.xlabel("Radial Position (cm)")
    plt.ylabel("Concentration")
    plt.legend()
    plt.grid(True)
    plt.show()
2

Could you write out the mathematical formulation behind your problem?

3

The PDE is: dc/dt = -s * omega^2 * (r * dc/dr + 2c). It represents the change in concentration with time for solute under centrifugal force. The initial condition is that c = 1 for all radial positions.

4

So that is an advection-reaction equation. A few things come to mind:

In doing integration by parts to obtain this term, it looks like you have a sign mistake, and are missing the grad(r) component.

What boundary conditions are you intending to specify? Due to the integration by parts, it looks like you are specifying homogeneous Dirichlet conditions on both sides. The stability of advection equations depends on prescribing the correct boundary conditions (e.g., for advection diffusion you’re not allowed to prescribe Neumann at the inflow boundary). I’m not exactly sure what is permitted for advection reaction.

In general, advection reaction equations are quite tricky from a stabilization perspective in FEM. I’d read up on “advection reaction stabilization in FEM”. You’ll find, for example SUPG.