Re: [Boost-users] multi-array and pdes

Subject: Re: [Boost-users] multi-array and pdes
From: Larry Evans (cppljevans_at_[hidden])
Date: 2011-10-27 14:13:28

On 05/24/11 00:58, pmamales_at_[hidden] wrote:
> Hi Larry,
> 
> Thank you very much for your extended response.
> I am not sure, will have to think about it, altough this seems right.
> In appreciation of your effort to help me, let me give you some color:
> Say I am trying to sove a 3d problem using splitting methods. Lets say that the original
> system f reference is xyz.
> One alays ends up to a system of equations in the vectorized reprezentation of the grid (very much like the
> array where the elements of the ma are stored). 
> Then, when trying to solve the problem in the x direction (while in fortran storage scheme),
> I obtain a nice tridiagonal system of equations which I can solve very efficiently (using Thomas algorithm which is O(N) ).
> When I go to the second dimension, the tridiagonal system is hidden (in the original vector). However, in the rotsted yzx system it is there!!
> So, I am "paying the price" of writting it in the rotated view and copying it into a different vector to solve efficiently
> (o/wise I end up with a much more expensive sparse system to solve..)
> Hence the root of my questions.
> One of the issues is that the rows that I need to permute is rather complicated and with a little of smoke and mirrors, I hope to bypass 
> the complications.
> Apologies, if all this is known to ou already or if went way too atray. Also, apologies for not catching your response earlier.
> Thank you agian for all your help-will look inth the references,
> Petros
[snip]
Hi Petros,

I've just uploaded:

http://svn.boost.org/svn/boost/sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/array_dyn.diff_pde.cpp

http://svn.boost.org/svn/boost/sandbox/variadic_templates/sandbox/stepper/libs/array_stepper/examples/solve_tridiag.cpp

The *pde.cpp demonstrates 2 methods for solving pde's:

  1) Alternating Directions Explicit (ADE)

  2) Alternating Directions Implicit (ADI)

References documenting the methods and many in source
comments hopefully make it fairly easy to see how the
code corresponds to the reference.  If not, please let
me know.

The ADI method uses the solve_tridiag and rotated axes.
The rotation does not rotate the axes of the array, it
simply changes 1 scalar member variable in the:

  index_stack_length_stride_crtp_indices<>

class.  I hope that's what you meant by rotating axes.
There are TRACE_* macros in the solve_tridiag.hpp file which
illustrate what axis is actually being solve for.

One thing that worries me is that the relative error for
ADI is about 10x what the relative error for ADE.  I've
not yet figured out why :(

HTH.

-regards,
Larry