Simple Heat Conduction

Lately, I was trying to get the flow running with H3D. Since it did not work (as usual), I was looking where the problem might be. So, I constructed a simple heat conduction problem:

where is chosen in such a way that

You can look at the resulting video.

Note: the purpose of this was to verify that the sequence of calculated solutions was correct.

Metis Support

In 3D, we have to deal with large domains. One way how to deal with it is domain decomposition. One package that is used for partitioning domains is Metis. I'm glad to say that we have the support of Metis in Hermes3D now. It was very easy to add it. Package itself compiled with not problems. Writing a thin layer around it was a question for couple of hours C++ coding. Very nice thing is that Metis provides standalone binaries that one can execute in some script (for example) and also a librbary that one can link and call it without running any external programs and struggling with paths, etc.

So how it looks in Hermes3D:

  Metis metis;
  metis.partition(mesh, n_partitions);            // partition the mesh into n_partition parts
  int *eparts = metis.get_element_partitions();   // get the element -> partition mapping

Of course, one has to enable Metis support by setting WITH_METIS in global CMake.vars.

For now, the code sits in my personal repository, since we do not have the other neccessary parts needed for doing the whole thing. However, this is the first step in this direction. And, of course, here are some pictures...

Call Stack Dump

When we were running some calculations using Hermes3D, it might happen that something went wrong and the library crashed. We could get segmentation fault for example. At this moment, we had no info what went wrong. If we were running debug version and some assert failed, we got the line of the code, that caused the problem but no additional info.

So I implemented call stack tracing (currently pushed into the main git repo). Now, if the program itself or the library crashes, users are provided with a callstack, so you can trace what exactly happened.

For your own programs, you need to #include < common/callstack.h > and put _F_ macro at the begining of each function that should appear in the call stack dump. When developing inside the Hermes3D, you should include the _F_ macros for newly created functions. Current ones are already patched.

Note: the call stack dump is printed to stderr.

Tests Time Length

I have recently ran the full test suite of Hermes3D and found out, that some hanging node related tests failed. So I started to dig into the code to find the problem. At first, I though the problem is in the baselist merging, where I did some patches, but continutity test of all base functions was ok. So I dug deeper in the code and discovered, that the right-hand side vector is corrupted. The problem was, that when I was adding Matrix and Vector data classes to support wider range of problems, I did not initialized the right hand side vector, thus the contributions were added to some garbagge.

Today, I came to the work and checked the test results. All hanging nodes tests are ok. It took 11 hours to run the full test suite (751 test cases - 8 processors @ 3GHz, 32 GB RAM). That is not so bad, I took 3 days on my laptop (dual core CPU @ 2.2GHZ). Still, distributed environment is the way we should go, since then everything should be done in couple of minutes.

Element Orders in VTK

Most of the time for visualization purposes, I use Gmsh, which has built-in post-processing features. The main reason is, that I can visualize the mesh (vertex, edge, face and element numbers) that helps when debugging hanging nodes. The big disadvantage for us is, that GMSH can not index elements from 0 (zero seems to has some special meaning), thus I implemented shift by one of all indices in our code (not very good, but it makes debugging easy, on do not have to think about shifting by one).

Thus, VTK output module was a little bit neglected. Recently, I do have some troubles with GMSH, so I'll probably move to paraview, that uses Vtk. Vtk seems to be very powerful so as paraview (one needs a larger screen for that, my 15" notebook screen is really small and graphic card really weak). There is a mayavi2 package that uses Vtk (mayavi2 is written in python(?) and should provide nice scripting feature to Vtk, which we should be able to use in future). Paraview should also be capable of scripting, but I had no time to look at it.

Another advantage of paraview is, that is should somehow run over the internet connecting. I did not explore that yet, but it if does all the calculation on the remote machine and send the data to visualize, that would be really nic, especially, when our meshes in 3D start to be huge and with high polynomial degrees.

Thus the reason of this blogspot is to announce out_orders() function for Vtk output module and hopefully, this module will thrive more and more.

API Change for Linear Solvers

Up to now, we were able to solve with h3d problems like: a(u,v) = l(v). But to solve an eigenvalue problem, one had to do a workaround which was not very convenient. Our current approach was to instantiate a LinearSolver class, pass it to Discretization class, assemble stiffness matrix and right-hand side, then solve the system and construct the solution (instance of a Solution class). The code could look like:

  UMFPackLinearSolver solver;
  Discretization d(&solver);
  ...
  // create and assemble stiffness matrix and right-hand side
  d.create_stiffness_matrix();
  d.assemble_stiffness_matrix_and_rhs();
  // solve the system
  Solution sln(&mesh);
  bool solved = d.solve_system(1, &sln);

This had two disadvantages:

1. It forced users to use some linear solver.
2. To use an eigenvalue solver, one had to construct two LinearSolver classes, which was stupid.

Thus, to address these issues we:

1. added new data classes to represent sparse matrices and vectors,
2. separated LinearSolver class from the Discretization class.

This greatly extended the modularity of the library. What one has to do now is this (same problem as above):

  // if using UMFPACK 
  UMFPackMatrix mat;   // stiffness matrix
  UMFPackVector rhs;   // right-hand side 

  Discretization d;
  ...
  // create and assemble stiffness matrix and right-hand side
  d.create(&mat, &rhs);
  d.assemble(&mat, &rhs);
  // solve the system
  UMFPackLinearSolver solver(mat, rhs);
  bool solved = solver.solve();
  // construct the Solution class
  Solution sln(&mesh);
  sln.set_space_and_pss(&space, &pss);
  sln.set_solution_vector(solver.get_solution(), false);

Calling create() is important if using sparse matrices - it will precalculate their structure. If using a dense matrix, that call would not be necessary.

All linear solver classes inherit API of LinearSolver class which consist of ony one method solve(). I'm not sure that is the most convenient way how to do stuff, but it works for now. One small problem here is that Solution class is using zeroth index of the solution vector for the Dirichlet DOF, but LinearSolver class do not need it (and do not know about it, of course), so the solution vector of LinearSolver starts with zero, but it corresponds to the first index in solution vector in the Solution class. To avoid the duplication of memory while adding this special coefficient for Dirichlet DOF, all LinearSolvers-derived classes are ment to return a solution vector where the zeroth entry is equal to 1.0 (which is the coefficient for Dirichlet DOF).

One can see that constructing the solution is not that easy as in previous implementation, but overriding the constructor of Solution class should help to squeeze last three lines to just one.

Welcome

Welcome to the Hermes3D blog. Some info about the development will be published here.