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