Nonlinear PDE PETSc Error 63

variational formulation
1

Good evening,

I am trying to solve a very simple 1D non linear PDE with constant Dirichlet boundary conditions (I’m a true beginner at FeniCS). After the first iteration, I always get the same error code:

Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /usr/local/miniconda/conda-bld/fenics-pkgs_1565954444213/work/dolfin/dolfin/la/PETScMatrix.cpp.

Solutions I’ve looked around the Internet involve split functions but I can’t really understand if those apply to my problem; by the way, I’ve installed FeniCS thanks to conda-forge. Anyway here is the code:

maillage = UnitIntervalMesh(50)
H_P1 = FunctionSpace(maillage, 'P', 1)
I = Function(H_P1)
phi = TestFunction(H_P1)

F =  (d/25) * div(grad(I)) * phi * dx + (beta/alpha) * phi * (I - (I**2)/K) * dx

class Dirichlet_0(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0,10**(-14)) and on_boundary

class Dirichlet_1(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],1,10**(-14)) and on_boundary
    
condition_0 = DirichletBC(H_P1, Constant(300), Dirichlet_0())
condition_1 = DirichletBC(H_P1, Constant(0), Dirichlet_1())

I = Function(H_P1)
J = derivative(F,I)

problem = NonlinearVariationalProblem(F, I, [condition_0, condition_1], J=J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

Could you help me fix that error?

Thank you in advance

2

You’ve redefined I which means you don’t have a bilinear operator when you compute the derivative J. So you’re assembling an empty LHS matrix.

Also, your FE formulation doesn’t look correct, but I’m assuming you can figure that out.

3

First of all, thank you very much for your quick response !
You’re completely right about the derivative J, I’ve moved that just after the definition of F.

On the topic of the variational formula, my PDE is -\frac{d}{L^2} \frac{\partial^{2}I_{\infty}}{\partial^2 x^{*}} = \frac{\beta}{\alpha}I_{\infty}(1-\frac{I_{\infty}}{K}) to be solved on ]0,1[. Have I misdefined F ?

4

That’s the so-called strong formulation. You need to find the weak formulation by integrating by parts the term with the second order derivative, and imposing appropriate boundary conditions.

For example
strong form: -\nabla^2 u = f, u = u_D on the Dirichlet boundary \partial\Omega_D and \nabla u \cdot n = g_N on the Neumann boundary \partial\Omega_N.

multiply by test function, v and integrate: -\int_\Omega \nabla^2 u v \mathrm{d}x = \int_\Omega f v \mathrm{d} x

and integrate by parts to get weak formulation: \int_\Omega \nabla u \cdot \nabla v \mathrm{d}x - \int_{\partial\Omega_N} g_N v \mathrm{d}s = \int_\Omega f v \mathrm{d} x

5

Let \phi \in \mathcal{H}^1(0,1) and I \in \mathcal{C}^2(0,1) solution of the PDE.
We have \frac{d}{L^2}(\int_{]0,1[} I' \phi'dx - \int_{\delta(]0,1[)} I' \phi dS) = \frac{\beta}{\alpha}\int_{]0,1[} I(1-\frac{I}{K}) \phi dx so

F = d/(L**2) * dot(grad(I),grad(phi)) * dx - grad(I)[0] * phi * ds - (b/a) * I * (1-I/K) * phi * dx

Do you agree with that expression?

If Neumann boundary conditions were to be respected (for example I'(0) = c and I'(1) = d), how would I define ds1 and ds2 in FeniCS to make the following code work?

F = d/(L**2) * dot(grad(I),grad(phi)) * dx - d * phi * ds1 + c * phi * ds2 - (b/a) * I * (1-I/K) * phi * dx