Dealii: Python binding

[Rombur]

This is a continuation of the discussion that happened in the PR #2834

@davydden

it depends... Without an intention to start a flame,

I think we always maintain a safe distance to a flamewar, don't we?

the original discussion about python bindings originated in this Google group thread. More specifically what i was thinking about is GUI preprocessor for triangulation and refinement. The current bindings are quite capable in this regard.

May I ask you to sum up the discussion with respect to the python bindings in a separate issue on github? I glimpsed over this discussion with the title "GUI preprocessor" only very briefly - at no point I had the impression we talked about introducing Python bindings.

@Rombur Do you mind if I ask you to also briefly outline the next steps you want to tackle and what the overall goal in terms of features is?

Introducing bindings for everything in the library (assembly, distribution of dofs, etc) is of course good, but my point here is that it is not a must.

But then I thoroughly misunderstood the whole purpose of this PR.

If the current bindings are primarily meant as a means to get a binding to jupyter (or similar pythonic code) in order to control computations interactively, then they should not be part of the library at all but maintained together with the GUI stuff for jupyter.

So, I hope (and this is how I interpreted the PR) that @Rombur 's effort in writing the python bindings is of more general value :smiley:

Python bindings are very fresh and are in the beginning so I don't think a lot of users will immediately jump using them after the next release. In other words, i don't think we need to be constrained in keeping backward compatibility w.r.t. bindings.

This is in essence an argument in the opposite direction: I have encountered numerous deal.II users that used older versions of the library due to _minor_ incompatibilities we introduced. We for example renamed our vector_tools.h and matrix_tools.h headers - already this created major problems (!). If we now tightly couple a deal.II release with the python bindings and happen to change them drastically from one release to the other we force users to use older deal.II versions despite the fact that this is not necessary at all. The current python bindings for example would probably work without any modification with the last 2, or 3 deal.II versions.

That said, I would prefer to relocate any further discussion to a separate issue (We can c&p essential comments made on this PR). :smiley:

@Rombur @davydden A separate issue would have the benefit that we could briefly discuss further work before someone wastes time writing code that has to be redesigned and rewritten another time. Thinking about the last years, I'd say we had this issue at least 3, or 4 times (myself included) where someone started to implement something just to have the discussion about correct design afterwards.

[Rombur]

@tamiko I am going to disappoint you but there is no next step for me. A few months ago, I was working on generating meshes using deal.II functions and I was frustrated with the process not being interactive. I wanted to implement python binding to quickly generate and check a Triangulation. However, I didn't do anything until about a month later when @davydden started to talk about writing a GUI using python (which would be simplified if we could manipulate the Triangulation directly from python). It then also happened that @luca-heltai could use python wrappers to manipulate the Triangulation inside Blender.

For me the only things missing are a way to set boundary ids and a way to deserialize a refined Triangulation to a refined distributed::Triangulation.

[tamiko]

@Rombur

I am going to disappoint you but there is no next step for me. […]

Then I thoroughly misunderstood you. I'm sorry.

For the time being I suggest we keep this issue open - more complete python bindings, or UML support are definitely something we should consider at some point.

[davydden]

thanks @Rombur for writing up and creating an issue. I organise your missing feature into a list so that they won't be buried within the discussion. I will update this comment if there will be further.

Python bindings missing features

• [x] a way to set boundary ids
• [ ] a way to deserialize a refined Triangulation to a refined distributed::Triangulation.

[bangerth]

In case anyone wants to look at such things, I learned yesterday that there are already some deal.II python bindings at
https://github.com/pymor/pymor-deal.II
and described here:
http://arxiv.org/abs/1506.07094

[drwells]

I talked to the pymor guy at PDEsoft: for what its worth, there is some general interest (from him and other people who were at PDEsoft) in making general purpose Python bindings since it allows for people to glue deal.II together with more or less any other scientific library that also has Python bindings.

[dangars]

I wrote a hdf5 wrapper. This improves the integration with python. It is compatible with h5py/numpy. Notably with python3 Unicode strings (h5py UTF8) and complex h5py/numpy numbers. See #7020 and step-62.ipynb

[Rombur]

@dangars very nice!

[agrayver]

Has there been any re-evalutation of this among @dealii/dealii? I use (and like) deal.II in my research, but I do not have a possibility to teach it, because the entering threshold is too high, mostly because of C++. Outside CS departments, you absolutely cannot assume students will know C++ and it is normally not part of their curriculum.

In our department (and I guess this is quite general) there is a strong trend towards integrating Python and its eco-system (JupyterLab, anaconda, etc) into educational process. There are many reasons for this: open-source, reproducibility, demand in industry.

I believe that creating fully functional python bindings for deal.II will greatly facilitate its adoption in communities other than CS and applied math.

What do you guys think about that?

[Rombur]

@agrayver The problem is that it's a lot of work and we don't have time to create fully functional python bindings. I wrote the python bindings for the mesh because it's self contained and you can do something useful with it but I don't see what could be the next step. To have something useful, we still need to port the vectors, the sparse matrices, the dof handler, the finite elements, the solver, and the data output. Each one of these elements alone is useless so we need to be sure that we will be able to do everything before we start to work on any of these. Either someone needs funding to do that, which I doubt will happen, or we need to find a way to split that huge project in smaller ones that we can slowly tackle. I think the way forward would be to do something similar to what @dangars is doing and slowly move functionalities to python and still use C++ for the core (with @dangars PR, we create the mesh and the input file in python, run the C++ code, and analyze the data in python) but we still need to find projects small enough that we can tackle one at the time.

[dangars]

Hi @agrayver. As @Rombur suggested the best approach is to tackle tackle one project at a time.

You can do most of your deal.ii work on a Jupyter notebook, indeed I do most of my work with python3. In order to do this you can use:

• Mesh with python
• HDF5 which is a powerful data exchange data format/library and is the natural way to exchange data between python and a deal.ii C++ program
• Paraview python interface to analyze the vtk/vtu files. It works great!
• I'll show in PR #6747 how to integrate HDF5 and Paraview on a Jupyter Notebook. Here you can find a preview step-62.ipynb
• With this approach you have the flexibility of python and the perfomance of deal.ii C++ :wink:

Before I finishing PR #6747, I want to finish PR #7020. The serial/MPI HDF5 interface is fully functional, although the API may change.

[agrayver]

@dangars @Rombur while I really like your work, I must admit it looks more like a pre- and post-processing with Python. You probably can call it "binding", but it is not what I meant.

My message above was more related to the fact that it is practically impossible to integrate deal.II into teaching process outside CS and math departments, simply because too few students know C++. What added to this is that I see multiple departments in my field which started to offer Python courses and also increasing number of publications, where people implement physics in Python.

Anyway, I understand there is problem with man power here.

[Rombur]

@dangars @Rombur while I really like your work, I must admit it looks more like a pre- and post-processing with Python. You probably can call it "binding", but it is not what I meant.

We do really have python binding for the Triangulation. If you use PyDealII, you really use the C++ object from python (you can see here (at the end of the file) that you can get an error from the C++ code) but of course since you are limited to the Triangulation, you will need to dump the mesh and then, read it from C++. Therefore, like you said, you are limited to pre-processing but the infrastructure is there to extend the bindings to more part of the library if someone wants to.

Anyway, I understand there is problem with man power here.

Yes, unfortunately. I would really love to try an idea using python and then, if it works create the "real" application in C++.

[bangerth]

We all agree that it's not a philosophical objections to Python but a limitation in the number of developers. But I would like to point out that I have been using deal.II in teaching for many years. It's true that few students know C++ well at the beginning, but most learn it passably well and get up to speed pretty quickly. My experience is that in a semester-long course, C++ vs Python is not the biggest obstacle; it is whether students have any kind of programming experience and know how to think like a programmer, not the particular syntax they have to use.

[dangars]

@Rombur I've been thinking about the python bindings. And I think that it should be possible to do python bindings building on the current capabilities of deal.ii.

I think that the requirements for the python bindings are:

• It should be easy to use. Somebody with minimal finite elements knowledge should be able to use it.
• There should not be a performance hit. It should be as fast as a pure C++ deal.ii program
• It should be possible to do refinement.
• The python bindings provide an easy interface to run simple linear simulations (ex elastic, maxwell...).
• It should be possible to extend the C++ part of the python bindings to run other complex simulations such as non-linear.

We already have

• Python bindings for the Triangulation by @Rombur
• HDF5 interface
• Symbolic Algebra by @jppelteret

This is how I would implement it. I would take the approach that I have used in step-62 but with a better structure. Let me know what you think!

• Proof of concept: https://github.com/dealii/dealii/blob/master/examples/step-62/step-62.ipynb
• The python program will be a pre-processor and a post-processor
• The python and the deal.ii programs are completely decoupled.
  
  
  
  • They can run in the same computer or in different computers.
  
  
  • For example, the python program can run in a desktop to prepare the simulation and then the deal.ii program can run in a large HPC system that does not have python.
  
  
  • Using python C bindings in the program that runs the simulation can be tricky. The interaction between the Python Global Interpreter Lock (GIL) and MPI can be hard to manage and if we use MPI4py then MPI will be managed by python. This is why I would prefer to have a clean separation between the python pre/post processor and the simulation deal.ii program.
  
  
• HDF5 is used to exchange data between the python program and the deal.ii program.
• In the python program the user defines:
  
  
  
  • The triangulation
  
  
  • The boundary conditions
  
  
  • The force
  
  
  • The weak form of the equation that has to be to solved. I think that it should be possible to write the weak form of the equation in python and then use Differentiation::SD::Expression in deal.ii to assemble the system.
  
  
  • Specify the solver to use.
  
  
• It should be easy to replace the Differentiation::SD::Expression by a C++ expression.
  
  
  
  • First you write the equation using python, then when you are happy with the results you can replace the Differentiation::SD::Expression by a C++ version.
  
  

If you think that this approach is good then I'm happy to work on this and provide a proof of concept :)

[Rombur]

I like the general idea but:

I think that it should be possible to write the weak form of the equation in python and then use Differentiation::SD::Expression in deal.ii to assemble the system.

I don't know how easy it is to do but if you want to give it a shot I think that would be great.

I think that there are different ways to use python. What you propose is great for prototyping a code or quickly try a different a slightly different equation. The way I usually it is for debugging purpose so it is useful to have a more tight coupling between deal.II and python. However, both cases are valid and I'd be happy to help you in case you need any help.

[jppelteret]

@dangars Its great that you're wanting to work more on the python bindings. I think that the possibilities on this front are quite exciting. I think that its important that, from the C++ side, we don't end up reinventing things that have already been implemented. The symbolic stuff actually has two ongoing projects. Mine is roadmapped here; everything you see in the assembly helper classes (except the GENERIC integrator) is already implemented but just not pushed upstream. That will happen _soon_ (TM)... As for the other project, I'll let @masterleinad decide whether or not he's willing to share what he's been up to.

What's good about these assemblers I've implemented is that once you've concretised the expression to be evaluated then its possible to compile them, so the performance during evaluation increases to (nearish) native speeds. However, I think that important to recognise that the Differentiation::SD::Expression is only good for scalar evaluation for now. This is because it (well, the library its built on) does not support non-associative and non-commutative operations at this point in time. For example, I have a wrapper for linear operators (so, sparse linear algebra) implemented, but its currently paused until the non-commutative operations are in place. So maybe it'd be good to set out a more concrete plan for what you're thinking on the assembly side, so that (1) there's no overlap in existing work and (2) you don't end up spending time working in a direction that might lead predictably (according to our current knowledge and experience) to a dead-end.

So my implementation of the assembly operations actually requires one to express on a per-dof basis the entire expression for the residual, which is automatically linearised, or an energy functional on a per-cell basis from which the residual or linearisation is automatically determined. So in my own work I actually write out an entire assembly operation that mimics what one usually does, but only using symbolic expressions. Later one performs a substitution operation that again requires a set of operations that very much look like those you'd do traditionally in deal.II. For example:

Pre-computed symbolic evaluation:

            std::vector<sd_type> residual_qp(n_dofs_per_cell, sd_type(0.0));

            // Symbolic expression for kinematic fields
            const Tensor<2, dim, sd_type> Grad_u = symb_data.solution_gradient(fe_values, u_fe, "u");
            const Tensor<2, dim, sd_type> F      = Physics::Elasticity::Kinematics::F(Grad_u);
            //        const Tensor<1,dim,sd_type> H = -symb_data.Grad_sp;

            // Express the contribution to each degree-of-freedom of the residual
            // in terms of the field variables that are "local" or natural to the
            // material law itself.
            const SymmetricTensor<2, dim, sd_type> S = lqph_q_point->get_symbolic_S();
            const Tensor<1, dim, sd_type>          B = lqph_q_point->get_symbolic_B();

            // Symbolic expressions for finite element data
            const sd_type                              JxW        = symb_data.JxW();
            const std::vector<Tensor<2, dim, sd_type>> Grad_Nx_u  = symb_data.shape_function_gradients(fe_values, u_fe, "u");
            const std::vector<Tensor<1, dim, sd_type>> Grad_Nx_sp = symb_data.shape_function_gradients(fe_values, sp_fe, "sp");
            Assert(Grad_Nx_u.size() == n_dofs_per_cell, ExcInternalError());
            Assert(Grad_Nx_sp.size() == n_dofs_per_cell, ExcInternalError());

            for (unsigned int I = 0; I < n_dofs_per_cell; ++I)
              {
                Assert((residual_qp[I] == sd_type(0.0)), ExcInternalError());

                // Determine the component and block associated with
                // the I'th local degree-of-freedom.
                const unsigned int block_I = fe.system_to_base_index(I).first.first;

                // We first compute the contributions from the internal forces.
                // Note that we assemble the residual contribution directly
                // as it is this vector that is to be automatically linearised.
                if (block_I == u_block) // u-terms
                  {
                    // Variation of the Green-Lagrange strain tensor
                    const SymmetricTensor<2, dim, sd_type> dE_I = symmetrize(transpose(F) * Grad_Nx_u[I]);
                    residual_qp[I] += (dE_I * S) * JxW;
                  }
                else if (block_I == sp_block)
                  {
                    // Variation of the magnetic field vector
                    const Tensor<1, dim, sd_type> dH_I = -Grad_Nx_sp[I];
                    residual_qp[I] -= (dH_I * B) * JxW;
                  }
                else
                  {
                    Assert(block_I <= sp_block, ExcInternalError());
                  }
              }

            // Construct a map to convert field variables that are "local" to
            // the constitutive law to their finite-element level description
            const sd_type::substitution_map_t sub_vals_local_to_fe =
              make_local_kinematic_to_fe_substitution_map(symb_data, lqph_q_point->get_symbolic_H(), lqph_q_point->get_symbolic_C(), fe_values, u_fe, sp_fe);

            // Register the dependent variables (i.e. the quadrature
            // point contribution to the total cell residual vector)
            symb_data.register_residual_vector(residual_qp, sub_vals_local_to_fe);

            // Construct the symbolic optimisation map, taking the following
            // symbols to be registered as independent variables:
            // - Constitutive law data
            //   - Constitutive parameters
            //   - Local / internal variables
            // - Cell FE data
            //   - Integration constant
            //   - DoF values
            //   - Shape function gradients
            typename sd_type::substitution_map_t sub_vals_optim = lqph_q_point->make_optimisation_symbol_map();
            SD::add_to_symbol_map(sub_vals_optim, make_fe_optimisation_symbol_map(symb_data, fe_values, u_fe, sp_fe));

            // And last of all we perform optimisation of all of the
            // collected symbolic expressions
            symb_data.optimize(sub_vals_optim);

Per-cell value substitution:

    data.reset();

    // Initialise all helper data structure for this cell
    scratch.reinit(cell);
    data.reinit(cell, scratch.hp_fe_values);
    Assert(data.cell_rhs.l2_norm() == 0.0, ExcInternalError());
    Assert(data.cell_matrix.frobenius_norm() == 0.0, ExcInternalError());

    // Define some aliases to make later code a little
    // bit more legible
    const FEValuesExtractors::Vector &u_fe  = scratch.u_fe;
    const FEValuesExtractors::Scalar &sp_fe = scratch.sp_fe;
    // const unsigned int &u_block             = scratch.u_block;
    //     const unsigned int &sp_block            = scratch.sp_block;
    const FEValues<dim> &fe_values = scratch.hp_fe_values.get_present_fe_values();
    //     const FiniteElement<dim> &fe            = fe_values.get_fe();
    const unsigned int &n_q_points      = fe_values.n_quadrature_points;
    const unsigned int &n_dofs_per_cell = fe_values.dofs_per_cell;

    // Retrieve the local data associated with this cell
    const std::vector<std::shared_ptr<const Continuum_Point<dim, number_type>>> lqph = quadrature_point_history.get_data(cell);
    Assert(lqph.size() == n_q_points, ExcInternalError());

    // Extract the local degree-of-freedom values
    // from the solution vector
    std::vector<double> local_dof_values(n_dofs_per_cell);
    for (unsigned int K = 0; K < n_dofs_per_cell; ++K)
      local_dof_values[K] = scratch.solution_total(data.local_dof_indices[K]);

    // Compute solution-related values at integration points
    fe_values[u_fe].get_function_gradients(scratch.solution_total, scratch.solution_grads_u_total);
    fe_values[sp_fe].get_function_gradients(scratch.solution_total, scratch.solution_grads_sp_total);

    // Use the optimiser to compute the residual and its
    // linearisation on this cell
    Vector<double>     residual_qp(n_dofs_per_cell);
    FullMatrix<double> linearisation_qp(n_dofs_per_cell, n_dofs_per_cell);
    for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
      {
        const Continuum_Point<dim, number_type> *const lqph_q_point = lqph[q_point].get();

        // Compute the kinematic values, with which we can update
        // data at the calculation point and ompute kinetic data.
        const Tensor<2, dim, double>                                       F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]);
        const Tensor<1, dim, double>                                       H = -scratch.solution_grads_sp_total[q_point];
        const Constitutive_Laws::Coupled_Material_Values<dim, number_type> values(F, H);

        // Update all data at the quadrature point and extract
        // that which is necessary to do future computations.
        {
          Continuum_Point<dim, number_type> *lqph_q_point_nc = const_cast<Continuum_Point<dim, number_type> *>(lqph[q_point].get());
          lqph_q_point_nc->solve_local_equilibrium(values);
          lqph_q_point_nc->compute_internal_data(values);
        }

        // Use the LQPH data to more quickly retrieve the symbolic data
        const auto it_symbolic_qp_hash = map_lqph_symbolic_qp_hash.find(lqph_q_point);
        Assert(it_symbolic_qp_hash != map_lqph_symbolic_qp_hash.end(), ExcInternalError());
        const std::size_t &hash_qp_energy      = it_symbolic_qp_hash->second;
        const auto         it_symbolic_qp_data = map_hash_symbolic_qp_data.find(hash_qp_energy);
        Assert(it_symbolic_qp_data != map_hash_symbolic_qp_data.end(), ExcInternalError());
        SD_Helper &symb_data = it_symbolic_qp_data->second;

        // Perform symbolic substitution, taking into consideration
        // the following symbols are are registered as independent
        // variables:
        // - Constitutive law data
        //   - Constitutive parameters
        //   - Local / internal variables
        // - Cell FE data
        //   - Integration constant
        //   - DoF values
        //   - Shape function gradients
        typename sd_type::substitution_map_t sub_vals = lqph_q_point->make_symbol_value_map(values);
        SD::add_to_symbol_value_map(sub_vals, make_fe_symbol_value_map(symb_data, local_dof_values, fe_values, q_point, u_fe, sp_fe));
        symb_data.substitute(sub_vals);

        // Accumulate the contribution from this quadrature point
        symb_data.compute_residual(residual_qp);
        symb_data.compute_linearization(linearisation_qp);
        data.cell_rhs -= residual_qp; // RHS = - residual
        data.cell_matrix.add(1.0, linearisation_qp);
      }

