Feed values from data file to a interpolate function
I was wondering that, in the demo of solving the nonlinear problem using picard iteration method in picard_np.py, how to provide initial values for the interpolate funtion u_k (previous known solution of u in picard iterations) from a data file.
To be specific, I first run the original picard_np.py and output the converged solution of u to a txt file 'u_sol.txt'. Then, I add two more lines following the interpolation of u_k,
u_k = interpolate(Constant(0.0), V) # previous (known) u u_k.vector().array()[:] = np.loadtxt('u_sol.txt') # provide converged solution of u as u_k print(u_k.vector().array())
which intends to use the converged solution as a good initial guess of u_k. However, when I print out the array of u_k and find it didn't change (all are still zero as interpolated).
So my question is that, how to initialize the u_k function with specified values (which in this case is the converged solution of u) from a data file.
The complete script is as follows:
""" FEniCS tutorial demo program: Nonlinear Poisson equation with Dirichlet conditions in x-direction and homogeneous Neumann (symmetry) conditions in all other directions. The domain is the unit hypercube in of a given dimension. -div(q(u)*grad(u)) = 0, u = 0 at x=0, u=1 at x=1, du/dn=0 at all other boundaries. q(u) = (1+u)^m Solution method: Picard iteration (successive substitutions). """ from dolfin import * import sys import numpy as np # Create mesh and define function space degree = int(sys.argv) divisions = [int(arg) for arg in sys.argv[2:]] d = len(divisions) domain_type = [UnitIntervalMesh, UnitSquareMesh, UnitCubeMesh] mesh = domain_type[d-1](*divisions) V = FunctionSpace(mesh, 'Lagrange', degree) # Define boundary conditions tol = 1E-14 def left_boundary(x, on_boundary): return on_boundary and abs(x) < tol def right_boundary(x, on_boundary): return on_boundary and abs(x-1) < tol Gamma_0 = DirichletBC(V, Constant(0.0), left_boundary) Gamma_1 = DirichletBC(V, Constant(1.0), right_boundary) bcs = [Gamma_0, Gamma_1] # Choice of nonlinear coefficient m = 2 def q(u): return (1+u)**m # Define variational problem for Picard iteration u = TrialFunction(V) v = TestFunction(V) u_k = interpolate(Constant(0.0), V) # previous (known) u u_k.vector().array()[:] = np.loadtxt('u_sol.txt') print(u_k.vector().array()) a = inner(q(u_k)*grad(u), grad(v))*dx f = Constant(0.0) L = f*v*dx # Picard iterations u = Function(V) # new unknown function eps = 1.0 # error measure ||u-u_k|| tol = 1.0E-5 # tolerance iter = 0 # iteration counter maxiter = 25 # max no of iterations allowed while eps > tol and iter < maxiter: iter += 1 solve(a == L, u, bcs) diff = u.vector().array() - u_k.vector().array() eps = np.linalg.norm(diff, ord=np.Inf) print 'iter=%d: norm=%g' % (iter, eps) u_k.assign(u) # update for next iteration convergence = 'convergence after %d Picard iterations' % iter if iter >= maxiter: convergence = 'no ' + convergence print """ Solution of the nonlinear Poisson problem div(q(u)*grad(u)) = f, with f=0, q(u) = (1+u)^m, u=0 at x=0 and u=1 at x=1. %s %s """ % (mesh, convergence) # Find max error u_exact = Expression('pow((pow(2, m+1)-1)*x + 1, 1.0/(m+1)) - 1', m=m, degree=1) u_e = interpolate(u_exact, V) diff = np.abs(u_e.vector().array() - u.vector().array()).max() print 'Max error:', diff
I would appreciate any help and suggestions on this. Thank you.
see Fenics Q&A: https://fenicsproject.org/qa/53/loading-an-initial-condition-from-a-vtu-file.
Basically, I removed the two added lines in the above code, save the converged solution of u in XML format and replace the line
u_k = interpolate(Constant(0.0), V) # previous (known) u
u_k = Function(V, 'saved_u.xml') # from converged solution of u
Now it works as I want.