Hermes3D weblog

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.

Powered by ikiwiki.