5-th Feel++ User Days

Past, Present and Future

Introduction

Versatile

Gallery feelpp 600x600

A large range of numerical methods to solve partial differential equations: cG, dG, hdG, crb,…​ in 1D, 2D and 3D

Powerful

supercomputer 1 600x600

Support for high performance computing up to thousands of cores for linear and non-linear problems using PETSc/SLEPc and InHouse solution strategies

Expressive

feelpp dsel 600x600

A language for Galerkin methods embedded into C++ for maximal mathematical expressivity: meshes, function spaces and elements, bilinear and linear forms, algebraic representation and post-processing…​

Toolboxes : Monophysics

CFD(left), CSM(middle), Heat transfer(right)

FlowAroundCylinder 400x200

torsionbarNeoHookIncompT2 400x200

heat transfer building 400x200

Toolboxes : Multiphysics

FSI(left), aerothermal(middle), thermoelectric(right)

wp3dP3P2G2 struct disp t2 400x200

cabine 400x200

peltiermodule electricpotential 400x200

A Wide Range of Numerical Methods

feelpp methods

Projects

Feelpp project 2017

Partners

AcademicEntreprises

CNRS/LNCMI

Plastic Omnium

Unistra/ICUBE

Airbus

Unistra/Inserm

EDF

Liphy

Socomec

CNR Pavia (It)

Atos (Spain)

Politecnico di Milano(It)

Holo3

U. Missouri (US)

See-d

Feel++ web ecosystem

http://www.feelpp.org

http://book.feelpp.org

http://github.com/feelpp/feelpp

http://twitter.feelpp.org

http://youtube.feelpp.org

http://plus.feelpp.org

http://publications.feelpp.org

Feel++ Team

prudhomm icon

vincentchabannes icon

szopos icon

guidoboni icon

bertoluzza icon

sacco icon

trophime icon

ismail icon

metivet icon

sala icon

hild icon

gdolle icon

wahl icon

francoisdh icon

ricka icon

Installation (simplified)

Docker

we are going to use Docker to demonstrate some Feel++ features during the next 4 days.

Installation

$ docker run feelpp/feelpp-toolboxes echo "Hello World!" (1)
Console Output
Unable to find image 'feelpp/feelpp-toolboxes:latest' locally  (1)
latest: Pulling from feelpp/feelpp-toolboxes
8e21f82d32cf: Pull complete
[...]
0a8dee947f9b: Pull complete
Digest: sha256:457539dbd781594eccd4ddf26a7aefdf08a2fff9dbeb1f601a22d9e7e3761fbc
Status: Downloaded newer image for feelpp/feelpp-toolboxes:latest
Hello World! (2)
More details on Docker during the Infrastructure talk.

From Source

  • Required: C++14 compiler, CMake, MPI(open-mpi or mpich) and boost-1.61.

  • Recommended: PETSc, SLEPc, HDF5, Gmsh{2,3}, Python,OPENTurns

  • Suggested: Paraview, Ensight, OpenModelica, Octave, MeshGems

Feel++ full build

git clone https://github.com/feelpp/feelpp.git
git submodule init (1)
git submodule update --recursive
mkdir build && cd build (2)
../configure -r (3)
1recent cleanup in submodules
2in-source builds forbidden
3configure an helper script (see --help)

Feel++ staged builds

# build only quickstart if Feel++ libs installed
<toplevel feel++ srcdir>/configure -r --root=<toplevel feel++ srcdir>/quickstart
make
# build only crb if Feel++ libs installed
<toplevel feel++ srcdir>/configure -r --root=<toplevel feel++ srcdir>/applications/crb
make
# build only toolboxes if Feel++ libs installed
<toplevel feel++ srcdir>/configure -r --root=<toplevel feel++ srcdir>/toolboxes
make

And for Feel++ projects/applications

#similarly for Feel++ projects
<toplevel feel++ srcdir>/configure -r --root=<project top srcdir>
<toplevel feel++ srcdir>/configure -r --root=<eye2brain top srcdir>
make

Feel++ dependencies

C++

  • We depend on a C++14 compiler with a C++14 standard C++ library: GCC (>=6) clang (>=3.4)

  • C++17 has just been formally approved. Many nice features like variable template folding, parallel STL, Structured binding declarations. Checkout C++17 on Wikipedia

  • GCC 7 is out

  • Clang 5 is out

Boost

Minimum (1.61) and Maximum(1.63) versions supported
New Boost Library : Boost.hana
  • C++14 compiler + C++14 stdlib

  • see product spaces, forms over product spaces, block vector/matrix

  • more efficient than Boost.fusion (see benchmarks)

Boost.hana

template<typename... SpaceList>
class ProductSpaces : public  hana::tuple<SpaceList...>, ProductSpacesBase
{
    using super = hana::tuple<SpaceList...>;
    using functionspace_type = ProductSpaces<SpaceList...>;
    ProductSpaces( SpaceList... l ) : super( l... ){}
    constexpr int numberOfSpaces() const { return hana::size( *this ); }

    //! \return the total number of degrees of freedom
    auto ndof = [&](size_type s, auto& e ) { return s + e->nDof(); };
    size_type nDof() const { return hana::fold_left( *this, 0, ndof  ); }
    //! \return the number of degrees of freedom owned by the process
    auto nlocaldof = [&](size_type s, auto& e ) { return s + e->nLocalDof(); };
    size_type nLocalDof() const { return hana::fold_left( *this, 0, nlocaldof ); }
...

Boost

From Boost to C++ std equivalent
  • boost::shared_ptr<> vs std::shared_ptr<> (C++14)

  • boost::optional<> vs std::optional<> (C++17)

  • others?

Keeping Boost provides (may be) more flexibility

Eigen

  • Tensor Library (used by TensorFlow)

  • Eigen objects within the algebraic core of Feel++

    • Geometric Mapping

    • `PolynomialSet`s data structure : interpolation, context…​

  • Eigen::Map is used for the transition from Boost::ublas

Large update upcoming in this area

GFLags/GLog

feelpp_qs_laplacian --config-file laplacian/triangle/triangle.cfg --no_log=0|1|2 (1)
1controls the generation of logfiles.
--no_log=[0,1,2]
  • 0 logfiles for all MPI processes,

  • 1(default) logfile only processor 0 and

  • 2 no logfile at all

See Best practices session

Feel++ Core

Multithreading and MPI

Multi-threading (multiple) enabled by default
  • MPI Multi-threading level support

    • single only one thread

    • funneled only main thread make MPI calls

    • serialized Only one thread at the time do MPI calls

    • multiple Multiple thread may do MPI calls.

  • Used in

    • Static condensation (See Lorenzo/Romain’s talk)

std::async

#include <future>
#include <feel/feel.hpp>
int main(int argc, char** argv )
{
  using namespace Feel;
  Environment env( _argc=argc, _argv=argv,
                   _treading=threading::multiple);
  auto mesh=loadMesh(_mesh=new Mesh<Simplex<2>>);
  auto i1 = std::async( std::launch::async,
    [&]() { return integrate(_range=elements(mesh), _expr=cst(1.0)).evaluate(); } );
  auto i2 = std::async( std::launch::async,
    [&]() { return integrate(_range=boundaryfaces(mesh),_expr=cst(1.0)).evaluate();} );
  std::cout << ”i1=” << i1.get() << ” i2=” << i2.get() << ”\n”;
}

Checker

feelpp_qs_laplacian.cpp
// compute l2 and h1 norm of u-u_h where u=solution
auto norms = [=]( std::string const& solution ) ->std::map<std::string,double>
{
   tic();
   double l2 = normL2(_range=elements(mesh), _expr=idv(u)-expr(solution) );
   toc("L2 error norm");
   tic();
   double h1 = normH1(_range=elements(mesh), _expr=idv(u)-expr(solution),
                                             _grad_expr=gradv(u)-grad<2>(expr(solution)) );
   toc("H1 error norm");
   return { { "L2", l2 }, {  "H1", h1 } };
};
int status = checker().runOnce( norms, rate::hp( mesh->hMax(), Vh->fe()->order() ) );

Checker

  • Check quickstart for other examples

