Hi everyone,
I am trying to implement some code able to solve the following Schrodinger equation
Do you have some advice on how to treat the variational formulation, since both the test function and the wave function psi are complex? I have read that dolfinx should be able to handle this task.
Thanks in advance
Have you considered the demo: dolfinx/demo_helmholtz.py at main · FEniCS/dolfinx · GitHub
Yes, I have checked it, even though I keep getting the following error on this part of the code:
a = ( ufl.inner( psi1 , v) - j * hbar * dt / 2.0 * ufl.inner(ufl.grad(psi1), ufl.grad(v) ) ) * ufl.dx
L = ufl.inner( psi1_old, v) * ufl.dx
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, psi1, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()
The error is:
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
Input In [4], in <cell line: 8>()
6 from dolfinx.fem.petsc import LinearProblem
7 problem = LinearProblem(a, L, psi1, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
----> 8 problem.solve()
File /usr/local/dolfinx-complex/lib/python3.8/dist-packages/dolfinx/fem/petsc.py:560, in LinearProblem.solve(self)
558 # Assemble lhs
559 self._A.zeroEntries()
--> 560 assemble_matrix(self._A, self._a, bcs=self.bcs)
561 self._A.assemble()
563 # Assemble rhs
File /usr/lib/python3.9/functools.py:877, in singledispatch.<locals>.wrapper(*args, **kw)
873 if not args:
874 raise TypeError(f'{funcname} requires at least '
875 '1 positional argument')
--> 877 return dispatch(args[0].__class__)(*args, **kw)
File /usr/local/dolfinx-complex/lib/python3.8/dist-packages/dolfinx/fem/petsc.py:338, in _(A, a, bcs, diagonal, constants, coeffs)
336 constants = _pack_constants(a) if constants is None else constants
337 coeffs = _pack_coefficients(a) if coeffs is None else coeffs
--> 338 _cpp.fem.petsc.assemble_matrix(A, a, constants, coeffs, bcs)
339 if a.function_spaces[0] is a.function_spaces[1]:
340 A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
File /usr/local/lib/python3.9/dist-packages/ufl/core/expr.py:390, in Expr.__len__(self)
388 if len(s) == 1:
389 return s[0]
--> 390 raise NotImplementedError("Cannot take length of non-vector expression.")
NotImplementedError: Cannot take length of non-vector expression.
Could you supply the full code to reproduce the error message?
Yes, of course
# Geometry generation
import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from dolfinx import mesh, fem, plot, io
# Scaled variable
L = 20
W = 0.5
H = 8
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0,0,0]), np.array([L, W, H])],
[512,12,208], cell_type=mesh.CellType.hexahedron)
n = ufl.FacetNormal(domain)
from dolfinx import io
with io.XDMFFile(domain.comm, "mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
print('Mesh Generated')
# Parameters
dt = 1/48
T = 20
radius = 0.3
uD = np.array([1, 0, 0]);
hbar = 0.03
c = 0.01
xcylinder = 2
zcylinder = 4
V = fem.FunctionSpace(domain, ("Lagrange",1))
psi1 = ufl.TrialFunction(V)
psi2 = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Psi initialization
from petsc4py import PETSc
j = PETSc.ScalarType(0 + 1j)
psi1_old = fem.Function(V)
psi1_old.interpolate(lambda x: 1/(1+c**2) + x[0] * 0 + x[1] * 0 + x[2] * 0 + j * 0)
psi2_old = fem.Function(V)
psi2_old.interpolate(lambda x: c/(1+c**2) + x[0] * 0 + x[1] * 0 + x[2] + j * 0)
# Schrodinger Solving
a = ( ufl.inner(psi1, v) - j * hbar * dt / 2.0 * ufl.inner(ufl.grad(psi1), ufl.grad(v))) * ufl.dx
# a = ( ufl.inner(psi1, v) - j * hbar * dt / 2.0 * ufl.inner(ufl.grad(psi1), ufl.grad(v))) * ufl.dx
L = ufl.inner( psi1_old, v) * ufl.dx
from dolfinx.fem.petsc import LinearProblem
problem = LinearProblem(a, L, psi1, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()
I’ll post the code in case someone needs it
import numpy as np
from time import process_time
import ufl
from dolfinx import plot
import pyvista
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form, locate_dofs_geometrical, dirichletbc, Expression
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, CellType
from ufl import dx, grad, inner, div, conj, real, imag, Index, Dx
from mpi4py import MPI
from petsc4py import PETSc
if np.issubdtype(PETSc.ScalarType, np.complexfloating):
imaginaryUnit = PETSc.ScalarType(0 + 1j)
else:
error('Complex Treatment is needed')
# Parameters
dt = 1/50;
hbar = 0.03;
inlet_velocity = np.array([1.0, 0.0]);
T = 50;
numsteps = 50/dt
# Mesh Discretization
dimensions = np.array([10,4]);
nx = 250;
ny = 100;
domain = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), dimensions],
[nx, ny], CellType.triangle)
tdim = domain.topology.dim;
# Definition of the functional space
deg = 1; # linear = 1, parabolic = 2 Lagrangian FE are used
V = FunctionSpace(domain, ("Lagrange", deg))
psi1 = Function(V)
psi2 = Function(V)
psi1.name = "psi1"
psi2.name = "psi2"
tmpPsi1 = Function(V)
tmpPsi2 = Function(V)
psi1old = Function(V)
psi2old = Function(V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Defining the inlet boundary condition
dofs_inlet = locate_dofs_geometrical(V, lambda x: np.isclose(x[0],0))
inletPlaneWave = Function(V)
kvec = inlet_velocity/hbar
inletPlaneWave.interpolate(lambda x: np.exp(imaginaryUnit * (kvec[0]*x[0] + kvec[1]*x[1])))
bc_Inlet = dirichletbc(inletPlaneWave, dofs_inlet)
# Initialization
def initial_condition1(x):
return 1 + 0j + 0 * x[0] + 1 * x[1]
def initial_condition2(x, c = 0.01 ):
return c + 0j + 0 * x[0] + 1 * x[1]
psi1old.interpolate(initial_condition1)
psi2old.interpolate(initial_condition2)
# Definition of the variational formulation
theta = 0.5 # Crank-Nicolson is used - theta method is implemented
a = inner(u, v) * dx - imaginaryUnit * dt * hbar / 2 * (1-theta) * inner(grad(u), grad(v)) * dx
L1 = inner(psi1old, v) * dx + imaginaryUnit * dt * hbar / 2 * theta * inner(grad(psi1old), grad(v)) * dx
L2 = inner(psi2old, v) * dx + imaginaryUnit * dt * hbar / 2 * theta * inner(grad(psi2old), grad(v)) * dx
problem1 = LinearProblem(a, L1, u = tmpPsi1, bcs = [bc_Inlet], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem2 = LinearProblem(a, L2, u = tmpPsi2, bcs = [bc_Inlet], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
cputime_start = process_time()
n = 0;
for n in range (0,int(numsteps)):
t = float(n+1)*dt
print("Time = %.2f" % t)
# Schrodinger Evolution
problem1.solve()
problem2.solve()
psi1old.interpolate(psi1)
psi2old.interpolate(psi2)
The issue here is that you send in psi1, a trial function in the place where the boundary conditions are expected. You can simply change it to:
problem = LinearProblem(a, L, petsc_options={
"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()