Yes it is, and using Constant(1/float(yseed)) now results in a fully working code! I originally tried Constant(float(1/yseed)) so I thought that I already eliminated that possibility.

Thanks!

On Wed, Oct 5, 2016 at 5:46 AM, Andrew McRae <A.T.T.McRae@bath.ac.uk> wrote:
Is yseed an integer?  If so, 1/yseed presumably gives 0.

On 5 October 2016 at 11:41, Justin Chang <jychang48@gmail.com> wrote:
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/~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


_______________________________________________
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