  • Simple, extensible, json based

{
   "x^2+y^2:x:y":{ "2": { "exact": true } },
   "sin(pi*x)*cos(pi*y):x:y":{
       "2": {
           "h": ["0.2","0.1","0.05"],
           "errors": {
               "L2": {
                  "order": 3,
                  "values": ["0.00245063","0.000340925","4.20359e-05"]
               },
               "H1": {
                   "order": 2,
                   "values": ["0.043385","0.0114881","0.00284834" ]
} } } } }

Checker

$ feelpp_qs_laplacian_2d --config-file laplacian/circle/circle.cfg
Console Output
$ ctest -R  feelpp_qs_laplacian_2d-circle-dirichlet-np-6
1: Test command: /usr/bin/mpiexec "-np" "6" "/home/feelpp/build/feelpp_qs_laplacian_2d" "--config-file" "/home/feelpp/src/feelpp/quickstart/laplacian/circle/circle-dirichlet.cfg" "--checker.tolerance.exact=5e-14"
[....]
1: Reading /home/feelpp/src/feelpp/quickstart/laplacian/circle/circle-dirichlet.cfg...
1: [ Starting Feel++ ] application qs_laplacian version 0.104.0-alpha.30 date 2017-Sep-10
1:  . qs_laplacian files are stored in /feel/qs_laplacian/circle-dirichlet/np_6
1:  .. exports :/feel/qs_laplacian/circle-dirichlet/np_6/exports
1:  .. logfiles :/feel/qs_laplacian/circle-dirichlet/np_6/logs
1: [loadMesh] Loading mesh in format geo+msh: "/home/feelpp/src/feelpp/quickstart/laplacian/circle/circle-dirichlet.geo"
1: [loadMesh] Use default geo desc: /home/feelpp/src/feelpp/quickstart/laplacian/circle/circle-dirichlet.geo 0.1
1: [loadMesh] Time : 0.0470395s
1: [Vh] Time : 0.00543124s
1: [l] Time : 0.0104013s
1: [a] Time : 0.025204s
1: [a.solve] Time : 0.0100693s
[...]
[exact verification success]||u-u_h_H1=1.30953e-14

TESTS and .tests.*

  • ctest is a powerful tool

  • Combined with Checker and we have a powerful combination to toroughly test our work

  • we want to check as many scenario as fast as possible without making changes to .cfg files

solution: use TESTS in feelpp_add_application and define a .tests.<yourapp>

TESTS and .tests.*

CMakeLists.txt
feelpp_add_application( laplacian_2d SRCS qs_laplacian.cpp ... TESTS)
.tests.laplacian_2d
triangle-x+y --config-file ${CMAKE_CURRENT_SOURCE_DIR}/laplacian/triangle/triangle.cfg   --functions.g=x+y:x:y --checker.solution=x+y:x:y
triangle-x^2+y^2 --config-file ${CMAKE_CURRENT_SOURCE_DIR}/laplacian/triangle/triangle.cfg --functions.g=x^2+y^2:x:y --checker.solution=x^2+y^2:x:y

TESTS and .tests.*

$ ctest -R laplacian
Test project /home/feelpp/build
      Start  1: feelpp_qs_laplacian_2d-circle-dirichlet-np-6
 1/30 Test  #1: feelpp_qs_laplacian_2d-circle-dirichlet-np-6 ...............   Passed    1.00 sec
      Start  2: feelpp_qs_laplacian_2d-circle-dirichlet-np-1
 2/30 Test  #2: feelpp_qs_laplacian_2d-circle-dirichlet-np-1 ...............   Passed    0.95 sec
      Start  3: feelpp_qs_laplacian_2d-circle-robin-np-6
.....
 29/30 Test #29: feelpp_qs_laplacian_3d-tetrahedron-quadratic-np-6 ..........   Passed    2.25 sec
      Start 30: feelpp_qs_laplacian_3d-tetrahedron-quadratic-np-1
30/30 Test #30: feelpp_qs_laplacian_3d-tetrahedron-quadratic-np-1 ..........   Passed    3.41 sec

100% tests passed, 0 tests failed out of 30

Total Test time (real) =  43.84 sec

Timers

feelpp_qs_laplacian_2d --config-file laplacian/triangle/triangle.cfg --display-stats=1
Console output
Timer                                                                Count    Total(s)      Max(s)      Min(s)     Mean(s)   StdDev(s)
  Backend::newMatrix:: build stencil                                     1    5.43e-03    5.43e-03    5.43e-03    5.43e-03    0.00e+00
  Backend::newMatrix:: initialize matrix                                 1    1.02e-03    1.02e-03    1.02e-03    1.02e-03    0.00e+00
  Backend::newMatrix:: zero out matrix + set split                       1    8.11e-05    8.11e-05    8.11e-05    8.11e-05    0.00e+00
      DofTable::build                                                    3    6.49e-03    3.96e-03    1.22e-04    2.16e-03    1.58e-03
[snip]
Timer                                                                Count    Total(s)      Max(s)      Min(s)     Mean(s)   StdDev(s)
 Vh                                                                       1    2.60e-01    2.60e-01    2.60e-01    2.60e-01    0.00e+00
  H1 error norm                                                          1    2.51e-01    2.51e-01    2.51e-01    2.51e-01    0.00e+00
 loadMesh                                                                 1    4.43e-02    4.43e-02    4.43e-02    4.43e-02    0.00e+00
 a                                                                        1    4.12e-02    4.12e-02    4.12e-02    4.12e-02    0.00e+00
 l                                                                        1    2.49e-02    2.49e-02    2.49e-02    2.49e-02    0.00e+00
  Exporter                                                               1    2.34e-02    2.34e-02    2.34e-02    2.34e-02    0.00e+00
 a.solve                                                                  1    1.17e-02    1.17e-02    1.17e-02    1.17e-02    0.00e+00
    Timeset::add uh                                                      1    8.25e-03    8.25e-03    8.25e-03    8.25e-03    0.00e+00
[snip]
use tic()/toc("<label>",print=true|false) to add new timer <label> in the table

Mesh

Mesh in Feel++

