DOLFIN

Segregated assembly of weak forms seems to be in trouble

Asked by Luis Linares

Hello,

I have wanted to experiment with the assembly process, and the following gets actually done:

wf_a = -fr*(dot(Er, vr)+dot(Ei, vi))*dx + fi*(dot(Ei, vr) - dot(Er, vi))*dx + dot(curl(Er), curl(vr))*dx + dot(curl(Ei), curl(vi))*dx
ds = Measure('ds')[boundaries]
n_v = FacetNormal(mesh)
d_exc = boundary_list['radiate']
r_v = Expression(['x[{0:d}]'.format(i) for i in range(3)], element=V_c.ufl_element())
k_dir_c = Expression(['{0:g}'.format(k_dir[i]) for i in range(3)], element=V_c.ufl_element() )
E0_c = Expression(['{0:g}'.format(E0[i]) for i in range(3)], element=V_c.ufl_element() )
wf_a += dot(k0*cross(n_v, cross(n_v, Ei)), vr)*ds(d_exc)
wf_L = k0*dot(sin( k0*dot(k_dir_c, r_v) )*cross(n_v, cross(E0_c, k_dir_c-n_v)), vr)*ds(d_exc)
wf_a -= dot(k0*cross(n_v, cross(n_v, Er)), vi)*ds(d_exc)
wf_L += k0*dot(cos( k0*dot(k_dir_c, r_v) )*cross(n_v, cross(E0_c, k_dir_c-n_v)), vi)*ds(d_exc)
A_all, b_rhs = assemble_system(wf_a, wf_L, exterior_facet_domains=boundaries)

whereas the following gets stuck in "Calling FFC just-in-time (JIT) compiler, this may take some time." while assembling the second matrix (it has been running for 10 minutes now and counting...):

wf_all = dot(curl(Er), curl(vr))*dx + dot(curl(Ei), curl(vi))*dx
wf_direc = -fr*(dot(Er, vr) + dot(Ei, vi))*dx
wf_cross = fi*(dot(Ei, vr) - dot(Er, vi))*dx
ds = Measure('ds')[boundaries]
n_v = FacetNormal(mesh)
d_exc = boundary_list['radiate']
r_v = Expression(['x[{0:d}]'.format(i) for i in range(3)], element=V_c.ufl_element())
k_dir_c = Expression(['{0:g}'.format(k_dir[i]) for i in range(3)], element=V_c.ufl_element() )
E0_c = Expression(['{0:g}'.format(E0[i]) for i in range(3)], element=V_c.ufl_element() )
wf_all += dot(k0*cross(n_v, cross(n_v, Ei)), vr)*ds(d_exc)
wf_rhs = k0*dot(sin( k0*dot(k_dir_c, r_v) )*cross(n_v, cross(E0_c, k_dir_c-n_v)), vr)*ds(d_exc)
wf_all -= dot(k0*cross(n_v, cross(n_v, Er)), vi)*ds(d_exc)
wf_rhs += k0*dot(cos( k0*dot(k_dir_c, r_v) )*cross(n_v, cross(E0_c, k_dir_c-n_v)), vi)*ds(d_exc)
A_all, b = assemble_system(wf_all, wf_rhs, exterior_facet_domains=boundaries)
A_direc = assemble(wf_direc)
A_cross = assemble(wf_cross)
A_full = A_all + A_direc + A_cross

For context, I list the declaration of relevant entities:

element_type = "Nedelec 1st kind H(curl)"
element_order = 2
V_c = FunctionSpace(mesh, element_type, element_order)
V = MixedFunctionSpace([V_c, V_c]) #Real part, imag. part
Ef = TrialFunctions(V)
vf = TestFunctions(V)
Er, Ei = Ef[0], Ef[1]
vr, vi = vf[0], vf[1]
# Mark boundaries
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
# Radiation boundaries
class RadiationBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
radiation_bnd = RadiationBoundary()
radiation_bnd.mark(boundaries, 1)
# Weak form Coefficient
fr = Function(FunctionSpace(mesh, "DG", 0))
fi = Function(FunctionSpace(mesh, "DG", 0))

Thank you!

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:

This question was reopened

Revision history for this message
Luis Linares (luis-linares) said : 		#1

Well, I have overcome the issue. Both codes had

parameters["form_compiler"]["optimize"] = True

The non-working code is now working by removing this line (I guess, effectively setting the option to False).

Revision history for this message
Luis Linares (luis-linares) said : 		#2

After sorting this out, I have faced another issue:

The non-segregated, working version was able to be solved with an LUSolver with a reasonable memory consumption. However, the segregated version is not. The LUSolver takes more memory (I guess, way more memory) than my system has and so the process gets killed. Is there anything I can do to remedy this? Perhaps direct addition of matrices is not supported?

