creating DG0 feilds computationally without floating point error.

[karimnorouzi]

when we want to create a DG0 field like:

assemble(sum(Constant(v)/CellVolume(mesh)*dx(i+1) for i, v in enumerate(val)))

the floating point error sometimes creates large distortion.
as per recommandation from Colin Cotter:

"This could be fixed by adding floor/ceil to UFL (and then handling in the form compiler).
6:03
If you raise this as an issue in Firedrake then someone can give you advise on how to do it."

How to do this ?

thanks.

[dham]

Can you not get what you want with conditional?

There's nothing wrong with adding floor and ceil to UFL other than it'll cause you a lot of legwork (the differentiation rules, for example). If you can get away with conditional then your life will be easier.

[karimnorouzi]

Is there a way to use something like this:

def Parameter(vals, FS):
    par = Function(FS)
    for j in range(len(vals)):
        par_loop(("{[i] : 0 <= i < f.dofs}", "f[i, 0] = {}".format(vals[j])),
             dx(j+1),
             {"f": (par, WRITE)})
    return par

mu = Parameter([ 1, 2] , fd.FunctionSpace(mesh, "DG", 0))

but correctly ? This obviously does not work for DG0

[pbrubeck]

Could please you send an example where the round-off error is noticeable?
You should only get large round-off errors if the values you are prescribing are astronomically large:

mesh = UnitSquareMesh(1000, 2)
V = FunctionSpace(mesh, "DG", 0)
val = 2E10
c = assemble(val*TestFunction(V)/CellVolume(mesh)*dx)
max(abs(c.dat.data[:] - val))  # --> np.float64(3.814697265625e-06)

[pbrubeck]

Here's how you assign int data on a subdomain:

from firedrake import *

def assign_subdomain(f, val, subdomain):
    mesh = f.function_space().mesh()
    f.assign(val, subset=mesh.cell_subset(subdomain))

mesh = UnitSquareMesh(1, 1)
q = Function(FunctionSpace(mesh, "DG", 0), dtype=int)
q.dat.data[1] = 1
rmesh = RelabeledMesh(mesh, (q,), (1,))
W = FunctionSpace(rmesh, "DG", 0)
w = Function(W, dtype=int)

assign_subdomain(w, 1E10, 1)

print(w.dat.data) # array([          0, 10000000000])

We should probably make this more user-friendly.

[pbrubeck]

@ksagiyam could we allow subset to be an iterable of subdomain_ids? And internally construct the subset via mesh.cell_subset?

[karimnorouzi]

dtype = int does not work with VTKFile

But the solution is very good. I like it.

[karimnorouzi]

def assign_subdomain(f, vals, subdomains):

mesh = f.function_space().mesh()
for v, s in zip(vals, subdomains):
            f.assign(v, subset=mesh.cell_subset(s))

is this ok for multiple valsues and subdomains ?

[ksagiyam]

@pbrubeck I think the potential problem is that, if we assign to a CG function with some cell subdomain IDs, what should happen at the subdomain interface is not clear. It is clear for DG functions.

[pbrubeck]

Wouldn't the obvious operation assign dofs on the closure of each cell in the subdomain?

[ksagiyam]

I think cg_assignee.assign(cg_assigner, subdomain) is well defined only when cg_assignee and cg_assigner have the same values at the subdomain interface, and we can not force users' code to satisfy this condition (unless we inspect data). Our current implementation might not be user friendly, but causes no ambiguity regarding this point.

[karimnorouzi]

this is Ok for DG functions right ?
for CG functions this should not work as the other element sharing the interface can be affected ?

[pbrubeck]

this is Ok for DG functions right ? for CG functions this should not work as the other element sharing the interface can be affected ?

Yes, for DG there's no confusion

[pbrubeck]

@ksagiyam I think for CG we can agree on a well defined operation, by considering the subspace of basis functions supported on the closure of cells in the subdomain.

DirichletBC uses this subset mechanism for continuous Functions, except that subdomains are constructed as subsets of the facets (we even allow them to be taken from the interior). That is why the ordering of bcs gives different results: [DirichletBC(V, g1, sub1), DirichletBC(V, g2, sub2)] writes g2 on the intersection of sub1 and sub2, while the reversed order would write g1 instead. So we are already exposing similar ambiguities to the user.