Is deal.II for me? Or, I got burned elsewhere

[Franco Milicchio]

Dear all, 

after a couple of years developing with Fenics on MacOS X, I stumbled upon two dreary problems: it won't provide binaries anymore, and it does not play well with multithreading. I was then thinking about switching to deal.II, but I hope you can help me in dispelling some doubts.

My analyses are non-linear and parametric, hence my problems with Fenics: I parallelized single analyses, but I found many hidden race conditions.

So, here are my doubts about moving to deal.II:

1 - Is it truly thread-safe? As in running a single analysis, or running several (no data sharing), without problems.

2 - Are there facilities of some kind that can help in generating the assembly code? In Fenics I just specified the weak formulations.

3 - Can I output some kind of time-dependant file, like a PVD? I see VTK, but no PVD files.

4 - I see Trilinos is included, is it possible to avoid completely PETSc/SLEPc, since they're not AFAIK thread-safe?

5 - Is there a multithreaded linear solver and vector operations, or should I write them myself?

6 - Can I use my MLK library without recompiling deal.II?

7 - Is eliminating binary distributions (DMGs and such) on the horizon?

Thank you very much for your help on any of these points!

    Franco

[Bruno Turcksin]

Franco,

I'll try to answer your questions but be aware that deal.II is a C++ library and that there are no python wrappers unlike Fenics:

On Friday, December 2, 2016 at 6:59:16 AM UTC-5, Franco Milicchio wrote:

1 - Is it truly thread-safe? As in running a single analysis, or running several (no data sharing), without problems.

deal.II is thread-safe but PETSc is not so if you are using deal.II with PETSc, your code won't be thread-safe.
 

2 - Are there facilities of some kind that can help in generating the assembly code? In Fenics I just specified the weak formulations.

You need to write the code in C++ and loop over cells, quadrature points, and basis functions yourself. So this will be more verbose than what you do in Fenics.
 

3 - Can I output some kind of time-dependant file, like a PVD? I see VTK, but no PVD files.

Yes, you can http://dealii.org/8.4.1/doxygen/deal.II/classDataOutInterface.html#a0d247c3780325546472c006d8ef3b12d
 

4 - I see Trilinos is included, is it possible to avoid completely PETSc/SLEPc, since they're not AFAIK thread-safe?

Yes, you can. I am not an expert of SLEPc but I think that the anasazi package from Trilinos https://trilinos.org/packages/anasazi has similar capabilities However, I am not sure what is in the Mac binary. I don't know if it comes with Trilinos and if it does, it may come with only the packages required to run deal.II.
 

5 - Is there a multithreaded linear solver and vector operations, or should I write them myself?

The linear solvers and the vector operations are multithreaded if you use deal.II own solvers and vectors. Many operations done inside deal.II are multithreaded but you are in charge of making sure that the code that you are writing is multithreaded. There is an infrastructure to help you with that though.
 

6 - Can I use my MLK library without recompiling deal.II?

No, I don't think that you can.

7 - Is eliminating binary distributions (DMGs and such) on the horizon?

No.

Best,

Bruno

[Franco Milicchio]

Thank you Bruno for your answers, I failed to mention that I use Fenics with C++, not with Python, sorry about that. I just have few points to ask, if you can bear with me.

1 - Is it truly thread-safe? As in running a single analysis, or running several (no data sharing), without problems.

deal.II is thread-safe but PETSc is not so if you are using deal.II with PETSc, your code won't be thread-safe.

So can I use solvers as in Tutorial 3, for instance, and I'll be safe and sound?

 

2 - Are there facilities of some kind that can help in generating the assembly code? In Fenics I just specified the weak formulations.

You need to write the code in C++ and loop over cells, quadrature points, and basis functions yourself. So this will be more verbose than what you do in Fenics. 

I don't mind at all, I started using FEMs way back, when there was no Fenics.

However, is someone planning some additions to deal.II to make it less verbose, as in, for example, providing a DSL for that?

 

4 - I see Trilinos is included, is it possible to avoid completely PETSc/SLEPc, since they're not AFAIK thread-safe?

Yes, you can. I am not an expert of SLEPc but I think that the anasazi package from Trilinos https://trilinos.org/packages/anasazi has similar capabilities However, I am not sure what is in the Mac binary. I don't know if it comes with Trilinos and if it does, it may come with only the packages required to run deal.II.

And just in case, is there a default eigensolver that is shipped for sure with the DMG?

 

5 - Is there a multithreaded linear solver and vector operations, or should I write them myself?

The linear solvers and the vector operations are multithreaded if you use deal.II own solvers and vectors. Many operations done inside deal.II are multithreaded but you are in charge of making sure that the code that you are writing is multithreaded. There is an infrastructure to help you with that though.
 

