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(a)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(a)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(a)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(a)imperial.ac.uk>
>>> wrote:
>>>
>>>>
>>>>
>>>> On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)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(a)imperial.ac.uk
>>>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>>>
>>>>
>>>
>>
>> _______________________________________________
>> firedrake mailing list
>> firedrake(a)imperial.ac.uk
>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>
>>
>
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(a)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(a)imperial.ac.uk>
> wrote:
>
>>
>>
>> On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)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(a)imperial.ac.uk
>> https://mailman.ic.ac.uk/mailman/listinfo/firedrake
>>
>>
>
On Thu, 29 Sep 2016 at 10:34 Justin Chang <jychang48(a)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
>
>
>