  • Supports 1D to 3D including curves and surfaces

  • Supports simplex and hypercube

  • Supports high order mesh

  • File formats :

    • geo format of GMSH

    • all mesh file formats supported by GMSH : msh, mesh, vtk, med, …​

    • Feel++ mesh file format : hdf5+json

Submesh

From a mesh of dim d and a range of entities, create a submesh by extracting:

  • Elements (submesh of dim d)

  • Faces (submesh of dim d-1)

  • Edges in 3d (submesh of dim d-2)

auto submesh = createSubmesh(mesh,markedelements("toto"));
auto submesh = createSubmesh(mesh,markedfaces("toto"));
auto submesh = createSubmesh(mesh,markededges("toto"));

Submesh

the partitioning is preserved in a submesh, which means that there might be load-balancing issues.

submeshinitialmeshpartitioning submeshpartitioning

Create a sequential submesh from a parallel mesh

auto submesh = createSubmesh(mesh,markedelements("toto"),
                             Environment::worldCommSeq());

Mesh markers

Feel++ had the possibility to use 3 kinds of marker : physical, marker2, marker3.

Now, we can define an arbitrary number of marker, defined by a tag (unsigned int)

for ( auto const& eltWrap : elements(mesh) )
{
   auto & eltModified = mesh->element(unwrap_ref( eltWrap ).id()):
   if ( eltModified.id() % 2 == 0 )
       eltModified.setMarker( 36, 3);
   else
       eltModified.setMarker( 36, 5);
}
the tag 1 is reserved for physical markers !

Iterating over marked entities

auto range = markedelementsByType(mesh,36,3);
markedelements(..), marked2elements(..) and marked3elements(..) are still available.

Mesh filter: range from expression

// select region from vertex coordinates
auto range = elements(mesh, pow(Px()-0.4,2)+pow(Py()-0.5,2) < pow(cst(0.23),2) );
// the same with Ginac
auto range = elements(mesh, expr( "(x-0.4)^2+(y-0.5)^2 < 0.23^2:x:y" );

Mesh Filter: concatenation of ranges

Let \(r_1\) and \(r_2\) be two range sets, compute \(r=r_1 \cup r_2\).

auto range=concatenate(r1,r2);


Let \((r_i)_{i=1,\ldots,N}\) be \(N\) range sets, compute \(r= \cup_{i=1}^N r_i\).

auto range=concatenate(r1,r2,...,rN);

Mesh filter: intersection of ranges

Let \(r_1\) and \(r_2\) be two range sets, compute \(r=r_1 \cap r_2\).

auto range=intersect(r1,r2);


Let \((r_i)_{i=1,\ldots,N}\) be \(N\) range sets, compute \(r= \cap_{i=1}^N r_i\).

auto range=intersect(r1,r2,...,rN);

Mesh Filter: range from entity ids

auto range=idelements(mesh,{2,35,67});


Mesh filters: an example

auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>> );
auto r1 = elements(mesh,pow(Px()-0.4,2)+pow(Py()-0.5,2) < pow(cst(0.23),2));
auto r2 = elements(mesh,pow(Px()-0.7,2)+pow(Py()-0.5,2) < pow(cst(0.15),2));
auto r3 = concatenate(r1,r2);
auto r4 = intersect(r1,r2);

auto Vh = Pdh<0>(mesh);
auto u = Vh->element(), v = Vh->element(), w = Vh->element();
u.on(_range=r3,_expr=cst(1.));
v.on(_range=r4,_expr=cst(1.));
w = u + v;

auto e = exporter( _mesh=mesh );
e->add( "u", u );
e->add( "v", v );
e->add( "w", w );
e->save();

Mesh filters: an example

meshfiltersconcatenante meshfiltersintersect meshfiltersconcatenanteintersect

Mesh filters: filter return type

bool useBoundary = true;
auto therange = (useBoundary)? boundaryfaces(mesh):internalfaces(mesh);
All mesh filters by entity return the same type

Range of entities

using elt_t = MeshTraits<MeshType>::elements_reference_wrapper_type;
auto myelts = boost::make_shared<elt_t>();

for( auto const& eltWrap : elements(mesh) )
{
     auto const& elt = unwrap_ref( eltWrap);
     if ( elt.id() % 2 == 0 )
          myelts->push_back( elt );
}
auto therange = boost::make_tuple(mpl::size_t<MESH_ELEMENTS>(),
                                  myelts->begin(),myelts->end(),
                                  myelts );
Perspectives

Definition of range types should be simplified

Mesh partitioner