That is reasonable. This is very good news to me. To what infrastructure are you referring to?

 

7 - Is eliminating binary distributions (DMGs and such) on the horizon?

No.

This is really awesome!

 

Thanks,

    Franco

[Bruno Turcksin]

Franco,

2016-12-02 10:57 GMT-05:00 Franco Milicchio <franco.m...@gmail.com>:
> So can I use solvers as in Tutorial 3, for instance, and I'll be safe and
> sound?

Yes.

> However, is someone planning some additions to deal.II to make it less
> verbose, as in, for example, providing a DSL for that?

There are no plan to add DSL. There has been some work on
LinearOperator that reduces the amount of work that needs to be
written but there is nothing planned in the future.

> And just in case, is there a default eigensolver that is shipped for sure
> with the DMG?

deal.II implements the power method and the inverse iteration method.
I have never used them, though

> That is reasonable. This is very good news to me. To what infrastructure are
> you referring to?

See tutorial 9 http://dealii.org/8.4.1/doxygen/deal.II/step_9.html

Best,

Bruno

[Wolfgang Bangerth]

On 12/02/2016 07:01 AM, Bruno Turcksin wrote:
>
> 2 - Are there facilities of some kind that can help in generating the
> assembly code? In Fenics I just specified the weak formulations.
>
> You need to write the code in C++ and loop over cells, quadrature points, and
> basis functions yourself. So this will be more verbose than what you do in Fenics.

True, but not very much more so. For example, the bilinear form for Stokes is
here:
https://github.com/dealii/dealii/blob/master/examples/step-22/step-22.cc#L649

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Franco Milicchio]

 

There are no plan to add DSL. There has been some work on
LinearOperator that reduces the amount of work that needs to be
written but there is nothing planned in the future.

Ok, I will look into it.

Thanks for all your answers!

    Franco 

[luca.heltai]

On 2 Dec 2016, at 15:01, Bruno Turcksin <bruno.t...@gmail.com> wrote:
>
>
> 2 - Are there facilities of some kind that can help in generating the assembly code? In Fenics I just specified the weak formulations.
> You need to write the code in C++ and loop over cells, quadrature points, and basis functions yourself. So this will be more verbose than what you do in Fenics.
>

However, there are projects that will help you prototype your application (in threadsafe/parallel mode), such as pi-domus

https://github.com/mathLab/pi-DoMUS

Here you can solve a non-linear, time dependent problem, by specifying just the energy of your system (similarly to what is done in the fenics). All system matrices are computed in parallel, using both MPI and TBB, by automatic differentiation.

See, for example, the interface

https://github.com/mathLab/pi-DoMUS/blob/master/include/interfaces/neo_hookean_two_fields.h

Interface files are the only thing needed to have the problem solve non-linear time dependent pdes in pi-domus, e.g.:



template <int dim, int spacedim,typename LAC>
template<typename EnergyType,typename ResidualType>
void NeoHookeanTwoFieldsInterface<dim,spacedim,LAC>::energies_and_residuals(const typename DoFHandler<dim,spacedim>::active_cell_iterator &cell,
FEValuesCache<dim,spacedim> &fe_cache,
std::vector<EnergyType> &energies,
std::vector<std::vector<ResidualType> > &,
bool compute_only_system_terms) const
{
EnergyType alpha = 0;
this->reinit(alpha, cell, fe_cache);

const FEValuesExtractors::Vector displacement(0);

auto &grad_us = fe_cache.get_gradients("solution", "u", displacement, alpha);
auto &Fs = fe_cache.get_deformation_gradients("solution", "Fu", displacement, alpha);

const FEValuesExtractors::Scalar pressure(dim);
auto &ps = fe_cache.get_values("solution","p", pressure, alpha);

auto &JxW = fe_cache.get_JxW_values();

const unsigned int n_q_points = ps.size();

energies[0] = 0;
for (unsigned int q=0; q<n_q_points; ++q)
{
Tensor <1, dim, double> B;

const EnergyType &p = ps[q];
const Tensor <2, dim, EnergyType> &F = Fs[q];
const Tensor<2, dim, EnergyType> C = transpose(F)*F;
auto &grad_u = grad_us[q];

EnergyType Ic = trace(C);
EnergyType J = determinant(F);

EnergyType psi = (mu/2.)*(Ic-dim)-p*(J-1.);
energies[0] += (psi)*JxW[q];
if (!compute_only_system_terms)
energies[1] += (scalar_product(grad_u,grad_u) + 0.5*p*p)*JxW[q];
}
}


with the above 60~ lines of codes you can solve a non-linear incompressible elasticity problem using two fields.

L.

[Franco Milicchio]

Thanks, this is really interesting!

    Franco