Yes. I think that's right. On Tue, 8 Dec 2015 at 17:54 Justin Chang <jychang48@gmail.com> wrote:

Lawrence,

Okay yeah I got it to work, at least for structured grids.

Now if I had unstructured grids, element_meas is not immediately known and varies among elements. In this case, would I use h = CellSize(mesh) and define element_meas = 1/h/h (I am using quadrilaterals for this)?

Thanks, Justin

On Tue, Dec 8, 2015 at 4:17 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:

...

Hi Justin,

On 08/12/15 09:19, Justin Chang wrote:

...

Hi all,

So our research group has this augmented lagrangian solver for the LSFEM in FEniCS. I have attached the code.

While the code works and does what we want, it is extremely slow! I suspect it is because of these lines of code:

1)

a_LMB, L_LMB = {}, {} element_meas = 0.5 / (XSeed * YSeed) for i in range(mesh.num_cells()): a_LMB[i] = (element_meas) * (1./element_meas**2) * div(v) * div(w) * dx(i) L_LMB[i] = (1./element_meas) * div(w) * dx(i)

So you mark each cell as its own subdomain. Then you build a weak form integrating over each cell.

...

...

3)

a, L = aStandard, LStandard for i in range(mesh.num_cells()): a += gamma * a_LMB[i] L += lmbda[i] * L_LMB[i]

Here you build the global weak form:

a = aStandard + \sum_i a_LMB[i]

surely this is just equivalent to saying:

a = a_Standard + gamma * (1/element_meas)*div(v)*div(w)*dx

and similarly:

L = LStandard + (1/element_meas) * lmbda * div(w) * dx

Where lmbda is just a DG0 Function.

...

What the FFC just-in-time compiler does is that it compiles each and every single one of these element-wise forms (e.g., if my mesh has over 1000 elements, the JIT compiler will be called 1000 times). This takes a stupidly long time to run. If I wanted to do something like this in Firedrake, how would I go about making this more efficient?

I think you don't have to define any of these tiny forms (unless I'm really missing something). Moreover, since they're all symbolic, they don't need to be redefined inside the loop over and over either.

...

For 2) I believe the solution is to do something like "Error.assign(assemble(div(v)*e*dx)/assemble(e*dx))" as suggested in a previous thread.

Yes, it looks like this is just the projection of div(mysol.sub(0)) into a DG0 space.

Does that all make sense?

Cheers,

Lawrence

_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake