### Feed values from data file to a interpolate function

148

views

0

Hi all,

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,

The complete script is as follows:

I would appreciate any help and suggestions on this. Thank you.

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[1])
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[0]) < tol
def right_boundary(x, on_boundary):
return on_boundary and abs(x[0]-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[0] + 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.

Community: FEniCS Project

### 1 Answer

3

I didn't notice that there already has a solution for this,

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

with

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`

with

`u_k = Function(V, 'saved_u.xml') # from converged solution of u`

.

Now it works as I want.

Please login to add an answer/comment or follow this question.