10
$\begingroup$

I am solving Turing's reaction-diffusion system with following C++ code. It is too slow: for 128x128 pixel texture, acceptable number of iterations is 200 – which results in 2.5 seconds of delay. I need 400 iterations to obtain interesting image – but 5 seconds of waiting is too much. Also, size of the texture should be in fact 512x512 – but this results in huge waiting time. The devices are iPad, iPod.

Is there any chance to do this faster? Euler method converges slowly (wikipedia) – having quicker method would allow to drop number of iterations?

EDIT: As Thomas Klimpel pointed out, the lines: "if( m_An[i][j] < 0.0 ) { ... }", "if( m_Bn[i][j] < 0.0 ) { ... }" are delaying convergence: after removing, meaningful image appears after 75 iterations. I have commented out the lines in code below.

void TuringSystem::solve( int iterations, double CA, double CB ) {
    m_iterations = iterations;
    m_CA = CA;
    m_CB = CB;

    solveProcess();
}

void set_torus( int & x_plus1, int & x_minus1, int x, int size ) {
    // Wrap "edges"
    x_plus1 = x+1;
    x_minus1 = x-1;
    if( x == size - 1 ) { x_plus1 = 0; }
    if( x == 0 ) { x_minus1 = size - 1; }
}

void TuringSystem::solveProcess() {
    int n, i, j, i_add1, i_sub1, j_add1, j_sub1;
    double DiA, ReA, DiB, ReB;

    // uses Euler's method to solve the diff eqns
    for( n=0; n < m_iterations; ++n ) {
        for( i=0; i < m_height; ++i ) {
            set_torus(i_add1, i_sub1, i, m_height);

            for( j=0; j < m_width; ++j ) {
                set_torus(j_add1, j_sub1, j, m_width);

                // Component A
                DiA = m_CA * ( m_Ao[i_add1][j] - 2.0 * m_Ao[i][j] + m_Ao[i_sub1][j]   +   m_Ao[i][j_add1] - 2.0 * m_Ao[i][j] + m_Ao[i][j_sub1] );
                ReA = m_Ao[i][j] * m_Bo[i][j] - m_Ao[i][j] - 12.0;
                m_An[i][j] = m_Ao[i][j] + 0.01 * (ReA + DiA);
                // if( m_An[i][j] < 0.0 ) { m_An[i][j] = 0.0; }

                // Component B
                DiB = m_CB * ( m_Bo[i_add1][j] - 2.0 * m_Bo[i][j] + m_Bo[i_sub1][j]   +   m_Bo[i][j_add1] - 2.0 * m_Bo[i][j] + m_Bo[i][j_sub1] );
                ReB = 16.0 - m_Ao[i][j] * m_Bo[i][j];
                m_Bn[i][j] = m_Bo[i][j] + 0.01 * (ReB + DiB);
                // if( m_Bn[i][j] < 0.0 ) { m_Bn[i][j]=0.0; }
            }
        }

        // Swap Ao for An, Bo for Bn
        swapBuffers();
    }
}
|cite|improve this question
$\endgroup$
  • $\begingroup$ Also, I want to mention that it is preferred that you don't cross-post questions, since it appears you have asked very similar questions both here and here. $\endgroup$ – Godric Seer Jul 22 '12 at 0:16
  • $\begingroup$ Have you already seen Greg Turk's work on this, by any chance? $\endgroup$ – J. M. Jul 22 '12 at 3:44
  • $\begingroup$ @J.M.: Not yet. I just tried run his code: it requires X server with PseudoColor, i.e. 8 bit color depth. I think I cannot provide this on OSX. I tried various VNC servers but no luck. $\endgroup$ – AllCoder Jul 22 '12 at 15:30
  • $\begingroup$ I think you should still be able to adapt Turk's approach to the matter at hand; reaction-diffusion patterns seem to be used a fair bit in computer graphics nowadays. $\endgroup$ – J. M. Jul 22 '12 at 15:32
  • 1
    $\begingroup$ I might be wrong, but the part with m_An[i][j] = 0.0; might actually add an element to this system that cannot be modeled by a differential equation with a continuous right hand side. This makes it a bit difficult to come up with a faster solver. $\endgroup$ – Thomas Klimpel Jul 22 '12 at 21:56
