1. FEniCS HPC Avatar FEniCS HPC
  2. FEniCS HPC
  3. unicorn

Wiki

Clone wiki

unicorn / Home

View History

Software distribution

For automatic installation we provide the following Unicorn/FEniCS-HPC Singularity container image:

singularity pull shub://sregistry.srv.cesga.es/mso4sc/fenics-unicorn:latest

There is also an experimental Docker image available here for e.g. Amazon EC2:

https://hub.docker.com/r/sputnikcem/fenicshpc/

you can use the following lines on tegner.pdc.kth.se to test the floating platform code within an interactive allocation:

cd /cfs/klemming/nobackup/${USER::1}/${USER}
module add singularity/2.5.1
module add gcc/5.1
module add openmpi/1.10-gcc-5.1
singularity pull shub://sregistry.srv.cesga.es/mso4sc/fenics-unicorn:latest
git clone https://ncde@bitbucket.org/fenics-hpc/unicorn.git
cd unicorn
git fetch --all
git checkout unicorn_var_den_ale
salloc -A 2018-3-206 -N 1 -n 32 -t 01:00:00
mpirun -n 1 singularity exec -B $PWD:/home --pwd /home ../mso4sc-fenics-unicorn-latest.simg make UNICORNLIBDIR=/home
mpirun -n 32 singularity exec -B $PWD:/home --pwd /home ../mso4sc-fenics-unicorn-latest.simg ./demo >  log1 2>log2
mpirun -n 1 singularity exec -B $PWD:/home --pwd /home ../mso4sc-fenics-unicorn-latest.simg /usr/local/bin/dolfin_post -M mesh -t vtk -n 40 -s velocity

you can use the following lines on tegner.pdc.kth.se to test the flow past wing code within an interactive allocation:

cd /cfs/klemming/nobackup/${USER::1}/${USER}
module add singularity/2.5.1
module add gcc/5.1
module add openmpi/1.10-gcc-5.1
singularity pull shub://sregistry.srv.cesga.es/mso4sc/fenics-unicorn:latest
git clone https://ncde@bitbucket.org/fenics-hpc/unicorn.git
cd unicorn
git fetch --all
git checkout next
cd wing_sim01
cp mesh0.bin mesh.bin
salloc -A 2018-3-206 -N 1 -n 32 -t 01:00:00
mpirun -n 1 singularity exec -B $PWD/../:/home --pwd /home/wing_sim01 ../../mso4sc-fenics-unicorn-latest.simg make  UNICORNLIBDIR=/home
mpirun -n 32 singularity exec -B $PWD/../:/home --pwd /home/wing_sim01 ../../mso4sc-fenics-unicorn-latest.simg ./wing >  log1 2>log2
mpirun -n 1 singularity exec -B $PWD/../:/home --pwd /home/wing_sim01 ../../mso4sc-fenics-unicorn-latest.simg /usr/local/bin/dolfin_post -m mesh_out.bin -t vtk -n 200 -s velocity

Validation cases

At the moment, two validation cases are implemented, the surface-mounted cube for viscous flow and the NACA0012 airfoil for high-Reynolds flow. In both cases we perform adaptive computation of the drag force comparing to reference values from the literature.

The drag is computed as

\begin{equation*} F_D = \int_\Gamma \sigma(u,p)(n,e_x) \end{equation*}

where \(\sigma(u,p) = 2 \nu \epsilon(u) - pI\) is the stress tensor, \(\Gamma\) is the surface of integration and \(e_x\) is the flow direction. In case of high-Reynolds flow we apply the approximation \(\nu = 0\).

Surface mounted cube

This test case reproduces the test case presented in the paper Hoffman, J. (2005). Computation of Mean Drag for Bluff Body Problems Using Adaptive DNS/LES. SIAM Journal on Scientific Computing, 27(1), 184–207. https://doi.org/10.1137/040614463.

The present test case involves the computation of the drag coefficient for a surface mounted cube at Reynolds approximately 40'000, and the drag coefficient is reported to be about 1.5 from literature sources. The code used for this simulation is in subfolder smCube of the source tree. The results for the code at revision e1108c3 were

iter_00/log1:drag coefficient: 0.61
iter_01/log1:drag coefficient: 0.504
iter_02/log1:drag coefficient: 0.576
iter_03/log1:drag coefficient: 0.688
iter_04/log1:drag coefficient: 0.75
iter_05/log1:drag coefficient: 0.676
iter_06/log1:drag coefficient: 0.7
iter_07/log1:drag coefficient: 0.88
iter_08/log1:drag coefficient: 1.032
iter_09/log1:drag coefficient: 1.152
iter_10/log1:drag coefficient: 1.248
iter_11/log1:drag coefficient: 1.332
iter_12/log1:drag coefficient: 1.398
iter_13/log1:drag coefficient: 1.442
iter_14/log1:drag coefficient: 1.446
iter_15/log1:drag coefficient: 1.538
iter_16/log1:drag coefficient: 1.582
iter_17/log1:drag coefficient: 1.544

