CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   Problems with biCGSTAB (http://www.cfd-online.com/Forums/main/88274-problems-bicgstab.html)

Petrov May 12, 2011 11:03

Problems with biCGSTAB
 
I wrote a biconjugate gradient method in C for my FEM solver, but it does not work. (I can tell because when I use a different method, the results are correct). Can you see if there is an obvious problem with the code below?
Code:

        r = addVectors(b, multiplyR(A, x0), n, -1);
        copyVector(r, rtilde, n);
        int i;
        for(i=0;i<n;i++)
        {
                rho_1 = dotProduct(rtilde, r, n);
                if(i==0)
                {
                        copyVector(r, p, n);
                } else {
                        beta = (rho_1/rho_2) * (alpha/omega);
                        double * tempA = addVectors(p, v, n, -omega);
                        p = addVectors(r, tempA, n, beta);
                        free(tempA);
                }
                v = multiplyR(A, p);
                alpha = rho_1 / (dotProduct(rtilde, v, n));
                s = addVectors(r, v, n, (-1)*alpha);
                t = multiplyR(A, s);
                omega = (dotProduct(t, s, n))/dotProduct(t, t, n);
                double* temp = addVectors(x, p, n, alpha);
                x = addVectors(temp, s, n, omega);
                free(temp);
                r = addVectors(s, t, n, -omega);
                rho_2 = rho_1;
        }

The code you see is the "meat" of the solver. Before comes the allocation, and after comes the free-ing of memory. The input of the function is A, the matrix; b the vector, and x0 the initial guess.

multiplyR(matrix, vector) is a function that multiplies a matrix on the left with a vector on the right.

addVectors(vector A, vector B, int n, double k) returns A+kB. (The n is length of each vector).

As you can see, I'm not using a pre-conditioner.

Any help would be much appreciated!!

:::UPDATE:::
I have recently tried using a gauss-siedel preconditioner. The results did not improve at all, but they did change. Does that mean that I should try more sophisticated preconditioners, that there something wrong with the algorithm above, or that I should move on to gmres or some other method?

(I think I'm just going to try gaussian elimination until I feel like tackling this problem again. I've spent far too long on it!!!)


All times are GMT -4. The time now is 09:08.