Simple exponential decay output wrong results

General
1

I’ve simplified the problem in my model to the model b_t=-b. More specifically, I’m using the nonlinear problem F = ((b - b_n) / dt) * v_1 * ufl.dx + b * v_1 * ufl.dx. When I test it I get exponential growth which is very weird. I would love to understand what I’m doing wrong.

My code:

import dolfinx as dx
import numpy as np
from mpi4py import MPI
import ufl
import pyvista as pv
import pyvistaqt as pvqt
dx.log.set_log_level(dx.log.LogLevel.ERROR)

# Simulation Constants
T = 100.0
dt = 1e-3
num_steps = T / dt


domain = dx.mesh.create_rectangle(MPI.COMM_WORLD, np.array(
    [[0, 0], [5, 5]]), [128, 128])

P1 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
element = ufl.MixedElement([P1])
V = dx.fem.FunctionSpace(domain, element)
v_1, = ufl.TestFunction(V)

u = dx.fem.Function(V)
u_n = dx.fem.Function(V)

b, = ufl.split(u)
b_n, = ufl.split(u_n)
b_n += 0.7

F = ((b - b_n) / dt) * v_1 * ufl.dx + b * v_1 * ufl.dx

problem = dx.fem.petsc.NonlinearProblem(F, u)
solver = dx.nls.petsc.NewtonSolver(domain.comm, problem)


# Plotting section - NOT INTERESTING!
# Create a VTK 'mesh' with 'nodes' at the function dofs
V0, dofs = V.sub(0).collapse()
topology, cell_types, x = dx.plot.create_vtk_mesh(V0)
grid = pv.UnstructuredGrid(topology, cell_types, x)

# Set output data
grid.point_data["b"] = u_n.sub(0).x.array[dofs].real
grid.set_active_scalars("b")

p = pvqt.BackgroundPlotter(title="biomass", auto_update=True)
p.add_mesh(grid, clim=[0, 100])
p.view_xy(True)
t = 0
p.add_text(f"time: {t:.2e}", font_size=12, name="timelabel")


def plot_state(t, u):
    p.add_text(f"time: {t:.2e}", font_size=12, name="timelabel")
    grid.point_data["b"] = u.sub(0).x.array[dofs].real
    p.app.processEvents()
# END of plotting section


for i in range(int(num_steps)):
    t += dt

    its, converged = solver.solve(u)
    assert converged
    if i % 100 == 0:
        print(f"{i}|\tt: {t}, iteration: {its}")
        plot_state(t, u)

    u_n.x.array[:] = u.x.array[:]
2

@Yotam_Ohad,

The issue comes from the line assigning directly to b_n += 0.7. You are messing with the Index b_n which is causing your exponential growth to 700 (0.7 / 1e-3). If you change the time step you should see the solution u scale to (0.7/dt). Instead you should initialize u_n with something like

u_n.sub(0).x.array[:] = np.ones(u_n.sub(0).x.array.size) * 0.7

-Sven

1 Like