  • Partitioning of mesh into Feel++ mesh format : hdf5+json

  • Adapted for the HPC

  • Reloading mesh from this format is faster than gmsh file (even in binary), and particularly in // with large number of process

Mesh partitioner Example

$ feelpp_mesh_partitioner --dim=2 --ifile toto.msh --odir mypartitions
                          --part 4 16 32 48

--part allows to generate a set of partitioning, here we generate 4 different partitioning with 4, 16, 32 and 48 partitions respectively

partitioning2d

aorta partitions32

peltiermodule partitioning

Partitioning by markers

$ feelpp_mesh_partitioner --dim=2 --ifile toto.msh --odir mypartitions
                          --part 4 16 32 48
                          --by-markers

option --by-markers enables partitioning by marker id

For each physical marker, generate a partitioning. If we have \(p\) processors and \(m\) markers then we have elements of each marker on \(p\) processors

partitioning2dbymarkergeo

partitioning2dbymarker

Partitioning by markers

Specify a list of markers to partition together or separately

  • : allows to split marked region partitioning

  • , allows to list marked region partitioned together

  feelpp_mesh_partitioner --dim=2 --ifile toto.msh --odir mypartitions
                          --part 4 16 32 48 --by-markers-desc mat1:mat2,mat3

partitioning2dbymarkergeo

partitioning2dbymarkerdesc

Reduce Mesh loading time

  • Objective: Load rapidly a 3d mesh with 200 millions of tethraedron

  • Loading time take a very long time

    • Internal data structure of mesh : boost::multi-index

    • Insertion of new element : re-ordering of all iterators (very expansive!)

  • Very high memory footprint

Faster mesh data structure

  • Internal data structure : replace boost::multi-index by std::unordered_map (container using hash-table)

    • Complexity : Average case: constant, Worst case: linear in container size.

    • References to elements in the container remains valid in all cases, even after a rehash.

  • Use a std::vector of reference to keep ordering (and have deterministic algorithm)

