dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

On 26 September 2011 15:21, Andre Massing <massing@xxxxxxxxx> wrote:
> On 09/26/2011 01:15 PM, Anders Logg wrote:
>> On Mon, Sep 26, 2011 at 12:06:55AM -0700, Johan Hake wrote:
>>> On Sunday September 25 2011 23:54:57 Anders Logg wrote:
>>>> As part of assembling systems on overlapping meshes (for Nitsche type
>>>> methods), my student Andre Massing is calling the regular DOLFIN
>>>> assemble function for part of the domain and then calling a special
>>>> purpose assembler (currently residing in a branch called dolfin-olm)
>>>> for terms that are not supported by the regular assembler.
>>>>
>>>> For this to work out well, I'd like to propose the following changes
>>>> to the assembler:
>>>>
>>>> 1. Make the assemble_cells, assemble_exterior/interior_facets public
>>>>
>>>> We need to access these separately.
>>>>
>>>> 2. Add options for controlling if and how we call A.apply()
>>>
>>> Wouldn't this be unnesseary if 1 is implemented? Then you just impement your
>>> own more flexible assemble routine?
>>
>> Good point.
>>
>> Andre: Would that be enough? Note that apply() is not called from the
>> assemble_cells/interior/exterior_facets functions, only from the main
>> (currently public) assemble function.
>
> In principle that would be okay. But it enforces the user
> to reimplement (or basically to copy paste) the entire
>
> void Assembler::assemble(GenericTensor& A,
>                         const Form& a,
>                         const MeshFunction<uint>* cell_domains,
>                         const MeshFunction<uint>* exterior_facet_domains,
>                         const MeshFunction<uint>* interior_facet_domains,
>                         bool reset_sparsity,
>                         bool add_values)
>
> even if the user only wants to first do an regular assemble (via that
> function) and then add/insert manually new entries.
> Top copy paste the entire assemble function just to get rid of one apply
> call at the end seems not to be really user friendly.
>
> One nice example is the use of the ident_zeros function which as of
> today is "damm" slow :) if there are many zeros along the diagonal.
> That is because the Matrix will be called with an apply() at the end
> of the assemble routine and hence adding values afterwards via
> ident_zeros is very expensive.

My understanding is that it does not add entries. Someone did add some
hacks which I never understood the need for. A good solution is to not
use Assembler but SystemAssembler. The latter should really be the
default in DOLFIN since it preserves symmetry.

I haven't tested, but I don't know if, for PETSc, MatZeroRows can be
called post set/insert but pre apply.

The function Assembler::assemble referred to above could have
something like 'bool finalise_matrix=true' added to the function
signature.

Garth

> Naturally a user want to assemble,  ident_zeros and then call apply.
> Basically copy pasting the the whole assemble function
> only to omit the postpone the one apply call at the end of assemble
> function seems to be an overkill.
>
> So one thing is the possibility to omit the apply at the end the
> assemble. The other thing is to make assemble_cells, etc. public or
> protected to ease the implementation of new assemblers. As you said
> there are related, but I can think of simple scenarios where  a
> completely new assembler is not needed but just a tiny function to add
> some elements. There, the missing option of postponing the apply makes
> it rather cumbersome to reuse the nice functionality which already is
> provided in DOLFIN.
>
> --
> Andre
>
>
>
>
>
>
>
>
>>
>> --
>> Anders
>>
>>
>>> Johan
>>>
>>>> A.apply() calls PETSc MatAssemblyBegin/End with MAT_FINAL_ASSEMBLY.
>>>> This makes it *very* expensive for us since we need to switch between
>>>> inserting and adding values to the matrix. We need to instead use the
>>>> option MAT_FLUSH_ASSEMBLY. It would further be useful to have an
>>>> option of not calling A.apply() at all.
>>>>
>>>> I'm not sure how this should be handled in the GenericTensor interface
>>>> since Trilinos does not seem to have a corresponding option.
>>>>
>>>> A first step would be to just add a boolean option to the assembler:
>>>>
>>>>   call_apply=false
>>>>
>>>> Then we could manually flush the values in between calls.
>>>>
>>>> Opinions?
>>>>
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help   : https://help.launchpad.net/ListHelp
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>