Hello,
I am solving a system of PDEs from here for a cylinder.
I have changed the code a little bit and for example, eliminated the advection terms and I have added another term myself ( last term in the weak form in the code below).
Now I have a problem solving this weak form when the power of u_1 is not an integer in the last term ww*(betaa*u_1**2.01)*v_1*dx. Here in the code when I use 2 as the power it works perfectly, but when I use for example 2.01 or 1.4 or even 2.0001 it gave me the error like:
*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
This is a sample code to run and if you can solve my problem I will be appreciated.
from dolfin import *
import numpy as np
from numpy import linalg as LA
import time
import sys
import argparse
import glob
import re
import math
from mshr import *
from scipy.integrate import ode
# import matplotlib.pyplot as plt
# import ODE
# import param
rank=MPI.rank(MPI.comm_world)
#---------------------------------
parser = argparse.ArgumentParser()
parser.add_argument("-pd", "--polynomial_degree",
help="polynomial degree of the ansatz functions",
default=1,
type=int)
args = parser.parse_args()
set_log_level(30)
# set_log_level(16)
#---------------------------------
# Optimization options for the form compiler
parameters['form_compiler']['cpp_optimize'] = True
parameters['form_compiler']['representation'] = "uflacs"
parameters['form_compiler']['quadrature_degree'] = 4
parameters['form_compiler']['optimize'] = True
ffc_options = {'optimize' : True,
'eliminate_zeros' : True,
'precompute_basis_const' : True,
'precompute_ip_const' : True}
x0=0
y0 = 0
z0 = -5
#End mesh values
x1 = .25
y1 = .25
z1 = 5
n_x = 20
n_y = 20
n_z = 50
geometry1 = Cylinder(Point(x0, y0, z0),Point(x0,y0,z1), .65,.65)
geo2=Cylinder(Point(x0, y0, z0),Point(x0,y0,z1),0.5,0.5)
geometry=geometry1-geo2
bmesh = generate_mesh(geometry, 120)
# if not Monsoon:
# File("cylinder.pvd") << bmesh
# bmesh = BoxMesh(Point(x0, y0, z0),Point(x1, y1, z1), n_x, n_y, n_z)
# left_x = CompiledSubDomain("x[0]**2+x[1]**2)**0.5-0.8<1e-6 && on_boundary")
left_x = CompiledSubDomain("abs(pow(x[0],2)+pow(x[1],2)-0.25)<1e-1 && on_boundary")
right_x = CompiledSubDomain("abs(pow(x[0],2)+pow(x[1],2)-.4225)<1e-1 && on_boundary")
top_z = CompiledSubDomain("near(x[2],-5) && on_boundary")
bottom_z = CompiledSubDomain("near(x[2],5) && on_boundary")
boundary_indices = {'left_x': 1,
'right_x': 2,
'top_z': 3,
'bottom_z': 4,
}
boundary_markers = MeshFunction("size_t", bmesh, dim=2, value=0)
left_x.mark(boundary_markers, boundary_indices["left_x"])
right_x.mark(boundary_markers, boundary_indices["right_x"])
top_z.mark(boundary_markers, boundary_indices["top_z"])
bottom_z.mark(boundary_markers, boundary_indices["bottom_z"])
# if not Monsoon:
# file_BC = File('facet_domains_main.pvd')
# file_BC << boundary_markers
V0=FunctionSpace(bmesh,'CG',1)
V00=FunctionSpace(bmesh,'DG',0)
element_numbers = bmesh.num_cells()
Vbb=FunctionSpace(bmesh,'CG',1)
GG = FunctionSpace(bmesh, 'P', 1)
node_numbers = bmesh.num_vertices()
element_numbers = bmesh.num_cells()
connectivity_matrix = bmesh.cells()
local_range = Vbb.dofmap().ownership_range()
local_dim = local_range[1] - local_range[0]
# print (normal_normalized)
vert=vertices(bmesh)
T = 2 # final time
num_steps = 50 # number of time steps
dt = T / num_steps # time step size
eps = 0.05 # diffusion coefficient
K = 5 # reaction rate
element = FiniteElement('CG', bmesh.ufl_cell(),1)
V = FunctionSpace(bmesh, MixedElement([element,element,element]))
(v_1, v_2, v_3 )= TestFunctions(V)
# Define test functions
u_D = Expression(('1*exp(-abs(pow(x[2],3)/(.3)))','0.5','0'),degree=3)
u = Function(V)
u_n = Function(V)
# Split system functions to access components
# growth_rate.vector().set_local(growth_rate_numpy)
u_n = project(u_D, V)
u_1, u_2, u_3 = split(u)
u_n1, u_n2, u_n3 = split(u_n)
# u_n1.vector().set_local(growth_rate_numpy)
# sys.exit()
# u_n1.assign(u_D)
ww=1
EC50=0.5
nn=1.4
betaa = ((EC50**nn)-1)/(2*EC50**nn-1)
KK = (betaa-1)**(1/nn)
# Define source terms
f_1 = Constant(0.0)
f_2 = Constant(0.0)
f_3 = Constant(0.0)
# Define expressions used in variational forms
k = Constant(dt)
K = Constant(K)
eps = Constant(eps)
# Define variational problem
F = ((u_1 - u_n1) / k)*v_1*dx\
+ eps*dot(grad(u_1), grad(v_1))*dx +K*u_1*u_2*v_1*dx \
+ ((u_2 - u_n2) / k)*v_2*dx \
+ eps*dot(grad(u_2), grad(v_2))*dx + K*u_1*u_2*v_2*dx \
+ ((u_3 - u_n3) / k)*v_3*dx \
+ eps*dot(grad(u_3), grad(v_3))*dx - K*u_1*u_2*v_3*dx\
- f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx-ww*(betaa*u_1**2.01)*v_1*dx
t = 0
xdmffile_u = XDMFFile( 'displacement_PDE9'+'.xdmf')
count=0
for n in range(num_steps):
solve(F == 0, u)
# Save solution to file (VTK)
_u_1, _u_2, _u_3 = u.split()
# for i in range(local_dim):
# np_u3[i]=_u_3.vector().get_local()[3*i+2]
_u_1.rename('A', 'Label')
_u_2.rename('B', 'Label')
_u_3.rename('C', 'Label')
# vtkfile_u_1 << (_u_1, t)
# vtkfile_u_2 << (_u_2, t)
# vtkfile_u_3 << (_u_3, t)
if count%5==0:
xdmffile_u.write(_u_1, t)
xdmffile_u.write(_u_2, t)
xdmffile_u.write(_u_3, t)
count+=1
t += dt
u_n.assign(u)
if rank ==0:
print(t)
thank you very much for your help