From my current perspective, this is sort-of what you'd need the python code to translate to / call. I guess that this is not so easy for multi-field problems but there might be ways to work around this (e.g. write some canonical symbolic assembly loops that are sequentially called to build up a per-cell and per-field expression, which is later linearised). Thoughts?

[dangars]

I think that there are different ways to use python. What you propose is great for prototyping a code or quickly try a different a slightly different equation. The way I usually it is for debugging purpose so it is useful to have a more tight coupling between deal.II and python. However, both cases are valid and I'd be happy to help you in case you need any help.

@Rombur I agree, both cases can be useful.

I think that its important that, from the C++ side, we don't end up reinventing things that have already been implemented.

@jppelteret I completely agree.

What's good about these assemblers I've implemented is that once you've concretised the expression to be evaluated then its possible to compile them, so the performance during evaluation increases to (nearish) native speeds. However, I think that important to recognise that the Differentiation::SD::Expression is only good for scalar evaluation for now.

I will start to play with Differentiation::SD::Expression, although probably I will wait until you and @masterleinad have made further progress. I could prepare a small proof of concept. With the proof of concept we can decide if it is worthy to continue or if it is a direction that might lead to a dead-end.

[jppelteret]

I could prepare a small proof of concept. With the proof of concept we can decide if it is worthy to continue or if it is a direction that might lead to a dead-end.

Sounds good to me. I'll be opening up a PR with a full implementation of the SD assembly mechanism in the next 2 weeks (although this would be merged in a systematic, piece-wise manner) so you won't have to wait too long to play around.

Thanks a lot for all your enthusiasm around this! I think that it'll be a very cool addition and useful feature to have. Everyone likes Juptyer notebooks nowadays :-D

[dangars]

Sounds good to me. I'll be opening up a PR with a full implementation of the SD assembly mechanism in the next 2 weeks (although this would be merged in a systematic, piece-wise manner) so you won't have to wait too long to play around.

Sounds good! I'll keep an eye in your PRs :)

[masterleinad]

I'll let @masterleinad decide whether or not he's willing to share what he's been up to.

I have been working on a domain-specific language for different backends, but in particular for MatrixFree, in C++. You can find some material from the last workshop here. This is likely a bit orthogonal to what you are trying to do here, though.

[dangars]

I have been working on a domain-specific language for different backends, but in particular for MatrixFree, in C++. You can find some material from the last workshop here. This is likely a bit orthogonal to what you are trying to do here, though.

I took a look to your presentation, maybe Differentiation::SD::Expression is more in the direction of what I'm trying to do. Although I'll still keep an eye also in your PRs!

[drwells]

@agrayver I think that your work on the python bindings has addressed nearly everything outlined above (except 'a way to deserialize a refined Triangulation to a refined distributed::Triangulation'). Do you think that we can now close this issue?

[agrayver]

@drwells probably yes, although I have not done much on the actual FE side (DoFHandler, FiniteElement), but I plan to do that. My long-term goal is to use python + deal.II for the teaching, which will require that.

For cross-reference, #9015 is the issue that summarizes current efforts in this direction.

[drwells]

The remaining feature seems very specialized - I don't think we can do this even in the main library yet. Since the current state of the bindings is tracked elsewhere lets close this.

[drwells]

Like the other old issues I have closed recently, if anyone believes that we should keep this open then please do so.

[agrayver]

except 'a way to deserialize a refined Triangulation to a refined distributed::Triangulation'

BTW, this is not true. You can create a refined p::d::t from a deserialized serial triangulation by using this constructor:

https://www.dealii.org/current/doxygen/deal.II/classparallel_1_1distributed_1_1Triangulation.html#ae09eb5f25fd81eca041851bd8eedec15