PETSc error code 63 for simple mixed element 1D simulation

variational formulation
1

Hello,

The simple Python code below shows my attempt to evolve two independent, uncoupled 1D equations. In the full problem, they are coupled, but this simple case demonstrates my difficulty.

The first equation is a envelope that propagates from left to right with constant velocity. The second equation is an initial condition that remains static. I can solve the single-equation cases with no problem, but when I try to solve them simultaneously, I get the following error message:

screnshot_fenics_error_63
screnshot_fenics_error_631208×721 82.2 KB

Here is the code:

Primary library for the FEniCS Python API (should not import *)

from dolfin import *

For plotting…

import matplotlib.pyplot as plt

Specify the BCs

def on_left(x, on_boundary):
return (on_boundary and near(x[0], 0.))

def on_right(x, on_boundary):
return (on_boundary and near (x[0], 1.))

V_x = 0.7212389380530974
n_ds = 226
ds = 0.004424778761061947
c_width = 0.4424778761061947
c_density = 1.e8
lp_width = 0.27876106194690264
lp_density = 1.

create a 1D mesh on the interval [0,1]

mesh = UnitIntervalMesh(n_ds)
V = VectorFunctionSpace(mesh, “Lagrange”, 1, dim=2)
v_1, v_2 = TestFunctions(V)

u = Function(V)
u_1, u_2 = split(u)
uprev = Function(V)
uprev_1, uprev_2 = split(uprev)

instantiate the corresponding FEniCS objects for Dirichlet BCs

bc_left = DirichletBC(V, [Constant(0),Constant(0)], on_left)
bc_right = DirichletBC(V, [Constant(0),Constant(0)], on_right)
bc = [bc_left, bc_right]

specify and apply the initial conditions

u0 = Expression((‘x[0]>2.ds && x[0]<(lp_width-2.ds) ?
lp_densitysin(pi(x[0]-2.ds)/(lp_width-4.ds)) : 0.',
'x[0]>=lp_width && x[0]<=(lp_width+c_width) ?
c_densitysin(pi(x[0]-lp_width)/c_width) : 0.’),
degree=2, lp_width=0.2, c_width=0.6, lp_density=1., c_density=1., ds=0.01)
u0.c_width = c_width
u0.lp_width = lp_width
u0.lp_density = lp_density
u0.c_density = c_density
u0.ds = ds

project the above expression onto the solution vector u

u = interpolate(u0, V)
uprev.assign(u)

n_steps = 100
time = 0.8
dt = time / (n_steps + 1)

DT = Constant(dt)
VX = Constant(V_x)

uplot = project(u, V)
uplot_1, uplot_2 = split(uplot)

plot some results

plot(uplot_1, title=(‘Photon density in laser pulse (initial)’))
plt.grid(True)
plt.show()
plt.close()

plot(uplot_2, title=(‘Density of crystal excited states (initial)’))
plt.grid(True)
plt.show()
plt.close()

for i_loop in range(0, n_steps):
F = ( (u_1-uprev_1)v_1/DT + VXu_1.dx(0)*v_1 ) * dx
+ ( (u_2-uprev_2)*v_2/DT ) * dx

solve(F==0, u, bc)
uprev.assign(u)

uplot = project(u, V)
uplot_1, uplot_2 = split(uplot)

plot some results

plot(uplot_1, title=(‘Photon density in laser pulse (final)’))
plt.grid(True)
plt.show()
plt.close()

plot(uplot_2, title=(‘Density of crystal excited states (final)’))
plt.grid(True)
plt.show()
plt.close()

2

Here you are redefining u, destroying its connection with

If you move the split after the last definition of u, your code should work.

For your next question, please make sure the code is actually possible to run after a copy paste, as there are several formatting issues in your current code block. Also remove all unnecessary plotting.

1 Like
3

Yes, that fixed the error in my code.
Thank you very much @dokken!

In the future, I will remove the plotting from my examples and will be more careful about how I paste it into the web form.

Sincerely,
David

1 Like