--- Begin Message ---
------------------------------------------------------------
revno: 4755
committer: Garth N. Wells <gnw20@xxxxxxxxx>
branch nick: dolfin-all
timestamp: Thu 2010-05-27 11:29:01 +0100
message:
Add SUPG stabilisation to adevction-diffusion Python demo.
I had to change the velocity field from the computed Stokes field to a constant vector. Also, the boundary conditions had to be changed to avoid a discontinuity. To use the Stokes field, a shock capturing terms needs to be added since the Stokes problem has some very high gradients.
C++ demo still needs to be updated.
modified:
demo/pde/advection-diffusion/python/demo.py
--
lp:dolfin
https://code.launchpad.net/~dolfin-core/dolfin/main
Your team DOLFIN Core Team is subscribed to branch lp:dolfin.
To unsubscribe from this branch go to https://code.launchpad.net/~dolfin-core/dolfin/main/+edit-subscription
=== modified file 'demo/pde/advection-diffusion/python/demo.py'
--- demo/pde/advection-diffusion/python/demo.py 2010-05-27 09:23:29 +0000
+++ demo/pde/advection-diffusion/python/demo.py 2010-05-27 10:29:01 +0000
@@ -1,9 +1,12 @@
# This demo solves the time-dependent convection-diffusion equation by
-# a least-squares stabilized cG(1)cG(1) method. The velocity field used
+# 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.
#
+# FIXME: Add shock capturing term and then revert back to the Stokes
+# velocity
+#
# Modified by Anders Logg, 2008
# Modified by Johan Hake, 2008
# Modified by Garth N. Wells, 2009
@@ -15,16 +18,24 @@
from dolfin import *
+def boundary_value(n):
+ if n < 10:
+ return float(n)/10.0
+ else:
+ return 1.0
+
# Load mesh and subdomains
mesh = Mesh("../mesh.xml.gz")
sub_domains = MeshFunction("uint", mesh, "../subdomains.xml.gz");
+h = CellSize(mesh)
# Create FunctionSpaces
Q = FunctionSpace(mesh, "CG", 1)
V = VectorFunctionSpace(mesh, "CG", 2)
# Create velocity Function
-velocity = Function(V, "../velocity.xml.gz");
+#velocity = Function(V, "../velocity.xml.gz");
+velocity = Constant( (-1.0, 0.0) )
# Initialise source function and previous solution function
f = Constant(0.0)
@@ -32,20 +43,30 @@
# Parameters
T = 5.0
-k = 0.1
-t = k
-c = 0.005
+dt = 0.1
+t = dt
+c = 0.00005
# Test and trial functions
v = TestFunction(Q)
u = TrialFunction(Q)
-# Variational problem
-a = v*u*dx + 0.5*k*(v*dot(velocity, grad(u))*dx + c*dot(grad(v), grad(u))*dx)
-L = v*u0*dx - 0.5*k*(v*dot(velocity, grad(u0))*dx + c*dot(grad(v), grad(u0))*dx) + k*v*f*dx
+# Resisual (split into LHS/RHS parts)
+r_a = u + 0.5*dt*dot(velocity, grad(u)) - 0.5*dt*div(grad(u))
+r_L = u0 - 0.5*dt*dot(velocity, grad(u0)) + 0.5*dt*div(grad(u0)) + dt*f
+
+# Galerkin variational problem
+a_g = v*u*dx + 0.5*dt*(v*dot(velocity, grad(u))*dx + c*dot(grad(v), grad(u))*dx)
+L_g = v*u0*dx - 0.5*dt*(v*dot(velocity, grad(u0))*dx + c*dot(grad(v), grad(u0))*dx) + dt*v*f*dx
+
+# Add SUPG stabilisation
+vnorm = sqrt(dot(velocity, velocity))
+supg = (h/2.0*vnorm)*dot(velocity, grad(v))
+a = a_g + supg*r_a*dx
+L = L_g + supg*r_L*dx
# Set up boundary condition
-g = Constant(1.0)
+g = Constant( boundary_value(0) )
bc = DirichletBC(Q, g, sub_domains, 1)
# Assemble matrix
@@ -76,8 +97,9 @@
# Save the solution to file
out_file << u
- # Move to next interval
- t += k
+ # Move to next interval and adjust boundary condition
+ t += dt
+ g.assign( boundary_value( int(t/dt) ) )
# Hold plot
interactive()
--- End Message ---