  • Improve algorithms used in mesh build thanks to new data structure

Faster mesh data structure

This type of code is now broken

for ( auto const& elt : elements(mesh) )
{
   double eltSize = elt.h();
}

it should now be written

for ( auto const& eltWrap : elements(mesh) )
{
   auto const& elt = unwrap_ref( eltWrap ):
   double eltSize = elt.h();
}

Reduce memory footprint

  • Do not store some quantities on each element but recompute on the fly : h(), hMin(), hMax(), G(), vertices(), hFaces()

  • Store some measures on each element only if used

  • Add possibility to build partialy sub data structure in the mesh (without faces not marked, no faces, no edges, no cache)

  • Add possibility to store nodes coordinates in float precision instead of double precision (warning, must be use with precaution!)

Reduce memory footprint next steps

  • Index type as a template argument to handle precision

    • Currently it’s 64 bit, 32 bit by default should be enough for all cases treated currently.

    • The memory footprint reduction gain should be important.

Vectors and Matrices

Algebraic representations

// Let Vh be a function space, e.g.
Vh=Pch<1>(mesh);
auto F = backend(_kind=<backend> )->newVector( Vh ); (1)
auto A = backend()->newMatrix( _test=Vh, _trial=Vh );
1Available <backend> : petsc : sequential and parallel, eigen : sequential .

Support various algebra operators

Files formats : ascii, binary (and soon hdf5)

with PETSc, the close() call is not necessary, expect after local assembly (e.g setValue(…​))

Static condensation

Three new classes
  • VectorCondensed<T>: stores local vector contributions with a dictionary and block wise

  • MatrixCondensed<T>: stores local matrix contributions with a dictionary and block wise

  • StaticCondensation<>: executes the static condensation within standard HDG and HDG/IBC, see Lorenzo/Romain talks

*Condensed<> inherit from *Block classes

Solver/Preconditioner

Solver/Preconditioner

Function space

Parallel Dof Table

  • Step 1 : build the local (sequential) dof table (including ghosts)

  • Step 2 : update the global (parallel) dof table at interprocess

partition dofghosts2

Parallel Dof Table

  • Step 3 (optional) : the extended dof table

extendeddoftablepartition
  • Discontinuous Galerkin method

  • Hybride Galerkin method

  • CIP Stabilisation

Parallel Dof Table

  • Step 4 : reordering of degrees of freedom

    • All active dofs before

    • All ghost dofs after

crutial step for the conversion of vectors

Parallel Dof Table

Take into account interprocess connection by a geometric point, edge or face.

helicespartitioning

Domain definition

These updates allow to build a functionspace on a range of mesh element, not necessary on all elements of the mesh.

Feature

Avoid building a submesh to reduce the memory footprint

Use Case: multiphysics applications

Each physic are define in different region that may or may not overlap

Domain definition

auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>> );
auto rangeElt = elements(mesh,expr(soption(_name="functions.a")));
auto Vh = Pch<2>( mesh,rangeElt );
auto u = Vh->element("u"), v = Vh->element("v");

auto l = form1( _test=Vh );
l = integrate( _range=rangeElt, _expr=id(v) );
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=rangeElt, _expr=gradt(u)*trans(grad(v)) );

Vh->dof()->meshSupport()->updateBoundaryInternalFaces();
auto rangeBoundaryFaces = Vh->dof()->meshSupport()->rangeBoundaryFaces();

a+=on(_range=rangeBoundaryFaces, _rhs=l, _element=u, _expr=cst(0.));
a.solve(_rhs=l,_solution=u);

auto e = exporter( _mesh=mesh );
e->add( "u", u ); e->save();

Domain definition

$ myprog --gmsh.hsize=0.01
  --functions.a="((x-0.4)^2+(y-0.5)^2<0.23^2)+((x-0.7)^2+(y-0.5)^2<0.15^2):x:y"

spaceonrangelaplacian

spaceonrangelaplacianwrap

Domain definition

Use in Toolbox

A physics section has been added in material property. It allow to define several physics with respect to the materials

Possible values in the fluid toolbox json specification file are :

  • aerothermal

  • fluid

  • heat-transfert

Domain definition

Example: geometry(left), velocity(middle), temperature(right)

spaceonrangetest2dgeo

spaceonrangetest2dvelocity