Use the following command to get the output above

grep -a "^mean drag:" iter_*/log1 | gawk '/drag/ { print $1 " drag coefficient " $3*200 } '

NACA0012 airfoil

Commit c39f006 introduced the validation case of the flow past a NACA0012 airfoil with angle of attack of 10º. The code is implemented in subfolder wing of the source tree.

After the updates to the computation of the weak normal introduced in commit 3cd8ae6, for the adaptive computation of the drag on the airfoil the following results were produced

iter_00/log1:mean drag coefficient 0.0180741
iter_01/log1:mean drag coefficient 0.0148148
iter_02/log1:mean drag coefficient 0.0136296
iter_03/log1:mean drag coefficient 0.0126667
iter_04/log1:mean drag coefficient 0.0122963
iter_05/log1:mean drag coefficient 0.0122963
iter_06/log1:mean drag coefficient 0.0123704
iter_07/log1:mean drag coefficient 0.010963
iter_08/log1:mean drag coefficient 0.0111852
iter_09/log1:mean drag coefficient 0.0111852

The reference value for this test case is reported in the Computational Technology Laboratory Report KTH-CTL-4023, titled Simulation of 3D unsteady incompressible flow past a NACA 0012 wing section and is reported in Table 2 on page 8. The reported values are 0.0119 and 0.0117. A third value of 0.0190 is reported that corresponds to a different test case with "wraparound" [but I don't know what it means]. The relative error with the above results is between 5 and 6 percent.

Use the following command to get the output above

grep -a "^mean drag:" iter_*/log1 | gawk '/drag/ { print $1 " drag coefficient " $3*2/(2.7*cos(10.0/360*2*pi)) } '

Weak normals

Method

With class NodeNormal we used to compute the node normal via a weighted average of the normals at the facets sharing a corner. The latest implementation, in commit 3cd8ae6, defines the node normal as the solution \(u_n\) to the variational problem

\begin{equation*} \int_\Omega u_n v_n + \frac{100}{h}\int_{\partial \Omega} (u_n - n) v_n = 0 \quad \forall v_n \end{equation*}

where \(h\) is the cell size, \(u_n\) is a P1 trial function, \(v_n\) is a P1 test function and n = FacetNormal(cell) is UFL's piecewise constant facet normal vector field. The first term is interpreted as a small perturbation of the problem of finding the projection of the normal defined on the boundary \(\partial\Omega\) of the domain. This is needed in order to generate a nonsingular linear system.

Results

Running the cube demo with the original implementation [reference] and the one implemented in the first version [new, commit e7a4828] gives the following results for the first eight iterations

Reference
iter_00/log1:mean drag: 0.464
iter_01/log1:mean drag: 0.438
iter_02/log1:mean drag: 0.406
iter_03/log1:mean drag: 0.417
iter_04/log1:mean drag: 0.417
iter_05/log1:mean drag: 0.417
iter_06/log1:mean drag: 0.417
iter_07/log1:mean drag: 0.417

New
iter_00/log1:mean drag: 0.429
iter_01/log1:mean drag: 0.422
iter_02/log1:mean drag: 0.399
iter_03/log1:mean drag: 0.414
iter_04/log1:mean drag: 0.414
iter_05/log1:mean drag: 0.414
iter_06/log1:mean drag: 0.414
iter_07/log1:mean drag: 0.414

where the difference is actually tiny.

Notes

  1. The variational problem above is a working one, but not necessarily the best one: others are welcome to be tried. One clear disadvantage is that a variational problem, however trivial, needs to be solved on the whole mesh even though we only need information that is defined on the boundary.

Boundaries of complex geometries

Many simulations have to be performed on geometries with very complex shapes and different boundary conditions need to be imposed on different parts of the boundary. In some special cases it is possible to define a boundary by subclassing SubDomain and overriding its virtual bool inside method to provide an analytical expression of the boundary, as in

class InflowBoundary : public SubDomain
{
  bool inside(const real* x, bool on_boundary) const
  {
    return x[0] <= xmin + DOLFIN_EPS and on_boundary;
  }
};

However, most of the times it is impossible, or not worth the effort, to identify a boundary in this way. The use of boundary markers is a possibility, but it is not reliable due to the fact that, by default, the mesh partitioner used by Unicorn is allowed to change the numbering of the elements to optimize parallel performance. The boundary markers usually store element indexes with the number they have in the mesh saved on disk, which might not correspond to the actual number after the mesh is read.

For this reason, the best option is to load a mesh of the boundary and use that to define a SubDomain instance to be passed to a DirichletBC class constructor. Suppose you want to simulate the flow past an airplane and want to impose a boundary condition on the airplane. You will have a mesh Domain.xml of the whole computational domain and a surface mesh Airplane.xml of the airplane surface; then you can define a SubDomain instance as follows:

Mesh mesh("Domain.xml");

dolfin_set("Mesh read in serial", true);

Mesh airplaneMesh("Airplane.xml");

dolfin_set("Mesh read in serial", false);

SubDomain airplane(airplaneMesh);

After the above snippet is executed, the airplane variable can now be used to define a Dirichlet boundary condition.

Notes:

  1. The two dolfin_set instructions are needed as of the moment of writing (2017-11-02) because of a limitation of the current design, which will be improved in the near future.
  2. The mesh Airplane.xml should be a submesh of the mesh Domain.xml. In other words, this procedure will likely not work if the triangles in the surface mesh are not the same triangles as in the mesh of the whole computational domain. Most meshing tools allow to export submeshes so this is not really a limitation.
  3. This method uses GTS, the GNU Triangulated Surface library, to perform the subdomain creation. In order to use this method you will then need to compile dolfin-hpc with the --with-gts flag.

Template policies

Due to the increase in the variety of uses of unicorn in recent times, it was necessary to refactor the code so that it would be more versatile and modular, as not all applications require all the functionalities unicorn offers. Also, we are moving in different research directions that require to explore a number of different models and implementations; having separate branches or, worse, repositories for each of them is bad practice, so now we have one code to rule them all.

The user is now able to choose at compile time what kind of solver he or she needs by selecting the correct template policies. Let's see an example.

Previously [before ecf830f], instantiating and running the solver was done like so

NSESolver solver(
    mesh,
    dbcs_m,
    sbcs_m,
    dbcs_c,
    dbcs_dm,
    &u_i0,
    outputQuantities,
    &sm,
    &psim,
    &psic,
    &bpsim,
    &f,
    cfl_target
    );

solver.run();

The new syntax is the following

typedef AdaptiveSolver<StrongSlip,
                       GlobalTimeStep,
                       UnitVelocityStabilization> NSESolver;
NSESolver solver(
    mesh,
    dbcs_m,
    dbcs_c,
    dbcs_dm,
    sbcs_m,
    &u_i0,
    outputQuantities,
    &psim,
    &psic,
    &bpsim,
    &f,
    cfl_target
    );

solver.run();

The main change is the typedef instruction, that allows you to select the features [or policies] of the solver you want to use. Currently there are three available policies, defined in three source files. The example above reproduces the default solver before this code refactoring was performed. The current policies are

  1. SlipPolicy: this can be NoSlip if you do not want to use slip boundary conditions, StrongSlip for the usual enforcement of slip conditions in the linear system and WeakSlip to have slip conditions enforced in the weak form.
  2. TimeStepPolicy: can be ConstantTimeStep, with obvious meaning, GlobalTimeStep to bound the time step with the minimum cell size and the maximum of the velocity or LocalTimeStep to compute a cell-wise time step bound and use the maximum allowed size.
  3. StabilizationPolicy: can be UnitVelocityStabilization to stabilize the flow assuming that the maximum velocity is 1 [in modulus] or DimensionalVelocityStabilization to actually compute the maximum velocity and stabilize accordingly.

Please note that your choice of policies determines the API of the resulting object. For example, if you choose ConstantTimeStep as your TimeStepPolicy, then you can call solver.setTimeStep(...), but otherwise this will generate a compile-time error. Another thing that changes is the fifth argument to the constructor. Its type is defined as typename SlipPolicy::SlipBCData slipData in the function declaration, which means that you will have to provide a different parameter depending on which SlipPolicy you chose. Check the code for more details; NoSlip::SlipBCData is defined to be void*, so you can pass 0, for example.

This enrichment of the API might look complicated at first, but it is one of the great strengths of template policy designs.

Finally, for those of you who do not want to do an adaptive simulation you can now define your solver to be of type PrimalSolver, like so

typedef PrimalSolver<StrongSlip,
                     GlobalTimeStep,
                     UnitVelocityStabilization> NSESolver;

and the code will run only the primal solver, without creating the many files that are only needed for adaptivity. The syntax is the same at the moment but it is likely to change, as some parameters are not needed for a simple solve.

Updated