Revision history for this message
Martin Sandve Alnæs (martinal) said : 		#3

Regarding the initial issue, can you post a full working script with the form that takes forever to compile? I need real work examples of complicated forms for profiling some new and faster form compiler algorithms.

Revision history for this message
Luis Linares (luis-linares) said : 		#4

Martin,

I wasn't sure about the best way to "post" the files. So I made a tar archive that you can download at this address:

http://db.tt/w6H1xdvL

It includes a parameter and mesh files, all in XML, and the two python scripts, one that works with the "optimize" option asserted, and the other that fails, everything else being equal.

Still, if there is anybody who can shine some light on the issue of LU solution of the segregated-then-added system matrix, it would be of invaluable help. Thank you.

Revision history for this message
Johan Hake (johan-hake) said : 		#5

On 06/17/2012 06:01 AM, Luis Linares wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Status: Open => Solved
>
> Luis Linares confirmed that the question is solved:
> Well, I have overcome the issue. Both codes had
>
> parameters["form_compiler"]["optimize"] = True
>
> The non-working code is now working by removing this line (I guess,
> effectively setting the option to False).
>

This is probably related to:

   https://bugs.launchpad.net/ffc/+bug/956979

The FFC optimization procedure can be time consuming. This is being
worked on.

Also you might want to try SFC (SyFi form compiler) instead of FFC. SFC
lack some features FFC supports (like internal facet integrals), but can
be more efficient on large and complex forms. Install sfc and put:

   parameters["form_compiler"]["name"] = "sfc"

at the top of your script.

Johan

Revision history for this message
Johan Hake (johan-hake) said : 		#6

On 06/17/2012 06:36 PM, Luis Linares wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Status: Solved => Open
>
> Luis Linares is still having a problem:
> After sorting this out, I have faced another issue:
>
> The non-segregated, working version was able to be solved with an
> LUSolver with a reasonable memory consumption. However, the segregated
> version is not. The LUSolver takes more memory (I guess, way more
> memory) than my system has and so the process gets killed. Is there
> anything I can do to remedy this?

What platform are you using and what linear algebra backend are you using.

> Perhaps direct addition of matrices is not supported?

what do you mean with this?

johan

Revision history for this message
Joachim Haga (jobh) said : 		#7

>
>
> > Perhaps direct addition of matrices is not supported?
>
> what do you mean with this?
>
>
Adding matrices together. It should work fine I think, as long as the
sparsity pattern matches. It usually does when using the same test/trial
functions.

I think there may be trouble when a form only uses a subdomain, but from a
quick look it seems like all your forms are global. So that shouldn't be
the case here.

If not, you may use the .axpy interface instead, where you can supply a
flag same_nonzero_pattern=False to indicate that they are different.

-j.

Revision history for this message
Johan Hake (johan-hake) said : 		#8

On 06/18/2012 09:56 AM, Joachim Haga wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Joachim Haga proposed the following answer:
>>
>>
>>> Perhaps direct addition of matrices is not supported?
>>
>> what do you mean with this?
>>
>>
> Adding matrices together. It should work fine I think, as long as the
> sparsity pattern matches. It usually does when using the same test/trial
> functions.
>
> I think there may be trouble when a form only uses a subdomain, but from a
> quick look it seems like all your forms are global. So that shouldn't be
> the case here.

The sparsity pattern is build using the functionspace of the trial and
test functions. subdomains should not play any part in that.

> If not, you may use the .axpy interface instead, where you can supply a
> flag same_nonzero_pattern=False to indicate that they are different.

same_nonzero_pattern is set to false if you use the included iadd and
add operators. If you know the matrices share the same sparsity patter,
which should be the case for you, you should use axpy when adding
matrices and same_nonzero_pattern=True. I think this flag is only used
for the PETSc backend.

Johan

>
> -j.
>

Revision history for this message
Joachim Haga (jobh) said : 		#9

>
> The sparsity pattern is build using the functionspace of the trial and
> test functions. subdomains should not play any part in that.
>

See https://bugs.launchpad.net/dolfin/+bug/927539 for a pathological case...

Revision history for this message
Johan Hake (johan-hake) said : 		#10

On 06/18/2012 11:46 AM, Joachim Haga wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Joachim Haga proposed the following answer:
>>
>> The sparsity pattern is build using the functionspace of the trial and
>> test functions. subdomains should not play any part in that.
>>
>
> See https://bugs.launchpad.net/dolfin/+bug/927539 for a pathological
> case...

Ok, sorry to misslead you.

Johan

Revision history for this message
Luis Linares (luis-linares) said : 		#11

Thank you all for your suggestions!

