Is it much work to add the shock capturing term? There's a point in having the velocity field from the Stokes demo as that illustrates that one can feed the output from one solution into another.
-- Anders
--- Begin Message --------------------------------------------------------------- revno: 4755 committer: Garth N. Wells <[email protected]> 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 ---
signature.asc
Description: Digital signature
_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
