1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
|
# Copyright (C) 2007 Kristian B. Oelgaard
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
#
# Modified by Anders Logg, 2008
# Modified by Johan Hake, 2008
# Modified by Garth N. Wells, 2009
#
# This demo solves the time-dependent convection-diffusion equation by
# a SUPG stabilized method. The velocity field used in the simulation
# is the output from the Stokes (Taylor-Hood) demo. The sub domains
# for the different boundary conditions are computed by the demo
# program in src/demo/subdomains.
from dolfin import *
import matplotlib.pyplot as plt
def boundary_value(n):
if n < 10:
return float(n)/10.0
else:
return 1.0
# Load mesh and subdomains
mesh = Mesh("../dolfin_fine.xml.gz")
sub_domains = MeshFunction("size_t", mesh, "../dolfin_fine_subdomains.xml.gz");
h = CellDiameter(mesh)
# Create FunctionSpaces
Q = FunctionSpace(mesh, "CG", 1)
V = VectorFunctionSpace(mesh, "CG", 2)
# Create velocity Function from file
velocity = Function(V);
File("../dolfin_fine_velocity.xml.gz") >> velocity
# Initialise source function and previous solution function
f = Constant(0.0)
u0 = Function(Q)
# Parameters
T = 5.0
dt = 0.1
t = dt
c = 0.00005
# Test and trial functions
u, v = TrialFunction(Q), TestFunction(Q)
# Mid-point solution
u_mid = 0.5*(u0 + u)
# Residual
r = u - u0 + dt*(dot(velocity, grad(u_mid)) - c*div(grad(u_mid)) - f)
# Galerkin variational problem
F = v*(u - u0)*dx + dt*(v*dot(velocity, grad(u_mid))*dx \
+ c*dot(grad(v), grad(u_mid))*dx)
# Add SUPG stabilisation terms
vnorm = sqrt(dot(velocity, velocity))
F += (h/(2.0*vnorm))*dot(velocity, grad(v))*r*dx
# Create bilinear and linear forms
a = lhs(F)
L = rhs(F)
# Set up boundary condition
g = Constant(boundary_value(0))
bc = DirichletBC(Q, g, sub_domains, 1)
# Assemble matrix
A = assemble(a)
bc.apply(A)
# Create linear solver and factorize matrix
solver = LUSolver(A)
# Output file
out_file = File("results/temperature.pvd")
# Set intial condition
u = u0
# Time-stepping, plot initial condition.
i = 0
plt.figure()
plot(u, title=r"t = {0:1.1f}".format(0.0))
i += 1
while t - T < DOLFIN_EPS:
# Assemble vector and apply boundary conditions
b = assemble(L)
bc.apply(b)
# Solve the linear system (re-use the already factorized matrix A)
solver.solve(u.vector(), b)
# Copy solution from previous interval
u0 = u
# Plot solution
if i % 5 == 0:
plt.figure()
plot(u, title=r"t = {0:1.1f}".format(t))
# Save the solution to file
out_file << (u, t)
# Move to next interval and adjust boundary condition
t += dt
i += 1
g.assign(boundary_value(int(t/dt)))
plt.show()
|
