← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/demo/pde/convection-diffusion/python/demo.py

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/demo/pde/convection-diffusion/python/demo.py
1
 
# This demo solves the time-dependent convection-diffusion equation by
2
 
# a least-squares stabilized cG(1)cG(1) method. The velocity field used
3
 
# in the simulation is the output from the Stokes (Taylor-Hood) demo.
4
 
# The sub domains for the different boundary conditions are computed
5
 
# by the demo program in src/demo/subdomains.
6
 
#
7
 
# Modified by Anders Logg, 2008
8
 
 
9
 
__author__ = "Kristian B. Oelgaard (k.b.oelgaard@tudelft.nl)"
10
 
__date__ = "2007-11-14 -- 2008-01-04"
11
 
__copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
12
 
__license__  = "GNU LGPL Version 2.1"
13
 
 
14
 
from dolfin import *
15
 
 
16
 
# Load mesh and create finite elements
17
 
mesh = Mesh("../../../../../data/meshes/dolfin-2.xml.gz")
18
 
scalar = FiniteElement("Lagrange", "triangle", 1)
19
 
vector = VectorElement("Lagrange", "triangle", 2)
20
 
 
21
 
# Load subdomains and velocity
22
 
sub_domains = MeshFunction("uint", mesh, "../subdomains.xml.gz");
23
 
velocity = Function(vector, "../velocity.xml.gz");
24
 
 
25
 
# Initialise source function and previous solution function
26
 
f = Function(scalar, mesh, 0.0)
27
 
u0 = Function(scalar, mesh, 0.0)
28
 
 
29
 
# Parameters
30
 
T = 1.0
31
 
k = 0.05
32
 
t = k
33
 
c = 0.005
34
 
 
35
 
# Test and trial functions
36
 
v = TestFunction(scalar)
37
 
u = TrialFunction(scalar)
38
 
 
39
 
# Functions
40
 
u0 = Function(scalar, mesh, Vector())
41
 
u1 = Function(scalar, mesh, Vector())
42
 
 
43
 
# Variational problem
44
 
a = v*u*dx + 0.5*k*(v*dot(velocity, grad(u))*dx + c*dot(grad(v), grad(u))*dx)
45
 
L = v*u0*dx - 0.5*k*(v*dot(velocity, grad(u0))*dx + c*dot(grad(v), grad(u0))*dx) + k*v*f*dx
46
 
 
47
 
# Set up boundary condition
48
 
g = Function(mesh, 1.0)
49
 
bc = DirichletBC(g, sub_domains, 1)
50
 
 
51
 
# Assemble matrix
52
 
A = assemble(a, mesh)
53
 
 
54
 
# Output file
55
 
out_file = File("temperature.pvd")
56
 
 
57
 
# Time-stepping
58
 
while ( t < T ):
59
 
 
60
 
    # Assemble vector and apply boundary conditions
61
 
    b = assemble(L, mesh)
62
 
    bc.apply(A, b, a)
63
 
    
64
 
    # Solve the linear system
65
 
    solve(A, u1.vector(), b)
66
 
    
67
 
    # Save the solution to file
68
 
    out_file << u1
69
 
 
70
 
    # Move to next interval
71
 
    t += k
72
 
    u0.assign(u1)
73
 
 
74
 
# Plot solution
75
 
plot(u1)