Pde with integral

General
1

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 :slight_smile:

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 !

2

I would introduce a scalar valued lagrange multiplier \lambda = \int_{\Omega} (...) dx. You can now do that in scifem GitHub - scientificcomputing/scifem: Scientific finite element toolbox by @dokken and @Henrik_Nicolay_Finsb

4

closing as per user request, but not deleting.