Re: [firedrake] DG for advection-diffusion
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
Okay thanks, I am attempting to do the diffusion part only at the moment under DG, but I am getting these errors: $ python 2D_OS_boundary_ex1.py 100 50 0 0 0 -D_pc_type bjacobi UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, creating DQ element. UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, creating DQ element. UFL:WARNING Discontinuous Lagrange element requested on quadrilateral, creating DQ element. UFL:ERROR Discontinuous type Coefficient must be restricted. Traceback (most recent call last): File "2D_OS_boundary_ex1.py", line 255, in <module> A_D = assemble(a_D,bcs=bc_D) File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", line 81, in assemble inverse=inverse, mat_type=mat_type, appctx=appctx) File "<decorator-gen-291>", line 2, in _assemble File "/home/justin/Software/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper return f(*args, **kwargs) File "/home/justin/Software/firedrake/src/firedrake/firedrake/assemble.py", line 120, in _assemble inverse=inverse) File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 189, in compile_form number_map).kernels File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line 198, in __new__ obj = make_obj() File "/home/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line 188, in make_obj obj.__init__(*args, **kwargs) File "/home/justin/Software/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 110, in __init__ tree = tsfc_compile_form(form, prefix=name, parameters=parameters) File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line 36, in compile_form fd = ufl_utils.compute_form_data(form) File "/home/justin/Software/firedrake/src/tsfc/tsfc/ufl_utils.py", line 47, in compute_form_data do_estimate_degrees=do_estimate_degrees, File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/compute_form_data.py", line 305, in compute_form_data form = apply_restrictions(form) File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", line 211, in apply_restrictions only_integral_type=integral_types) File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 58, in map_integrand_dags form, only_integral_type) File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 39, in map_integrands for itg in form.integrals()] File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 46, in map_integrands return itg.reconstruct(function(itg.integrand())) File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/map_integrands.py", line 57, in <lambda> return map_integrands(lambda expr: map_expr_dag(function, expr, compress), File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 37, in map_expr_dag result, = map_expr_dags(function, [expression], compress=compress) File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/map_dag.py", line 84, in map_expr_dags r = handlers[v._ufl_typecode_](v) File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", line 102, in reference_value g = self(f) File "/home/justin/Software/firedrake/src/ufl/ufl/corealg/multifunction.py", line 95, in __call__ return self._handlers[o._ufl_typecode_](o, *args) File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", line 129, in coefficient return self._require_restriction(o) File "/home/justin/Software/firedrake/src/ufl/ufl/algorithms/apply_restrictions.py", line 59, in _require_restriction "Discontinuous type %s must be restricted." % o._ufl_class_.__name__) File "/home/justin/Software/firedrake/src/ufl/ufl/assertions.py", line 51, in ufl_assert error(*message) File "/home/justin/Software/firedrake/src/ufl/ufl/log.py", line 167, in error raise self._exception_type(self._format_raw(*message)) ufl.log.UFLException: Discontinuous type Coefficient must be restricted. I can't seem to figure out why. Attached is the code. Run as: python 2D_OS_boundary_ex1.py 100 50 0 0 0 Thanks, Justin On Thu, Sep 29, 2016 at 4:49 AM, David Ham <David.Ham@imperial.ac.uk> wrote:
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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
David Ham
-
Justin Chang