Newton solver diverges
I'm struggling with mesh convergence of the problem described here: https://www.allanswered.com/post/xnnob/bc-for-a-hyperbolic-1st-order-system-of-pde/
I'm using the solution proposed by pf4d.
When I increase the number of elements in the mesh I would expect the solutions to converge. However, the min and max values fluctuate widely.
I have solved the problem for increasing mesh densitiy N and printed the min and max of the solution. This is the output:
N= 32 min = -14.213930736628392 max = 9.046530167293112 N= 64 min = -2.2292666922532374 max = 2.5227212758820046 N= 128 min = -499.57598020660737 max = 460.7371101576128 N= 256 min = -6.909612880745566 max = 6.33138530677423 N= 512 min = -5.53777342589658 max = 3.98577741119527 N= 1024 min = -30.5893553449149 max = 23.895255532377046 N= 2048 .... *** Error: Unable to successfully call PETSc function 'KSPSolve'. *** Reason: PETSc error code is: 76 (Error in external library). *** Where: This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
If I use econd order elements instead of first order elements,
the solutions clearly diverge:
P1 = FiniteElement('CG', triangle,2)
I guess that the error message comes because the mesh is too fine and /or numbers are getting too big.
N= 32 min = -2.5901446402149064 max = 1.6248883001122068 N= 64 min = -61.72313541884393 max = 98.2433484043985 N= 128 min = -133.96435807556034 max = 176.622721145695 N= 256 min = -229.42673119333665 max = 238.81759600560238 N= 512 min = -992.7646840741731 max = 968.4059946564718 N= 1024 ... *** Error: Unable to successfully call PETSc function 'KSPSolve'. *** Reason: PETSc error code is: 76 (Error in external library). *** Where: This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
Should I use another solver? Which one might give a better convergence? Any other proposals what I could do?
from fenics import * from mshr import generate_mesh, Rectangle # Define domain and create mesh H = 0.1 L = 1.0 domain = Rectangle(Point(0, 0), Point(L, H)) mesh = generate_mesh(domain, 400) # Define the function space P1 = FiniteElement('CG', triangle, 1) element = MixedElement([P1, P1]) V = FunctionSpace(mesh, element) # Define test and trial functions psi = TestFunction(V) phi = TrialFunction(V) u = Function(V) h = CellSize(mesh) f = Constant((0,0)) # Define boundaries def surface(x, on_boundary): return on_boundary and near(x, 0) def glacier(x, on_boundary): return on_boundary and near(x, L) def front(x, on_boundary): return on_boundary and near(x, 0) # Define front stress front_normal_stress = Expression('-x', degree=1) # Define boundary conditions bc_surface = DirichletBC(V, Constant((0, 0)), surface) bc_glacier = DirichletBC(V, Constant((0, 0)), glacier) bc_front = DirichletBC(V.sub(0), front_normal_stress, front) bc = [bc_glacier, bc_surface, bc_front] # SUPG intrinsic time parameter : tau = 10*h # the deviatoric stress tensor : def sigma(u): return as_matrix([[u, u], [u, -u]]) # the differential operator : def Lu(u): return div(sigma(u)) # Define the variational problem : F = inner(Lu(u), psi)*dx - inner(f, psi)*dx \ + inner(Lu(psi), tau*Lu(u)) * dx \ - inner(Lu(psi), tau*f) * dx # define the Jacobian : J = derivative(F, u, phi) # compute solution using Newton's method : solve(F==0, u, bc, J=J) print u.vector().min(), u.vector().max() File('u.pvd') << u
Note that I have adjusted "tau" to be dependant on the mesh size, and found by experimentation that a factor of 10 produces a smooth result ... You should experiment with this value. For mathematical background regarding the technique, look up "Galerkin/least-squares" stabilization, or read the paper:
Thomas J. R. HUGHES, Leopoldo E. FRANCA and Gregory M. HULBERT (1988): A NEW FINITE ELEMENT FORMULATION FOR COMPUTATIONAL FLUID DYNAMICS: VIII. GALERKIN/LEAST-SQUARES METHOD FOR ADVECTIVE-DIFFUSIVE EQUATIONS.
Usually, water pressure is applied over the terminus, but as you have specified hydrostatic internal pressure and eliminated the forcing term \(f\) from the momentum balance, you are specifying deviations from the mean pressure.