On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48@gmail.com> wrote:
Hi all,
I have written here the Discontinuous Galerkin method for 1D advection-diffusion equations using backward euler time discretization:
mesh = UnitIntervalMesh(100) V = FunctionSpace(mesh, 'DG', 1) u = TrialFunction(V) v = TestFunction(V)
u0 = Function(V, name="initial") sol = Function(V, name="solution")
# Diffusivity D = Constant(0.005)
# Velocity vel = as_vector((Constant(1.0),))
# Forcing function f = interpolate(Expression(....),V)
# Penalty k = 1 d = 1 gamma = Constant((k+1)*(k+d)/d)*FacetArea(mesh)/avg(CellVolume(mesh)) # based on the formula in http://webpages.sdsmt.edu/~kshahbaz/shahbazi_05.pdf
# Other parameters n = FacetNormal(mesh) h = CellSize(mesh) vn = 0.5*(dot(vel,n) + abs(dot(vel,n)))
# Weak form a = v*u/dt*dx - dot(grad(v),vel*u)*dx + dot(jump(v),vn('+')*u('+')-vn('-')*u('-'))*dS + dot(v, vn*u)*ds \ + dot(grad(v),D*grad(u))*dx - dot(jump(v,n),avg(D*grad(u)))*dS - dot(avg(D*grad(v)),jump(u,n))*dS \ + gamma*dot(jump(v,n),jump(u,n))*dS L = v*f*dx + v*u0/dt*dx
...
I have a couple questions:
1) If say I am solving for the transport of 4 chemical species, Normally I would create a mixed function space like V_mixed = V*V*V*V, split u into u1, u2, u3, and u4, write bilinear and linear forms for each of these individual components and add them up. That is, a = a1 + a2 + a3 + a4, L = L1 + L2 + L3 + L4. This can become quite a lot of code if I end up having to solve ~20 chemical species eventually. The catch is that the velocity and diffusivity will be the same for all species but the forcing function and boundary conditions are going to be different. Is there a more efficient of writing the bilinear and linear form instead of having to copy and paste the a = ... and L = ... multiple times and changing the trial/test functions?
a) If you are independently advecting lots of species (ie no reaction) then you probably want to just solve the equations separately, rather than making a massive mixed function space. b) You can avoid lots of retyping by just using normal Python functions to construct forms. You could, for example, write a function which takes a test and trial function, and maybe some BC information, and returns the appropriate forms.
2) I see that in FEniCS's undocumented dg-advection-diffusion example, they enforce the DirichletBC's strongly. Is it possible to enforce them weakly like the SIPG or DG upwinding techniques for the pure diffusion or advection equations?
Yes, just substitute boundary conditions for surface integrals in the usual way.
Thanks, Justin