9
$\begingroup$

You seem to be limited by stability, which is expected since diffusion is stiff as you refine the grid. Good methods for stiff systems are at least partly implicit. It will take some effort, but you can implement a simple multigrid algorithm (or use a library) to solve this system with a cost of less than ten "work units" (essentially the cost of one of your time steps). When you refine the grid, the number of iterations will not increase.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ If only diffusion would be stiff here, he could use an ADI method like Douglas-Gunn and everything would be fine. However, in my own experience, the reaction part is often much worse with respect to stiffness in addition to being badly nonlinear. $\endgroup$ – Thomas Klimpel Jul 22 '12 at 14:27
  • 1
    $\begingroup$ ADI unfortunately has terrible memory locality. Also note that reaction can be treated implicitly regardless of whether diffusion is. Under grid refinement, diffusion will eventually become dominant, but we can't tell where the threshold is without knowing the constants. $\endgroup$ – Jed Brown Jul 22 '12 at 15:05
  • $\begingroup$ Example code implementing backward Euler for this (in Python) is here: scicomp.stackexchange.com/a/2247/123 $\endgroup$ – David Ketcheson Jul 22 '12 at 20:08
  • $\begingroup$ @DavidKetcheson: Using implicit methods requires solving an equation? This is why there is linalg.spsolve() in the code? $\endgroup$ – AllCoder Jul 22 '12 at 22:00
  • 1
    $\begingroup$ @AllCoder Yes, it requires a solve, but the solve can be done much faster than all time steps required for an explicit method to be stable. $\endgroup$ – Jed Brown Jul 22 '12 at 22:05
2
$\begingroup$

From a practical point of view: the A5 processor is not that much powerful, so you can wait a few HW iterations, or if your ipod / ipad are going to be connected to the internet, solve your problem remotely or in the cloud.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ I am surprised how small power the A5 offers. How can Pages, Safari and other large applications work so well? I need to generate random, abstract images, thought that the morphogenesis will be enough simple.. $\endgroup$ – AllCoder Jul 22 '12 at 14:57
  • $\begingroup$ Well, A5 is a energy-efficient processor optimised for web and video (Pages, Safari, etc). In contrast, most numerical workloads perform tons of floating-point operations and data-movements, these features are not the focus of a low-power mobile processors. $\endgroup$ – fcruz Jul 22 '12 at 16:25
0
$\begingroup$

Euler does converge slowly in relation to other methods, however I don't think that is what you are interested. If you are just looking for "interesting" images, increase the size of your time step and take less iterations. The problem, as Jed points out, is that the explicit euler method has stability issues with large time steps in relation to grid size. the smaller your grid is (i.e. the higher resolution your image is), the smaller your time step must be to account for it.

For instance, by using implicit euler instead of explicit, you don't gain any orders of convergence, but the solution will have unconditional stability, allowing much larger time steps. Implicit methods are more complicated to implement and take more computation per time-step, but you should see gains well beyond that by taking less steps in total.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ This problem is limited by stability so simply increasing the time step size will not work. $\endgroup$ – Jed Brown Jul 22 '12 at 0:49
  • $\begingroup$ If I change 0.01 to e.g. 0.015, then I get "concentration of chem. sp. near zero" at all points – i.e. a gray square. Here is origin of my code: drdobbs.com/article/print?articleId=184410024 $\endgroup$ – AllCoder Jul 22 '12 at 0:50
  • $\begingroup$ Yes, that would be a result of the stability issues Jed mentioned. As he mentions in his answer, using an implicit method characterized by better stability performance will resolve this issue for you. I will update my answer to remove the irrelevant information. $\endgroup$ – Godric Seer Jul 22 '12 at 0:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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