Boundary Condition Error in MethodOfLines

Specific Domains Modelling & Simulations
, ,
1

I’m trying to use MethodOfLines to solve the time-independent Schrodinger equation for hydrogen, as a prelude to more complicated systems. The equation is as follows:

\nabla^2 \psi(\vec{r}) - \frac{1}{r}\psi(\vec{r}) = E\psi(\vec{r})

This problem has analytically known solutions and eigenvalues. Following the tutorials for the steady-state heat equation and 1D Schrodinger equation in MethodOfLines.jl, I’ve written the following:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, NonlinearSolve

# Define parameters (indep. vars)
@parameters x, y, z, E 

# Define variables (known and unknown functions of params)
@variables psi(..)

# Define operators
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dzz = Differential(z)^2

V(x, y, z) = -1/sqrt(x^2 + y^2 + z^2)

# Define PDE
eq = [E*psi(x, y, z) ~ Dxx(psi(x, y, z)) + Dyy(psi(x, y, z)) + Dzz(psi(x, y, z)) + V(x, y, z)*psi(x, y, z)]

# set domain ranges and construct domain
xmin, ymin, zmin = 0, 0, 0
xmax, ymax, zmax = 10, 10, 10 

domains = [x in Interval(xmin, xmax), y in Interval(ymin, ymax), z in Interval(zmin, zmax)]

# set boundary conditions
bcs = [psi(0, 0, 0) ~ 0,
       psi(xmax, ymax, zmax) ~ 0]

# construct system, discretize
@named sys = PDESystem(eq, bcs, domains, [x, y, z], [psi(x,y,z)])
disc = MOLFiniteDifference([x => 100, y => 100, z => 100])
prob = discretize(sys, disc)

# solve
sol = NonlinearSolve.solve(prob, NewtonRaphson())

println(sol[E])

(Note: the operators and variables are specified in Cartesian coordinates, since the documentation says its restricted to Cartesian grids. The same page says it can handle spherical Laplacians but then says nothing else and provides no examples, so go figure.) As written, this currently errors on the boundary condition psi(0, 0, 0) ~ 0, with the following message:

ERROR: AssertionError: Boundary condition psi(0, 0, 0) ~ 0 is not on a boundary of the domain, or is not a valid boundary condition
Stacktrace:
 [1] parse_bcs(bcs::Vector{…}, v::PDEBase.VariableMap, orders::Dict{…})
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\parse_boundaries.jl:370
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; checks::Bool)
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\symbolic_discretize.jl:27
 [3] symbolic_discretize
   @ C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\symbolic_discretize.jl:9 [inlined]
 [4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, checks::Bool, kwargs::@Kwargs{})
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\discretization_state.jl:74
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\discretization_state.jl:69
 [6] top-level scope
   @ c:\Users\njsjr\Documents\attosecond-research\helium-tdse\_research\pde_test.jl:36
Some type information was truncated. Use `show(err)` to see complete types.

I’m very confused by this, because I set the boundaries of my domain to be 0 and 10, so I don’t see why there should be an issue here.

2

Bear with me. I hope to be useful even if it turns out that I’m misguided.

The immediate problem is that the MethodOfLines requires boundary conditions on faces, not vertices. psi(0,0,0) ~ 0 defines a point in 3D, not a surface.

But, once past that E is an unknown Eigenvector and as used you’ll just get ψ = 0 which is not very informative. What I think is needed is to take an initial whack at a trial value for E or add a normalization constraint and treat E as a free variable but be careful to avoid singularity at the origin by a slight shift. in one dimension, say xmin=0.1. Take care not to start with a large grid size or you may find yourself wondering if it will ever finish. Better yet explolit spherical symmetry and reduce to one dimension in r,

-u''(r) + V(r)·u(r) = E·u(r)

which avoids the 3D grid.

But please don’t follow this rabbit too far into Wonderland, because I haven’t tested it. Offered just as a potentially useful insight.