Johan, I have built the dev version of fenics, but I didn't find a repository for the SyFi source, as it does not appear on
https://launchpad.net/fenics-project

Where can I find this? Thank you.

Revision history for this message
Luis Linares (luis-linares) said : 		#12

Well I found it, but it is not present on the fenics-project list. Has it been unsuspectingly dropped from the list?

Revision history for this message
Johan Hake (johan-hake) said : 		#13

On 06/18/2012 05:06 PM, Luis Linares wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Luis Linares posted a new comment:
> Well I found it, but it is not present on the fenics-project list. Has
> it been unsuspectingly dropped from the list?

That is a bug. Thanks for reporting.

Someone with admin rights need to fix that :)

Johan

Revision history for this message
Luis Linares (luis-linares) said : 		#14

Johan

After building SyFi I have recompiled dolfin, ufc, ufl and yet it complains:

No module named sfc
*** Warning: Could not import sfc form compiler, falling back to FFC.
Calling FFC just-in-time (JIT) compiler, this may take some time.

Any tips for this issue?

Revision history for this message
Luis Linares (luis-linares) said : 		#15

Ok I have solved the sfc module issue by sourcing the file syfi.conf (I would have thought the install step would have taken care of that, since mine is a system-wide installation...). But now the ultimate issue is:

  File "/usr/local/lib/python2.7/site-packages/sfc/representation/elementrepresentation.py", line 216, in __init__
    self.syfi_element = create_syfi_element(self.ufl_element, self.polygon, default_order=self.options.code.finite_element.default_order_of_element)
  File "/usr/local/lib/python2.7/site-packages/sfc/representation/elementrepresentation.py", line 78, in create_syfi_element
    raise NotImplementedError("Not implemented element family '%s'." % f)

NotImplementedError: Not implemented element family 'Nedelec 1st kind H(curl)'.

So I'm out of luck here. But thanks alot anyway for the support!

Revision history for this message
Luis Linares (luis-linares) said : 		#16

Johan,

On the issue of adding matrices, I confirm that using axpy works, with either value of same_nonzero_pattern. Thanks a lot for the suggestion!

Revision history for this message
Johan Hake (johan-hake) said : 		#17

On 06/18/2012 08:31 PM, Luis Linares wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Luis Linares posted a new comment:
> Ok I have solved the sfc module issue by sourcing the file syfi.conf (I
> would have thought the install step would have taken care of that, since
> mine is a system-wide installation...). But now the ultimate issue is:
>
> File "/usr/local/lib/python2.7/site-packages/sfc/representation/elementrepresentation.py", line 216, in __init__
> self.syfi_element = create_syfi_element(self.ufl_element, self.polygon, default_order=self.options.code.finite_element.default_order_of_element)
> File "/usr/local/lib/python2.7/site-packages/sfc/representation/elementrepresentation.py", line 78, in create_syfi_element
> raise NotImplementedError("Not implemented element family '%s'." % f)
>
> NotImplementedError: Not implemented element family 'Nedelec 1st kind
> H(curl)'.

Oups! Sorry to misslead you. I thought these elements were supported.

> So I'm out of luck here. But thanks alot anyway for the support!

You're welcome

Johan

Revision history for this message
Anders Logg (logg) said : 		#18

On Mon, Jun 18, 2012 at 04:41:27PM -0000, Johan Hake wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Johan Hake proposed the following answer:
> On 06/18/2012 05:06 PM, Luis Linares wrote:
> > Question #200665 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/200665
> >
> > Luis Linares posted a new comment:
> > Well I found it, but it is not present on the fenics-project list. Has
> > it been unsuspectingly dropped from the list?
>
> That is a bug. Thanks for reporting.
>
> Someone with admin rights need to fix that :)

Which list?

--
Anders

Revision history for this message
Johannes Ring (johannr) said : 		#19

On Tue, Jun 19, 2012 at 12:25 PM, Anders Logg
<email address hidden> wrote:
> Which list?

I think it was this list:

  https://launchpad.net/fenics-project

I have added SyFi to the list now.

Johannes

Revision history for this message
Anders Logg (logg) said : 		#20

On Tue, Jun 19, 2012 at 10:36:01AM -0000, Johannes Ring wrote:
> Question #200665 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200665
>
> Johannes Ring proposed the following answer:
> On Tue, Jun 19, 2012 at 12:25 PM, Anders Logg
> <email address hidden> wrote:
> > Which list?
>
> I think it was this list:
>
> https://launchpad.net/fenics-project
>
> I have added SyFi to the list now.

ok, great.

--
Anders

Can you help with this problem?

Provide an answer of your own, or ask Luis Linares for more information if necessary.

To post a message you must log in.