spaceonrangetest2dtemperature

Domain definition

Aerothermal/thermo-electric/heat-transfert
nste temp int11

Product of Spaces: same mesh

auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
int n = 10; // product of 10 function spaces
auto X = ProductSpace<decltype(Pch<1>(mesh))>( n, mesh );
auto U = X.element(); //  element of X
auto u = U[0];      // view on element in product
auto ex = exporter(_mesh=mesh);
for(int i = 0; i < n; ++i )
{  ex->add(nameU[i],U[i]); }
ex->save();
build the space resulting from the product of \(n\) spaces of the same type

Product of Spaces : same instance of same mesh type

  • only one function space(and dof table) is actually built internally

  • many applications such as handling chemicals or different species

Different meshes of the same type

auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
int n = 10; // product of 10 function spaces
std::vector<mesh_ptr_t<Simplex<2>>> meshes ( n );
for( int i : irange( 0, n ) )
{ meshes[i] = createSubmesh( mesh, markedelements(mesh, std::to_string(i+1) ) ); }

auto X = ProductSpace<decltype(Pch<1>(mesh))>( meshes );
auto U = X.element(); //  element of X
// requires n exporters for each mesh or interpolate on mesh
auto ex = exporter(_mesh=mesh);
for(int i = 0; i < n; ++i )
{  ex->add(nameU[i],U[i]); }
ex->save();

Element of a function space

auto Vh = Pch<2>( mesh );
auto v = Vh->element();
auto Wh = THch<1>( mesh );
auto U = Wh->element();
auto u = U.element<0>(); // view
auto p = U.element<1>(); // view
auto ux = u[Component::X] // view
Very powerful framework, important to understand this!

Element of a function space

  • Algebraic representations : a vector

  • Internal storage : boost uBLAS

  • Support views : component, sub-element

  • Support various algebraic operators dot(), max(), …​

  • Files formats : ascii, binary and hdf5 (recommended)

Vector Conversions between uBLAS and PETSc

PETSc to uBLAS
auto backendPetsc = backend(_kind="petsc");
auto v_petsc = backendPetsc->newVector( Vh );
auto v_ublas = Vh->element( v_petsc );
uBLAS to PETSc
auto v_ublas = Vh->element();
auto v_petsc = toPETScPtr( v_ublas );
No copy, share internal storage

Parallel synchronisation

auto v = Vh->element();
// do some operations on v
sync( v, "=" );
sync( v, "+" );
sync( v, "min" );
sync( v, "max" );
sync assigns the dofs values shared by several proc, it requires communication
Current operators:
  • = : value equal to one of them (default active dof)

  • + : sum all values for each shared dof

  • min/max : min or max value for each shared dof

Parallel synchronisation

If the vector v is modified only locally (e.g. on a marker), the values synchronized with operator "=" can be incorrect. The active dof is not necessarily defined on this marker.
This is the case with strong Dirichlet boundary conditions in Newton method

Parallel synchronisation

Solution : specify the correct dof values
v.on(_range=markedfaces(mesh,"toto"),_expr=cst(3.1415));
// get dof used
std::set<size_type> dofUsed;
for ( auto const& faceWrap : markedfaces(mesh,"toto") )
{
   auto const& face = unwrap_ref( faceWrap );
   auto facedof = Xh->dof()->faceLocalDof( face.id() );
   for ( auto it= facedof.first, en= facedof.second ; it!=en;++it )
      dofUsed.insert( it->index() );
}
// apply sync
sync(v,"=",dofUsed );

Variational formulation and integration

Forms on product spaces

  • build (bi)linear forms on product spaces of

    • the same type (dynamic)

    • different types (static)

    • a mix of same and different types

  • construction of a block wise algebraic representation, including

    • stencil (sparsity allocation)

    • block information for PETSc (fieldsplit)

Forms on product space (same space)

auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
int n = int(doption("parameters.n"));
auto ps = ProductSpace<decltype(Pch<1>(mesh))>( n, mesh );
auto U = ps.element();
auto u = U[0];
auto b = blockform2( ps );
auto l = blockform1( ps );
for( int i = 0; i < n; ++i )
{
   b(i,i) += integrate( _range=elements(mesh), _expr=idt(u)*id(u));
   l(i) += integrate( _range=elements(mesh), _expr=expr(soption("functions."+alphabet[i]))*id(u));
}
b.solve( _rhs=l, _solution=U );
auto ex = exporter(_mesh=mesh);
for(int i = 0; i < n; ++i ) {
  ex->add(std::to_string(i),U[i]);
}
ex->save();

