1.) Is this the way to read the ChangeLog?
For example, in the following website
https://fenics.readthedocs.io/projects/dolfin/en/latest/ChangeLog.html#id3
Deprecate subclassing of Expression in Python; new Python class UserExpression introduced for user overloads
Hence, if I use the fenics version 2019.2.0.dev0
I replace Expression with UserExpression
Am I correct to how to read the changelog?
2.a) Let’s say I a matrix A; I write in matlab syntax.
A = [2, 3, 4; 1, 3, 4, 0, 2, 4]
What is the difference between dolfin matrix A and the numpy converted matrix A?
In page 34. of fenics book, it is mentioned
We may convert the matrix and vector data to numpy arrays by calling the array() method as shown before.
2.b) How to print the A matrix and b vector in AU = b. What I read so far, the .vector() is deprecated, but somehow it shows in the matrix in a terminal. The vector b is not showing at all. How do I look for the proper syntax in a change log?
If I change b to b.get_local() and then print it, I can see the vector. I am confused by the difference among b, b.vector().get_local(), and b.get_local()
Here is a full code
“”"
FEniCS tutorial demo program:
Poisson equation with Dirichlet and Neumann conditions.
As dn2_p2D.py, but the linear system is explicitly formed and solved.-Laplace(u) = f on the unit square.
u = 1 + 2y^2 on x=0.
u = 2 + 2y^2 on x=1.
-du/dn = g on y=0 and y=1.
u = 1 + x^2 + 2y^2, f = -6, g = -4y.
“”"from dolfin import *
import numpyCreate mesh and define function space
mesh = UnitSquareMesh(2, 1)
V = FunctionSpace(mesh, ‘Lagrange’, 1)Define Dirichlet conditions for x=0 boundary
u_L = Expression(‘1 + 2*x[1]*x[1]’, degree = 2)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0]) < tolGamma_0 = DirichletBC(V, u_L, LeftBoundary())
Define Dirichlet conditions for x=1 boundary
u_R = Expression(‘2 + 2*x[1]*x[1]’, degree = 2)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0] - 1) < tolGamma_1 = DirichletBC(V, u_R, RightBoundary())
bcs = [Gamma_0, Gamma_1]
Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
g = Expression('-4x[1]', degree = 2)
a = inner(grad(u), grad(v))dx
L = fvdx - gvdsAssemble and solve linear system
A = assemble(a)
b = assemble(L)if mesh.num_cells() < 16:
print(‘A = assemble(a); b = assemble(L)’)
print((‘A before incorporation of essential BC:\n’, A.array()))
print((‘b before incorporation of essential BC:\n’, b.array()))for bc in bcs:
bc.apply(A, b)if mesh.num_cells() < 16:
print((‘A after incorporation of essential BC:\n’, A.array()))
print((‘b after incorporation of essential BC:\n’, b.array()))Alternative creation of the linear system
(symmetric modification of boundary conditions)
A, b = assemble_system(a, L, bcs)
if mesh.num_cells() < 16:
print(‘\nA, b = assemble_system(a, L, bcs)’)
print((‘A after incorporation of essential BC:\n’, A.array()))
print((‘b after incorporation of essential BC:\n’, b.array()))Compute solution
u = Function(V)
U = u.vector()
solve(A, U, b)#plot(u)
print((“”"
Solution of the Poisson problem -Laplace(u) = f,
with u = u0 on x=0,1 and -du/dn = g at y=0,1.
%s
“”" % mesh))Dump solution to the screen
u_nodal_values = u.vector()
u_array = u_nodal_values.array()
coor = mesh.coordinates()
for i in range(len(u_array)):
print((‘u(%8g,%8g) = %g’ % (coor[i][0], coor[i][1], u_array[i])))Exact solution:
u_exact = Expression(‘1 + x[0]x[0] + 2x[1]*x[1]’, degree = 2)
Verification
u_e = interpolate(u_exact, V)
u_e_array = u_e.vector().array()
print((‘Max error:’, numpy.abs(u_e_array - u_array).max()))Compare numerical and exact solution
center = (0.5, 0.5)
print((‘numerical u at the center point:’, u(center)))
print((‘exact u at the center point:’, u_exact(center)))#interactive()
2c.) In principle, I should be able to print and see the mesh with the followign command in the following section of the above code
print(“”"
Solution of the Poisson problem -Laplace(u) = f,
with u = u0 on x=0,1 and -du/dn = g at y=0,1.
%s
“”" % mesh)
but I get
<dolfin.cpp.generation.UnitSquareMesh object at 0x7f230c0dde50>
I am not sure what is going on?
A minimal working code is
from dolfin import *
mesh = IntervalMesh(5, -2, 0)
plot(mesh)
plt.show()
print('%.2f' % mesh)
3.) How to plot the unknown in 1D problem? Like the familiary y = sin(x), problem, how do I plot here? I input the dimension of the problem via command line.
For 3D, I type (works great):
python3 d1_p2D_GKC.py 2 10 3 4
For 2D I type (works great):
python3 d1_p2D_GKC.py 2 10 3
For 1D I type (Does not plot):
python3 d1_p2D_GKC.py 2 10
Here is the minimal working code.
from __future__ import print_function
from dolfin import *
import matplotlib.pyplot as plt
import sys
degree = int(sys.argv[1])
print(degree)
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)
# Create mesh and define function space
#mesh = UnitCubeMesh(6, 4, 5)
#V = FunctionSpace(mesh, 'Lagrange', 1)
# Define boundary conditions
u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree = 2)
tol = 1E-14 # tolerance for coordiante comparisons
def Dirichlet_boundary0(x, on_boundary):
return on_boundary and abs(x[0]) < tol
def Dirichlet_boundary1(x, on_boundary):
return on_boundary and abs(x[0] - 1) < tol
bc0 = DirichletBC(V, Constant(0), Dirichlet_boundary0)
bc1 = DirichletBC(V, Constant(1), Dirichlet_boundary1)
bcs = [bc0, bc1]
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-2.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx
# Compute solution
#u = Function(V)
#solve(a == L, u, bcs)
u = Function(V)
solve(a == L, u, bcs, solver_parameters = {'linear_solver': 'cg', 'preconditioner': 'ilu'})
#File('poisson/parameters.xml') << parameters
from vedo.dolfin import plot as plot
plot(u)
#plt.show()
#plot(mesh)
#plt.show()
# Dump solution to file in VTK format
#file = File("poisson/solution.pvd")
#file << u
4.) What does the ; mean in F(u;v)?
5.) As per previous discussion, I think I have correctly set the iterative solver, however, I get this error
File “vp2_np.py”, line 78, in
prm[“linear_solver”] = ‘gmres’
RuntimeError: Parameter linear_solver not found in Parameters object
here is the minimal working code. I typed
python3 vp2_np.py a g 1 3 4
"""
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: automatic, i.e., by a NonlinearVariationalProblem/Solver
(Newton method), with automatic UFL computation of the derivative.
"""
from dolfin import *
import numpy, sys
# Usage: ./vp2_np.py m|a |g|l degree nx ny nz
# Example: ./vp2_np.py m l 1 3 4
J_comp = sys.argv[1] # m (manual) or a (automatic) computation of J
answer = sys.argv[2] # g (GMRES) or l (sparse LU) solver
iterative_solver = True if answer == 'g' else False
# Create mesh and define function space
degree = int(sys.argv[3])
divisions = [int(arg) for arg in sys.argv[4:]]
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
def Dq(u):
return m*(1+u)**(m-1)
# Define variational problem
v = TestFunction(V)
du = TrialFunction(V)
u_ = Function(V) # most recently computed solution
F = inner(q(u_)*grad(u_), grad(v))*dx
# J must be a Jacobian (Gateaux derivative in direction of du)
if J_comp == 'm':
J = inner(q(u_)*grad(du), grad(v))*dx + \
inner(Dq(u_)*du*grad(u_), grad(v))*dx
else:
J = derivative(F, u_, du)
# Compute solution
problem = NonlinearVariationalProblem(F, u_, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
#info(prm, True)
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['maximum_iterations'] = 25
prm['newton_solver']['relaxation_parameter'] = 1.0
if iterative_solver:
prm['linear_solver'] = 'gmres'
prm['preconditioner'] = 'ilu'
prm['krylov_solver']['absolute_tolerance'] = 1E-9
prm['krylov_solver']['relative_tolerance'] = 1E-7
prm['krylov_solver']['maximum_iterations'] = 1000
prm['krylov_solver']['monitor_convergence'] = True
prm['krylov_solver']['nonzero_initial_guess'] = False
prm['krylov_solver']['gmres']['restart'] = 40
prm['krylov_solver']['preconditioner']['same_nonzero_pattern'] = True
prm['krylov_solver']['preconditioner']['ilu']['fill_level'] = 0
PROGRESS = 16
set_log_level(LogLevel.PROGRESS)
solver.solve()
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
""" % mesh)
# Find max error
u_exact = Expression('pow((pow(2, m+1)-1)*x[0] + 1, 1.0/(m+1)) - 1', m=m, degree = 2)
u_e = interpolate(u_exact, V)
import numpy
diff = numpy.abs(u_e.vector().get_local() - u_.vector().get_local()).max()
print ('Max error:', diff)
from vedo.dolfin import plot
plot(u_)
# Define variational problem
#v = TestFunction(V)
#du = TrialFunction(V)
#u = Function(V) # the unknown
#F = inner(q(u)*grad(u), grad(v))*dx
#J = derivative(F, u, du)
6.) I came across this, while doing the problem from the fencis book. I have a minimum working example at the end of the problem. I don’t know how to install scitools following instructions here
https://scitools.org.uk/iris/docs/v1.9.2/installing.html
In brief, I
apt-get install python3-numpy
python3 setup.py install
ModuleNotFoundError: No module named ‘pyke’
error: command ‘/usr/bin/python3’ failed with exit status 1
sudo apt-get install pyke
Reading package lists… Done
Building dependency tree
Reading state information… Done
E: Unable to locate package pyke
Minimum working example is
"""Temperature variations in the ground."""
from dolfin import *
import sys, numpy, time
# Usage: sin_daD.py degree D nx ny nz
# Example: sin_daD.py 1 1.5 4 40
# Create mesh and define function space
degree = int(sys.argv[1])
D = float(sys.argv[2])
W = D/2.0
divisions = [int(arg) for arg in sys.argv[3:]]
print('degree=%d, D=%g, W=%g, %s cells' % \
(degree, D, W, 'x'.join(sys.argv[3:])))
d = len(divisions) # no of space dimensions
if d == 1:
mesh = IntervalMesh(divisions[0], -D, 0)
elif d == 2:
mesh = RectangleMesh(Point(-W/2, -D), Point(W/2, 0), divisions[0], divisions[1])
elif d == 3:
mesh = BoxMesh(Point(-W/2, -W/2), Point(-D, W/2), Point(W/2, 0),
divisions[0], divisions[1], divisions[2])
V = FunctionSpace(mesh, 'Lagrange', degree)
# Define boundary conditions
T_R = 0
T_A = 1.0
omega = 2*pi
# soil:
T_R = 10
T_A = 10
omega = 7.27E-5
T_0 = Expression('T_R + T_A*sin(omega*t)',
T_R=T_R, T_A=T_A, omega=omega, t=0.0, degree = 2)
def surface(x, on_boundary):
return on_boundary and abs(x[d-1]) < 1E-14
bc = DirichletBC(V, T_0, surface)
period = 2*pi/omega
t_stop = 5*period
#t_stop = period/4
dt = period/14
theta = 1
kappa_0 = 2.3 # KN/s (K=Kelvin, temp. unit), for soil
kappa_0 = 12.3 # KN/s (K=Kelvin, temp. unit)
kappa_1 = 1E+4
kappa_1 = kappa_0
kappa_str = {}
kappa_str[1] = 'x[0] > -D/2 && x[0] < -D/2 + D/4 ? kappa_1 : kappa_0'
kappa_str[2] = 'x[0] > -W/4 && x[0] < W/4 '\
'&& x[1] > -D/2 && x[1] < -D/2 + D/4 ? '\
'kappa_1 : kappa_0'
kappa_str[3] = 'x[0] > -W/4 && x[0] < W/4 '\
'x[1] > -W/4 && x[1] < W/4 '\
'&& x[2] > -D/2 && x[2] < -D/2 + D/4 ?'\
'kappa_1 : kappa_0'
# Alternative way of defining the kappa function
class Kappa(Expression):
def eval(self, value, x):
"""x: spatial point, value[0]: function value."""
d = len(x) # no of space dimensions
material = 0 # 0: outside, 1: inside
if d == 1:
if -D/2. < x[d-1] < -D/2. + D/4.:
material = 1
elif d == 2:
if -D/2. < x[d-1] < -D/2. + D/4. and \
-W/4. < x[0] < W/4.:
material = 1
elif d == 3:
if -D/2. < x[d-1] < -D/2. + D/4. and \
-W/4. < x[0] < W/4. and -W/4. < x[1] < W/4.:
material = 1
value[0] = kappa_0 if material == 0 else kappa_1
# Physical parameters
kappa = Expression(kappa_str[d],
D=D, W=W, kappa_0=kappa_0, kappa_1=kappa_1, degree = 2)
#kappa = Kappa(V)
# soil:
rho = 1500
c = 1600
print('Thermal diffusivity:', kappa_0/rho/c)
# Define initial condition
T_1 = interpolate(Constant(T_R), V)
# Define variational problem
T = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = rho*c*T*v*dx + theta*dt*kappa*\
inner(nabla_grad(v), nabla_grad(T))*dx
L = (rho*c*T_1*v + dt*f*v -
(1-theta)*dt*kappa*inner(nabla_grad(v), nabla_grad(T)))*dx
A = assemble(a)
b = None # variable used for memory savings in assemble calls
# Compute solution
T = Function(V) # unknown at the new time level
# hack for plotting (first surface must span all values (?))
dummy = Expression('T_R - T_A/2.0 + 2*T_A*(x[%g]+D)' % (d-1),
T_R=T_R, T_A=T_A, D=D, degree = 2)
# Make all plot commands inctive
import scitools.misc
plot = scitools.misc.DoNothing(silent=True)
# Need initial dummy plot
viz = plot(dummy, axes=True,
title='Temperature', wireframe=False)
viz.elevate(-65)
#time.sleep(1)
viz.update(T_1)
import scitools.BoxField
start_pt = [0]*d; start_pt[-1] = -D # start pt for line plot
import scitools.easyviz as ev
def line_plot():
"""Make a line plot of T along the vertical direction."""
if T.ufl_element().degree() != 1:
T2 = interpolate(T, FunctionSpace(mesh, 'Lagrange', 1))
else:
T2 = T
T_box = scitools.BoxField.dolfin_function2BoxField(
T2, mesh, divisions, uniform_mesh=True)
#T_box = scitools.BoxField.update_from_dolfin_array(
# T.vector().array(), T_box)
coor, Tval, fixed, snapped = \
T_box.gridline(start_pt, direction=d-1)
# Use just one ev.plot command, not hold('on') and two ev.plot
# etc for smooth movie on the screen
if kappa_0 == kappa_1: # analytical solution?
ev.plot(coor, Tval, 'r-',
coor, T_exact(coor), 'b-',
axis=[-D, 0, T_R-T_A, T_R+T_A],
xlabel='depth', ylabel='temperature',
legend=['numerical', 'exact, const kappa=%g' % kappa_0],
legend_loc='upper left',
title='t=%.4f' % t)
else:
ev.plot(coor, Tval, 'r-',
axis=[-D, 0, T_R-T_A, T_R+T_A],
xlabel='depth', ylabel='temperature',
title='t=%.4f' % t)
ev.savefig('tmp_%04d.png' % counter)
time.sleep(0.1)
def T_exact(x):
a = sqrt(omega*rho*c/(2*kappa_0))
return T_R + T_A*exp(a*x)*numpy.sin(omega*t + a*x)
n = FacetNormal(mesh) # for flux computation at the top boundary
flux = -kappa*dot(nabla_grad(T), n)*ds
t = dt
counter = 0
while t <= t_stop:
print('time =', t)
b = assemble(L, tensor=b)
T_0.t = t
help = interpolate(T_0, V)
print('T_0:', help.vector().array()[0])
bc.apply(A, b)
solve(A, T.vector(), b)
viz.update(T)
#viz.write_ps('diff2_%04d' % counter)
line_plot()
total_flux = assemble(flux)
print('Total flux:', total_flux)
t += dt
T_1.assign(T)
counter += 1
viz = plot(T, title='Final solution')
viz.elevate(-65)
viz.azimuth(-40)
viz.update(T)
#interactive()