Hi, im working in a 3d linear elasticity model for a underground mining. In the model, i am using a anisotropic material, whose elastic constants are:
G=5e3
V1=0.3
V2=0.27
E1=30e3
E2=20e3
Where G and E are in MPa. For error, in some point i use G=5, and the model run correctly, but, when i fix this error for G=5e3, the model show an error about the solver:
File "verde_ani.py", line 176, in <module>
solver.solve(u.vector(), b);
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 77 (Petsc has generated inconsistent data).
*** Where: This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics_1555852928859/work/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2018.1.0
*** Git changeset: 9e1777c0df7c762c7230189c87c13fd5ad33a4f0
I think that the magnitude order for the elastic constants makes the constitutive equations tensor so heavy, but i dont know how fix it and havent found any way to fix it. The script is:
from dolfin import*
from mshr import*
from fenics import*
from ufl import nabla_div
import numpy as np
# Set backend to PETSC
parameters["linear_algebra_backend"] = "PETSc"
def build_nullspace(V, x):
"""Function to build null space for 3D elasticity"""
# Create list of vectors for null space
nullspace_basis = [x.copy() for i in range(6)]
# Build translational null space basis
V.sub(0).dofmap().set(nullspace_basis[0], 1.0);
V.sub(1).dofmap().set(nullspace_basis[1], 1.0);
V.sub(2).dofmap().set(nullspace_basis[2], 1.0);
# Build rotational null space basis
V.sub(0).set_x(nullspace_basis[3], -1.0, 1);
V.sub(1).set_x(nullspace_basis[3], 1.0, 0);
V.sub(0).set_x(nullspace_basis[4], 1.0, 2);
V.sub(2).set_x(nullspace_basis[4], -1.0, 0);
V.sub(2).set_x(nullspace_basis[5], 1.0, 1);
V.sub(1).set_x(nullspace_basis[5], -1.0, 2);
for x in nullspace_basis:
x.apply("insert")
# Create vector space basis and orthogonalize
basis = VectorSpaceBasis(nullspace_basis)
basis.orthonormalize()
return basis
mesh = Mesh("f3_nomanto.xml")
rho = Constant(0.027) #Densidad
alfa=90
a=alfa*np.pi/180
G=5e3
V1=0.3
V2=0.27
E1=30e3
E2=20e3
a11 = 1/E1
a12 = -((V1/E1)*np.cos(a)**2+(V2/E2)*np.sin(a)**2)
a13 = -((V1/E1)*np.sin(a)**2+(V2/E2)*np.cos(a)**2)
a14 = 0
a15 = 2*((V2/E2)-(V1/E1))
a16 = 0
a21 = -((V1/E1)*np.cos(a)**2+(V2/E2)*np.sin(a)**2)
a22 = (1/E1)*np.cos(a)**4 + ((1/G)-2*(V2/E2))*np.sin(a)**2*np.cos(a)**2 + (1/E2)*np.sin(a)**4
a23 = ((1/E1) + (1/E2) + 2*(V2/E2) - (1/G))*np.sin(a)**2*np.cos(a)**2 - (V2/E2)
a24 = 0
a25 = 2*((1/E1)*np.cos(a)**2-(1/E2)*np.sin(a)**2)-((1/G)-2*(V2/E2))*(np.cos(a)**2-np.sin(a)**2)*np.sin(a)*np.cos(a)
a26 = 0
a31 = -((V1/E1)*np.sin(a)**2+(V2/E2)*np.cos(a)**2)
a32 = ((1/E1) + (1/E2) + 2*(V2/E2) - (1/G))*np.sin(a)**2*np.cos(a)**2 - (V2/E2)
a33 = (1/E1)*np.sin(a)**4 + ((1/G)-2*(V2/E2))*np.sin(a)**2*np.cos(a)**2 + (1/E2)*np.cos(a)**4
a34 = 0
a35 = 2*((1/E1)*np.sin(a)**2-(1/E2)*np.cos(a)**2)+((1/G)-2*(V2/E2))*(np.cos(a)**2-np.sin(a)**2)*np.sin(a)*np.cos(a)
a36 = 0
a41 = 0
a42 = 0
a43 = 0
a44 = (1/E1)*2*(1+V1)*np.cos(a)**2 + (1/G)*np.sin(a)**2
a45 = 0
a46 = ((1/E1)*2*(1+V1) - (1/G))*np.sin(a)*np.cos(a)
a51 = 2*((V2/E2)-(V1/E1))
a52 = 2*((1/E1)*np.cos(a)**2-(1/E2)*np.sin(a)**2)-((1/G)-2*(V2/E2))*(np.cos(a)**2-np.sin(a)**2)*np.sin(a)*np.cos(a)
a53 = 2*((1/E1)*np.sin(a)**2-(1/E2)*np.cos(a)**2)+((1/G)-2*(V2/E2))*(np.cos(a)**2-np.sin(a)**2)*np.sin(a)*np.cos(a)
a54 = 0
a55 = 4*((1/E1) + (1/E2) + 2*(V2/E2) - (1/G))*np.sin(a)**2*np.cos(a)**2 + (1/G)
a56 = 0
a61 = 0
a62 = 0
a63 = 0
a64 = ((1/E1)*2*(1+V1) - (1/G))*np.sin(a)*np.cos(a)
a65 = 0
a66 = (1/E1)*2*(1+V1)*np.sin(a)**2 + (1/G)*np.cos(a)**2
S = np.array([[a11,a12,a13,a14,a15,a16],[a21,a22,a23,a24,a25,a26],[a31,a32,a33,a34,a35,a36],[a41,a42,a43,a44,a45,a46],[a51,a52,a53,a54,a55,a56],[a61,a62,a63,a64,a65,a66]])
C = np.linalg.inv(S)
D = as_tensor(C)
def eps(v):
return sym(grad(v))
def strain2voigt(e):
#"""e is a 2nd-order tensor, returns its Voigt vectorial representation"""
return as_tensor([e[0,0],e[1,1],e[2,2],2*e[1,2],2*e[0,2],2*e[0,1]])
def voigt2stress(vec):
#"""
#s is a stress-like vector (no 2 factor on last component)
#returns its tensorial representation
#""" sxx sxy sxz syx syy syz szx szy szz
return as_tensor([[vec[0], vec[5], vec[4]],
[vec[5], vec[1], vec[3]],
[vec[4], vec[3], vec[2]]])
def sigma(v):
return voigt2stress(dot(D,strain2voigt(eps(v))))
b_z = -rho
b = Constant((0.0,0.0,b_z))
#Function Spaces
V = VectorFunctionSpace(mesh, "Lagrange", 1)
u_tr = TrialFunction(V)
u_test = TestFunction(V)
a = inner(sigma(u_tr), eps(u_test))*dx #bilineal form
l = inner(b,u_test)*dx
x_max = 1216.0
y_max = 929.0
tol=1e-8
def bottom(x, on_boundary):
return (on_boundary and near(x[2],0,tol))
bc1 =DirichletBC(V.sub(2), Constant(0.),bottom)
def left(x, on_boundary):
return (on_boundary and near(x[0], 0,tol))
bc2 =DirichletBC(V.sub(0), Constant(0.),left)
def right(x, on_boundary):
return (on_boundary and near(x[0], x_max,tol))
bc3 =DirichletBC(V.sub(0), Constant(0.),right)
def north(x, on_boundary):
return (on_boundary and near(x[1], y_max,tol))
bc4 =DirichletBC(V.sub(1), Constant(0.),north)
def south(x, on_boundary):
return (on_boundary and near(x[1], 0,tol))
bc5 =DirichletBC(V.sub(1), Constant(0.),south)
bc=[bc1,bc2,bc3,bc4,bc5]
A, b = assemble_system(a, l, bc)
u = Function(V,name="Displacement")
# Create near null space basis (required for smoothed aggregation
# AMG). The solution vector is passed so that it can be copied to
# generate compatible vectors for the nullspace.
null_space = build_nullspace(V, u.vector())
# Attach near nullspace to matrix
as_backend_type(A).set_near_nullspace(null_space)
# Create PETSC smoothed aggregation AMG preconditioner and attach near
# null space
pc = PETScPreconditioner("petsc_amg")
# Use Chebyshev smoothing for multigrid
PETScOptions.set("mg_levels_ksp_type", "chebyshev")
PETScOptions.set("mg_levels_pc_type", "jacobi")
# Improve estimate of eigenvalues for Chebyshev smoothing
PETScOptions.set("mg_levels_esteig_ksp_type", "cg")
PETScOptions.set("mg_levels_ksp_chebyshev_esteig_steps", 50)
# Create CG Krylov solver and turn convergence monitoring on
solver = PETScKrylovSolver("cg", pc)
solver.parameters["monitor_convergence"] = True
# Set matrix operator
solver.set_operator(A);
# Compute solution
solver.solve(u.vector(), b);
Vsig = TensorFunctionSpace(mesh, "Discontinuous Lagrange",0)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u),V=Vsig, solver_type="cg"));
file_results = XDMFFile("Verde_ANi2.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)
I really hope you can help me.