Forms on product space

Different spaces
auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
auto Xh=Pch<1>(mesh); auto Wh=Pchv<2>(mesh);
auto u = Xh->element(); auto v = Wh->element();
auto bbf = blockform2( ps );
bbf( 0_c, 0_c ) = integrate( _range=elements(mesh), _expr=id(u)*idt(u) );
bbf( 0_c, 1_c ) += integrate( _range=elements(mesh), _expr=id(u)*(trans(idt(v))*one() ) ) ;
bbf( 1_c, 0_c ) += integrate( _range=elements(mesh), _expr=1./2*(trans(id(v))*one())*idt(u) );
bbf( 1_c, 1_c ) += integrate( _range=elements(mesh), _expr=trans(id(v))*idt(v) );
bbf.close();
auto blf = blockform1( ps );
blf( 0_c) = integrate( _range=elements(mesh), _expr=id(u) );
blf( 1_c) += integrate( _range=elements(mesh), _expr=trans(id(v))*one() );
blf.close();
auto U=ps.element();
bbf.solve( _solution=U, _rhs=blf );

Forms on product space

Different spaces and same spaces
auto p = product( ps, Xh, Yh, Zh );
auto U = p.element();
auto u = U(0_c);
auto w = U(1_c);
auto z = U(2_c);
auto l = blockform1( p );
auto b = blockform2( p );
l(0_c) = integrate( _range=elements(mesh), _expr=expr(soption("functions."+alphabet[0]))*id(u));
b(0_c,0_c) += integrate( _range=elements(mesh), _expr=idt(u)*id(u));
l(1_c) = integrate( _range=elements(mesh), _expr=expr(soption("functions."+alphabet[1]))*trans(id(w))*one());
b(1_c,1_c) += integrate( _range=elements(mesh), _expr=trans(idt(w))*id(w));
l(2_c) = integrate( _range=elements(mesh), _expr=expr(soption("functions."+alphabet[0]))*id(z));
b(2_c,2_c) += integrate( _range=elements(mesh), _expr=idt(z)*id(z));
for( int i = 0; i < n; ++i )
{
   auto v = U(3_c,i);
   b(3_c,3_c,i,i) += integrate( _range=elements(mesh), _expr=idt(v)*id(v));
   l(3_c,i) += integrate( _range=elements(mesh), _expr=expr(soption("functions."+alphabet[2+i]))*id(v));
}                                                                                                                                             b.solve( _rhs=l, _solution=U );

Lambda expressions

test_vf_integrals.cpp (See Romain and JB talks)
auto mesh = loadMesh( _mesh=new Mesh<Simplex<3>> );
auto Xh = Pchv<2>( mesh );
auto u = Xh->element();
u.on(_range=elements(mesh),_expr=vec(cst(1.),cst(1.),cst(1.)));
Eigen::Matrix<double,3,1> n; n << 1,2,3;
std::vector<Eigen::Matrix<double,3,1>> x( ioption("N"),  n);
tic();
// batch
auto v = integrate(_range=elements(mesh),_expr=_e1v,_quad=_Q<1>()).template evaluate<double,3,1>( x );
toc("lambda integral _e1=vec(x,y,z)",FLAGS_v>0);
// one by one
for( auto e : x )
{
    tic();
    auto vv = integrate(_range=elements(mesh),_expr=vec(cst(e(0)),cst(e(1)),cst(e(2))),_quad=_Q<1>()).evaluate(false);
    toc("integral vec(x,y,z)",FLAGS_v>0);
}
Use in EIM and Biosavart

Future

Next PR

Ongoing Developments in Feel++ core
  1. Geomap/PolynomialSet update (internal change)

  2. Dynamic quadratures feature/quad (will break code)

  3. Fix finite element framework for any order feelpp/ho (internal change)

  4. Dynamic Polynomial Order (api change, existing code should compile)