Re: [firedrake] DG for advection-diffusion
On 5 October 2016 at 11:25, Justin Chang <jychang48@gmail.com> wrote:
So I should have noted this earlier, but the code fails only with quadilateral meshes. The code works fine on triangular meshes. So if I did this:
h = CellSize(mesh) h_avg = 0.5*(h('+')+h('-'))
on a quadrilateral mesh, i get this error:
UFL:WARNING Only know how to compute the circumradius of an affine cell. UFL:WARNING Only know how to compute the circumradius of an affine cell. Traceback (most recent call last): File "DG_test.py", line 174, in <module> solve(a_D==L_D,u0_A,bcs=bcs,solver_parameters=OS_ parameters,options_prefix="D_") File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", line 119, in solve _solve_varproblem(*args, **kwargs) File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", line 146, in _solve_varproblem appctx=appctx) File "/home/justin/Software/firedrake/src/firedrake/ firedrake/variational_solver.py", line 251, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/home/justin/Software/firedrake/src/firedrake/ firedrake/variational_solver.py", line 126, in __init__ appctx=appctx) File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", line 178, in __init__ appctx=appctx) 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 43, in compile_form kernels.append(compile_integral(integral_data, fd, prefix, parameters)) File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line 217, in compile_integral facetarea=facetarea) File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 405, in compile_ufl return map_expr_dags(translator, expressions) 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/tsfc/tsfc/ufl_utils.py", line 107, in _modified_terminal return self.modified_terminal(o) File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 254, in modified_terminal return translate(mt.terminal, mt, self.parameters) File "/home/justin/Software/firedrake/local/lib/python2.7/ site-packages/singledispatch.py", line 210, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 312, in translate_geometricquantity return geometric.translate(terminal, mt, params) File "/home/justin/Software/firedrake/local/lib/python2.7/ site-packages/singledispatch.py", line 210, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/justin/Software/firedrake/src/tsfc/tsfc/geometric.py", line 29, in translate raise AssertionError("Cannot handle geometric quantity type: %s" % type(terminal)) AssertionError: Cannot handle geometric quantity type: <class 'ufl.geometry.Circumradius'>
So I resorted to:
h_avg = Constant(1/yseed)
which also works for triangular meshes but not work quadrilateral meshes.
What do you mean by 'not work'?
Is there a more appropriate way to obtain the cell size for quads?
Justin
On Wed, Oct 5, 2016 at 4:41 AM, David Ham <David.Ham@imperial.ac.uk> wrote:
Hi Justin,
A couple of things to try:
1. Work with the weak BCs. The strong BC formulation is rather dubious from a maths perspective (not that we expect it to be wrong, but still...)
2. Cut down to single material in the first instance.
3. Drop back to 1D so you can look at actual matrices, and convince yourself that the matrix you are producing has the right invariants. EG positive definite, known couplings correct etc.
On Wed, 5 Oct 2016 at 10:14 Justin Chang <jychang48@gmail.com> wrote:
Sorry, here's a simpler version. Still has the same errors regardless of whether bcs are strongly or weakly imposed. Direct solver (-D_ksp_type preonly -D_pc_type lu) gives a complete nonsensical answer. Run problem as:
python DG_test.py 100 50
Any help appreciated, thanks!
Justin
On Wed, Oct 5, 2016 at 3:28 AM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
Hi Justin,
I don't think it's a good use of time for us to debug your full program (which runs an external python file at line 74, for example!)
I suggest you strip down the program to its absolute core. If even this doesn't run then get back in touch with a minimal failing example.
Thanks, Andrew
On 5 October 2016 at 09:22, Justin Chang <jychang48@gmail.com> wrote:
Thanks, so the problem was that the diffusivity tensor was not inside the avg() terms. That problem is fixed.
However, when the code was modified as such, i get this error:
Traceback (most recent call last): File "2D_OS_boundary_ex1.py", line 313, in <module> u0 = ADsolve(solver_D,tao_D,opt_D,u0) File "2D_OS_boundary_ex1.py", line 292, in ADsolve solver_O.solve(u0_O,b) File "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", line 141, in solve raise ConvergenceError("LinearSolver failed to converge after %d iterations with reason: %s", self.ksp.getIterationNumber(), solving_utils.KSPReasons[r]) firedrake.exceptions.ConvergenceError: ('LinearSolver failed to converge after %d iterations with reason: %s', 0, 'DIVERGED_NANORINF')
Is something wrong with my weak formulation? I even tried switching back to weakly applied Boundaries but I still get the same error.
Thanks, Justin
On Wed, Oct 5, 2016 at 3:05 AM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
Because all TrialFunctions, TestFunctions and Functions in interior facet integrals need to include restrictions (or use operators such as jump and avg)
On 5 October 2016 at 08:32, Justin Chang <jychang48@gmail.com> wrote:
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/~ksh ahbaz/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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Correction, h_avg = Constant(1/yseed) fails regardless if quad or triangle. Using h_avg = 0.5*(h('+')+h('-')) seems to be the only way to have a working code, but how do I get this equivalent for quadilateral meshes? Attached are the DG problems, one with triangular elements and one with quadrilateral elements On Wed, Oct 5, 2016 at 5:27 AM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
On 5 October 2016 at 11:25, Justin Chang <jychang48@gmail.com> wrote:
So I should have noted this earlier, but the code fails only with quadilateral meshes. The code works fine on triangular meshes. So if I did this:
h = CellSize(mesh) h_avg = 0.5*(h('+')+h('-'))
on a quadrilateral mesh, i get this error:
UFL:WARNING Only know how to compute the circumradius of an affine cell. UFL:WARNING Only know how to compute the circumradius of an affine cell. Traceback (most recent call last): File "DG_test.py", line 174, in <module> solve(a_D==L_D,u0_A,bcs=bcs,solver_parameters=OS_parameters, options_prefix="D_") File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", line 119, in solve _solve_varproblem(*args, **kwargs) File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving.py", line 146, in _solve_varproblem appctx=appctx) File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", line 251, in __init__ super(LinearVariationalSolver, self).__init__(*args, **kwargs) File "/home/justin/Software/firedrake/src/firedrake/firedrake/variational_solver.py", line 126, in __init__ appctx=appctx) File "/home/justin/Software/firedrake/src/firedrake/firedrake/solving_utils.py", line 178, in __init__ appctx=appctx) 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 43, in compile_form kernels.append(compile_integral(integral_data, fd, prefix, parameters)) File "/home/justin/Software/firedrake/src/tsfc/tsfc/driver.py", line 217, in compile_integral facetarea=facetarea) File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 405, in compile_ufl return map_expr_dags(translator, expressions) 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/tsfc/tsfc/ufl_utils.py", line 107, in _modified_terminal return self.modified_terminal(o) File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 254, in modified_terminal return translate(mt.terminal, mt, self.parameters) File "/home/justin/Software/firedrake/local/lib/python2.7/site-packages/singledispatch.py", line 210, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/justin/Software/firedrake/src/tsfc/tsfc/fem.py", line 312, in translate_geometricquantity return geometric.translate(terminal, mt, params) File "/home/justin/Software/firedrake/local/lib/python2.7/site-packages/singledispatch.py", line 210, in wrapper return dispatch(args[0].__class__)(*args, **kw) File "/home/justin/Software/firedrake/src/tsfc/tsfc/geometric.py", line 29, in translate raise AssertionError("Cannot handle geometric quantity type: %s" % type(terminal)) AssertionError: Cannot handle geometric quantity type: <class 'ufl.geometry.Circumradius'>
So I resorted to:
h_avg = Constant(1/yseed)
which also works for triangular meshes but not work quadrilateral meshes.
What do you mean by 'not work'?
Is there a more appropriate way to obtain the cell size for quads?
Justin
On Wed, Oct 5, 2016 at 4:41 AM, David Ham <David.Ham@imperial.ac.uk> wrote:
Hi Justin,
A couple of things to try:
1. Work with the weak BCs. The strong BC formulation is rather dubious from a maths perspective (not that we expect it to be wrong, but still...)
2. Cut down to single material in the first instance.
3. Drop back to 1D so you can look at actual matrices, and convince yourself that the matrix you are producing has the right invariants. EG positive definite, known couplings correct etc.
On Wed, 5 Oct 2016 at 10:14 Justin Chang <jychang48@gmail.com> wrote:
Sorry, here's a simpler version. Still has the same errors regardless of whether bcs are strongly or weakly imposed. Direct solver (-D_ksp_type preonly -D_pc_type lu) gives a complete nonsensical answer. Run problem as:
python DG_test.py 100 50
Any help appreciated, thanks!
Justin
On Wed, Oct 5, 2016 at 3:28 AM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
Hi Justin,
I don't think it's a good use of time for us to debug your full program (which runs an external python file at line 74, for example!)
I suggest you strip down the program to its absolute core. If even this doesn't run then get back in touch with a minimal failing example.
Thanks, Andrew
On 5 October 2016 at 09:22, Justin Chang <jychang48@gmail.com> wrote:
Thanks, so the problem was that the diffusivity tensor was not inside the avg() terms. That problem is fixed.
However, when the code was modified as such, i get this error:
Traceback (most recent call last): File "2D_OS_boundary_ex1.py", line 313, in <module> u0 = ADsolve(solver_D,tao_D,opt_D,u0) File "2D_OS_boundary_ex1.py", line 292, in ADsolve solver_O.solve(u0_O,b) File "/home/justin/Software/firedrake/src/firedrake/firedrake/linear_solver.py", line 141, in solve raise ConvergenceError("LinearSolver failed to converge after %d iterations with reason: %s", self.ksp.getIterationNumber(), solving_utils.KSPReasons[r]) firedrake.exceptions.ConvergenceError: ('LinearSolver failed to converge after %d iterations with reason: %s', 0, 'DIVERGED_NANORINF')
Is something wrong with my weak formulation? I even tried switching back to weakly applied Boundaries but I still get the same error.
Thanks, Justin
On Wed, Oct 5, 2016 at 3:05 AM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
Because all TrialFunctions, TestFunctions and Functions in interior facet integrals need to include restrictions (or use operators such as jump and avg)
On 5 October 2016 at 08:32, Justin Chang <jychang48@gmail.com> wrote:
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/~ksh ahbaz/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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Andrew McRae
-
Justin Chang