Hello everyone,
I am trying to solve the following PDE:
\partial_t\psi = -\left[ (1-\varepsilon)\psi +\psi^3+2\Delta \psi + \Delta^ 2 \psi - \frac{1}{\Omega} \int_\Omega ((1-\varepsilon)\psi +\psi^3)dx\right]
I am solving the problem under periodic boundary conditions, using a mixed element and a Crank–Nicolson scheme.
However, I am worried about the integral term. I think it is not updated in my non linear problem thoughout the time loop.
Is there a better way to approach such pde ?
The integral term should ensure that the average \psi over the domain is constant over time.
Here is my code
from dolfinx import fem, mesh, plot
from ufl import dx, grad, inner
import numpy as np
from mpi4py import MPI
import ufl
from utils import *
import ufl.integral
# import pytest
from MyMPCSolver import *
import time
import pyvista
#### Model parameter
dt=1e-2
eps=+0.3
avg = -0.3
################################Generate Mesh###################
xell_ = np.pi/np.sqrt(3)
yell_ = np.pi
L = 20*xell_; H = 20*yell_
cell_sizex = xell_/5
cell_sizey = yell_/5
nx = int(L/cell_sizex)
ny = int(H/cell_sizey)
comm = MPI.COMM_WORLD
domain = mesh.create_rectangle(comm, [(0.0, 0.0), (L, H)], [nx, ny],mesh.CellType.triangle)
ndim = domain.geometry.dim
############################## FInite element definition#####################""""
P = basix.ufl.element("Lagrange", domain.basix_cell(), 1)
space = fem.functionspace(domain, P)
ME = fem.functionspace(domain, basix.ufl.mixed_element([P, P]))
q, v = ufl.TestFunctions(ME)
Chi = fem.Function(ME)
Chi0 = fem.Function(ME)
psi, mu = ufl.split(Chi)
psi0, mu0 = ufl.split(Chi0)
#############################################Initial Condition
Chi.x.array[:] = 0.0
rng = np.random.default_rng(13)
initialCpsi = lambda x : (rng.random(x.shape[1])-0.5)*0.2-0.3
Chi.sub(0).interpolate(initialCpsi)
Chi.x.scatter_forward()
##################### Define non linear problem
psi_mid = (1/2)*(psi0 + psi)
mu_mid = (1/2)*(mu0 + mu)
Fpsi = inner(psi,q)*dx -inner(psi0,q)*dx+dt((1-eps)*inner(psi_mid,q)*dx
+(1/2)*inner(psi**3+psi0**3,q)*dx
-inner(grad(mu_mid),grad(q))*dx
+2*inner(mu_mid,q)*dx
-(1/(L*H))* fem.assemble_scalar(fem.form(((1-eps)*psi_mid+(1/2)*(psi**3+psi0**3))*dx))*
inner(1,q)*dx
)
Fmu = inner(mu,v)*dx+inner(grad(psi),grad(v))*dx
F=Fpsi+Fmu
###################### Boundary conditions and solver
mpc = periodic_bcs(domain, ME)
problem = NonlinearMPCProblem(F, Chi, mpc, bcs=[])
solver = NewtonSolverMPC(comm, problem, mpc)
##### Defined output
file = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "ConservedSH.xdmf", "w")
file.write_mesh(domain)
V0, dofs = ME.sub(0).collapse()
topology, cell_types, x = plot.vtk_mesh(V0)
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
######################################## Run main loop
psi = Chi.sub(0)
Chi0.x.array[:] = Chi.x.array
t1=time.time()
t=0
file.write_function(psi, t)
while t < 100:
t += dt
r = solver.solve(Chi)
print(f"Step {int(t / dt)}: num iterations: {r[0]}","time=",t)
Chi0.x.array[:] = Chi.x.array
file.write_function(psi, t)
t2=time.time()
print("Job done in ", t2-t1)
I really appreciate your help
Thank you very much !