I am solving the unsteady heat equation in 3D using an alternating direction implicit (ADI) method. This means that I am solving three different tridiagonal systems within a single timestep (once for each cartesian axis). The size of these matrices is on the order of 10^5 x 10^5.

Although each direction's A matrix has all non-zeros along the three diagonals, there are situations where one direction will have a zero value along one of the diagonals where another direction does not.

My original plan was to set up an AIJ matrix and preallocate for 3 non-zero values and reuse the same memory allocation for each of the three directions. This works fine for solving the first direction, but when I hit the second direction PETSc complains because I am attempting to add a new non-zero value where there wasn't one before (related to differences in boundary conditions).


What is the most efficient way to proceed:

  • Is there a way to allocate differently that will allow me to have values change between zero and non-zero?
  • OR should I create a new matrix and new allocation for each direction? (This I know works, but I'm worried about the overhead of destroying and creating a new matrix for every solution.)
  • Other options?

In pseudo code, my program looks like this:

// Create and preallocate A matrix:
// Create x vector

// Create b vector

// Solve for each timestep:
for (each timestep)
    for (X,Y,Z)

        KSPSetOperators(ksp, Amat, Amat, DIFFERENT_NONZERO_PATTERN);


1 Answer 1

  1. ADI is not a very good parallel algorithm. You should seriously consider formulating the problem in 3D and solving with multigrid. You could get a start with src/ksp/ksp/examples/tutorials/ex45.c.

  2. If you insist on using ADI, you should seriously consider allocating the three matrices separately. It's more memory, but then you won't have to reassemble on each iteration, so you amortize setup costs. This uses more memory.

  3. If you insist on ADI and on reusing the same matrix for each direction, you should allocate for three nonzeros per row and fill in explicit zeros even when the direction is missing an entry. Otherwise the first assembly compacts out the unused entries and they can't be cheaply inserted later.

  • $\begingroup$ I'm actually building several different solution techniques, ADI being one of them, so that I can compare relative performance. This being the case, I take it that you answers are in order of recommendation? If so, how is answer 2 a better solution than answer 3? I have to reset the values in the matrices (or at last some of them) due to changing boundary conditions anyway. $\endgroup$
    – Neal Kruis
    Commented Nov 25, 2013 at 22:46
  • $\begingroup$ How often do the boundary conditions change? What kind of boundary conditions and how do they change? (It may be possible to formulate so the matrix entries do not change.) Do you converge the ADI or only do one step? $\endgroup$
    – Jed Brown
    Commented Nov 26, 2013 at 4:27
  • $\begingroup$ I have Neumann type boundary conditions where the derivative at the boundary is dependent on the temperature at the boundary and a time varying constant (i.e. thermal convection: dT/dx = f(T)*(T - C)). I don't know what you mean by converging the ADI, but the method I use is based on this reference: tandfonline.com/doi/abs/10.1080/… $\endgroup$
    – Neal Kruis
    Commented Nov 26, 2013 at 17:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.