On Fri, 13 Nov 2015 at 16:29 Justin Chang <jychang48@gmail.com> wrote:

Yeah I found out about errornorm recently. Have been using that now.

But as for the original issue... would D.assign(assemble(div(q) * e * dx) suffice?

I'm afraid not. assemble(div(q)*e*dx) is \int_\Omega div(q) e dx \forall e \in E where E is the basis for DG0. This is not the same thing as the projection of div(q) into DG0. In fact, for the very specific case of DG0, this will differ from the correct answer in each element by a factor of the volume of that element. This is why dividing by \int_\Omega e dx \forall e \in E is necessary to achieve the correct answer. To go into slightly more detail, note that the e_i \in E is: / 1 if x \in T_i (the i-th element) e_i(x) = | \ 0 otherwise. Hence: \int_\Omega \div(q) e_i dx == \int_{T_i} div(q)) dx and \int_\Omega e_i dx == \int_{T_i} dx == vol(T_i) but it must be the case that: \int_T_i D_i e_i dx == \int_{T_i} div(q)) dx and note: \int_T_i D_i e_i dx == D_i \int_T_i dx (since D_i is just a scalar, not a function of x) hence by rearranging: D_i = \int_\Omega \div(q) e_i dx/ \int_\Omega e_i dx and hence: D = \int_\Omega \div(q) e_i dx/ \int_\Omega e_i dx \forall e_i \in E

On Fri, Nov 13, 2015 at 8:27 AM, David Ham <David.Ham@imperial.ac.uk> wrote:

...

On Thu, 12 Nov 2015 at 21:48 Justin Chang <jychang48@gmail.com> wrote:

...

Hi guys,

So I was attempting to do something similar:

L2_error_norm = norm(assemble(v*(u-u_exact)*dx)/assemble(v*dx))

where v is test function, u is FE solution, and u_exact is the analytical solution. All of which are CG1 space. I get an error saying "fl.log.UFLException: Division by non-scalar is undefined"

Know what's up?

Yes, that's either a missing feature or a bug depending on how you look at it. However I think you are just confused in what you are trying to do.

What you probably meant was:

errornorm(u, u_exact)

(errornorm is a Firedrake-provided function).

What errornorm is basically doing is:

assemble((u-u_exact)**2*dx))**0.5

Note that this is a zero-form (no test or trial function; the answer is a scalar) so there are no solves, or divisions by lumped mass.

errornorm actually does somewhat more sophisticated things such as increasing the polynomial order of the functions in order to avoid certain sorts of error, so it's typically better to use errornorm than to do it yourself.

Regards,

David

...

Thanks, Justin

On Thu, Nov 12, 2015 at 4:33 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:

...

On 12/11/15 11:27, Justin Chang wrote:

...

David,

So if D.assign(assemble(div(q) * e * dx)/assemble(e * dx)) returns cell-wise div(q), what's the denominator "/assemble(e * dx)" for?

It's just the normal FE L2 projection:

you want:

u = div(q)

So you hit both sides with a test function and integrate:

u*e*dx = div(q)*e*dx

Where u is a trial function in DG0, e is a test function in DG and q is in whereever.

But, as David points out, the DG0 mass-matrix is completely diagonal. The values on the diagonal are just obtained by assembling the lone test function:

diag = assemble(e*dx)

But now, because the matrix is diagonal, you can solve the linear system pointwise by doing a pointwise division:

rhs = assemble(div(q)*e*dx)

solution = assemble(rhs/diag)

Cheers,

Lawrence

_______________________________________________ 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