Restricting trial function to converge the solution

variational formulation
1

Hi everyone,
I am trying to solve a nonlinear problem in a mixed space. My solution doesn’t converge due to the presence of tan(ksi) in the integrand4. I want to know if we can somehow restrict ksi to ignore the values pi/2. If yes, how can we do this? e.g. if we can substitute ksi by arctan(z) or any other trigonometric function to make this work.
I’ll highly appreciate if you can find some time to look at my code and provide me with your valuable suggestion.
Here’s my code

import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
from ufl.operators import min_value
import scipy.integrate as integrate
import ufl as ufl
from dolfin import *
import matplotlib.pyplot as plt

mesh = IntervalMesh(2000, 0, 10)
phiP_element = fe.FiniteElement('CG', mesh.ufl_cell(), 2)
ksi_element  = fe.FiniteElement('CG', mesh.ufl_cell(), 2)
pk_element = fe.MixedElement([phiP_element, ksi_element])
pk_space = fe.FunctionSpace(mesh, pk_element) 
Real_scalar = fe.FunctionSpace(mesh, 'CG', 1)
Real_vector = fe.VectorFunctionSpace(mesh, 'CG', 1)
def left(x, on_boundary):
    return on_boundary and fe.near(x[0], 0.0) 

def right(x, on_boundary):
    return on_boundary and fe.near(x[0], 10.0)

bc_1 = fe.DirichletBC(pk_space.sub(0), fe.Constant(0.0), left) 
bc_2 = fe.DirichletBC(pk_space.sub(0), fe.Constant(0.0), right)
bc_3 = fe.DirichletBC(pk_space.sub(1), fe.Constant(0.0), left) 
bc_4 = fe.DirichletBC(pk_space.sub(1), fe.Constant(np.pi/2-0.00001), right)
bcs = [bc_1, bc_2, bc_3, bc_4] 
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for vertex in entities(mesh, 0):
    if (vertex.midpoint().x() < 5.0):
        marker[vertex] = 1
    elif (vertex.midpoint().x() 5.0): 
        marker[vertex] = 2
    else:
        marker[vertex] = 3
xa = Expression("x[0]",degree=1)
r=interpolate(xa,Real_scalar)
e_phiP = Expression('0.0', degree=1)
e_ksi = Expression(('atan((x[0])/(sqrt((r_0 * r_0)- (x[0]*x[0]))))'), degree=1, r_0  = 10) # expression
phiP_0 = interpolate(e_phiP, pk_space.sub(0).collapse())
ksi_0 = interpolate(e_ksi, pk_space.sub(1).collapse())
pk = fe.Function(pk_space)
assign(pk, [phiP_0, ksi_0])
pk = fe.Function(pk_space) 
phiP, ksi = fe.split(pk)
ksi=ksi_0
phiP=phiP_0
vn = fe.TestFunction(pk_space)
v1, v2 = fe.split(vn)
dx = Measure("dx", domain=mesh)
ds = Measure("dS", domain=mesh, subdomain_data=marker)

K_c = fe.Constant(1.0)
k_o = fe.Constant(0.0)
K_0 = fe.Constant(20.0)
gamma_L = fe.Constant(1.0)
delta = fe.Constant(4.0)
del_P = fe.Constant(1.0)
r_p = r + delta * sin(ksi)

integrand1 = (1 + (float(delta) *cos(ksi) *fe.Dx(ksi,0))) *2 * np.pi * r_p * (1/cos(ksi)) * ((K_c/2) * (((cos(phiP + ksi) + cos(phiP - ksi))/2) * fe.Dx(phiP + ksi,0) * ((1 + (delta * cos(ksi) * fe.Dx(ksi,0))) ** (-1)) + (sin(phiP + ksi)/r_p) - k_o) ** 2 - (k_o ** 2))
integrand2 = (K_0/2) * np.pi * r_p * (1 - cos(2 * phiP)) * (1/cos(ksi))* (1 + (float(delta) *cos(ksi) *fe.Dx(ksi,0)))
integrand3 = gamma_L * 2 * np.pi * r_p * (1/cos(ksi)) *(1 + (float(delta) *cos(ksi) *fe.Dx(ksi,0)))
integrand4 = del_P * fe.pi * (r_p * r_p) * tan(ksi) * (1 + (float(delta) *cos(ksi) *fe.Dx(ksi,0))) 

functional = integrand1 * fe.dx + integrand2 * fe.dx+ integrand3 * fe.dx  + integrand4 * fe.dx
Integ=functional

weak_form = fe.derivative(Integ , pk, fe.TestFunction(pk_space)) 
jacobian = fe.derivative(weak_form, pk, fe.TrialFunction(pk_space))
problem = fe.NonlinearVariationalProblem(weak_form, pk, bcs, jacobian)
solver = fe.NonlinearVariationalSolver(problem)
solver.solve()

and the error is:

*** Error: Unable to successfully call PETSc function ‘MatSetValuesLocal’.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside ./dolfin/la/PETScMatrix.cpp.
*** Process: 0

*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.571e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
2

This is not correct. Here you are overwriting your solution function by the initial guess. You should use a FunctionAssigner.

3

Dear Prof. Dokken,
Thanks alot for your response and for your valuable suggestion. I have overwritten these functions because the solver completely ignored them. I will be using a FunctionAssigner to invoke these values. If you have any suggestions for a better approach to solving this problem, I would greatly appreciate them.

4

Im not sure what you mean by this.

You can have a look at:

5

Using FunctionAssigner returns the following error:

TypeError Traceback (most recent call last)
/tmp/ipykernel_75412/1312915568.py in
41 ksi_0 = interpolate(e_ksi, pk_space.sub(1).collapse())
42 pk = fe.Function(pk_space)
—> 43 assigner = fe.FunctionAssigner(pk, [pk_space(0), pk_space(1)])
44 assigner.assign(pk, [phiP_0, ksi_0])
45 pk = fe.Function(pk_space)

TypeError: ‘FunctionSpace’ object is not callable

6

Nowhere in the post does it state that you can do pk_space(0), pk_space(1).
These should be fe.FunctionAssigner(pk_space, [phiP_0.function_space(), ksi_0.function_space()])

7

Thank you very much. It’s working great. However, the Newton solver failed to converge and returns the following error.

*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0

8

That is now a very different error, this means that the non-linear solver was not able to solve the problem with the set maximum iteration bound.

You can increase the maximum number of iterations, or change the underlying solver, to see if it converges.

I would strongly suggest you try to simplify your example, so that you can pinpoint what makes the problem fail.

9

I sincerely appreciate your valuable suggestions and I will try to implement them. I am optimistic that following your advice will lead to successful outcomes.