CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

LUDecompose, LUBacksubstitue in class scalarMatrix

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By jaswi

Reply
 
LinkBack Thread Tools Display Modes
Old   July 1, 2009, 04:43
Default LUDecompose, LUBacksubstitue in class scalarMatrix
  #1
Senior Member
 
Join Date: Mar 2009
Posts: 248
Rep Power: 9
jaswi is on a distinguished road
Dear Developers

Good Day to All !

While using the scalarMatrix class (/src/OpenFOAM/matrices/scalarMatrix) I tried
the LUDecompose and LUBacksubstitute functions. They work fine but when i tried to compare the results with the MATLAB answers , it seems the decomposed matrix has some special format to store the L and U parts. Here is what I did:


Code:
int main(int argc, char *argv[])
{
    scalarMatrix hmm(3);

    hmm[0][0] = -3.0;
    hmm[0][1] = 10.0;
    hmm[0][2] = -4.0;
    hmm[1][0] = 2.0;
    hmm[1][1] = 3.0;
    hmm[1][2] = 10.0;
    hmm[2][0] = 2.0;
    hmm[2][1] = 6.0;
    hmm[2][2] = 1.0;

    Info<< hmm << endl;

    labelList pivotIndices(hmm.n());

    hmm.LUDecompose(hmm, pivotIndices);

    scalarField source(hmm.n(), 0);

    source[0] = 1;
    source[1] = 3;
    source[2] = 4;

    hmm.LUBacksubstitute(hmm, pivotIndices, source);

    Info << hmm << endl;

    Info << source << endl;

    Info << "End\n" << endl;

    return 0;
}
the matrix is a 3 x 3 square matrix

hmm = ((-3 10 -4)
(2 3 10)
(2 6 1))

According to the class documentation LUDecompose will decompose the matrix with pivoting, which in my understanding implies MATLAB command [L,U] = lu(A).
I called hmm.LUDecompose(hmm, pivotIndices), and assuming that the decomposition is saved back in the scalarMatrix hmm , the Info << hmm << endl;
after decomposition gives:
3 3 ((2 6 1)
(-1.5 19 -2.5)
(1 -0.157895 8.60526))

and Info << source << endl; gives the solution

3(0.883792 0.370031 0.0122324)

Now on the other hand using MATLAB one gets:
Code:
Lower Matrix L:
1            0             0
-0.6667   1             0
-0.6667   1.31034   1

Upper Matrix U:

-3      10         -4
0       9.667     7.333
0       0          -11.275

Solution :
0.8837    0.370036   0.0122324
The solution matches perfectly in both cases but the matrix returned by OpenFOAM has special storage format which is not clear to me. I hope Gurus might help me to understand this . I tried some additions and subtractions but could not find the storage pattern.

Thanks for attention
BR
jaswi
Tushar@cfd likes this.
jaswi is offline   Reply With Quote

Old   July 2, 2009, 02:52
Default
  #2
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Jaswinder,

just a quick pointer: Have a look into the ODE solver classes and how scalarMatrix is used there.

Henrik
henrik is offline   Reply With Quote

Old   July 2, 2009, 05:53
Default
  #3
Senior Member
 
Join Date: Mar 2009
Posts: 248
Rep Power: 9
jaswi is on a distinguished road
Hi Henrik

Thanks for the tip. I will follow that clue and post back .

Best Regards
jaswi
jaswi is offline   Reply With Quote

Old   July 2, 2009, 06:52
Default
  #4
Super Moderator
 
niklas's Avatar
 
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 19
niklas will become famous soon enough
a third variant
i tried it in octave for fun...
[l,u,p] = lu(hmm)

and I got this

l =

1.00000 0.00000 0.00000
-0.66667 1.00000 0.00000
-0.66667 0.76316 1.00000

u =

-3.00000 10.00000 -4.00000
0.00000 12.66667 -1.66667
0.00000 0.00000 8.60526

p =

1 0 0
0 0 1
0 1 0

This is what 'help lu' says:
Compute the LU decomposition of A, using subroutines from LAPACK.
niklas is offline   Reply With Quote

Old   July 8, 2009, 03:37
Default
  #5
Senior Member
 
Henrik Rusche
Join Date: Mar 2009
Location: Braunschweig, Niedersachsen, Germany
Posts: 275
Rep Power: 9
henrik is on a distinguished road
Dear Jaswi,

you need to understand what pivotIndices is for, which takes the value (2 2 2) in my case.

Have a look into Numerical Recipes (Press et al.).

Henrik
henrik is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
STAR-CCM+ Java class Claudene CD-adapco 1 July 17, 2008 05:49
Errors running allwmake in OpenFOAM141dev with WM_COMPILE_OPTION%3ddebug unoder OpenFOAM Installation 11 January 30, 2008 21:30
How to add a new class locally ville OpenFOAM 4 December 11, 2006 14:20
About UList and List class leosding OpenFOAM Running, Solving & CFD 1 December 2, 2005 00:52
Expanding a class fabianpk OpenFOAM 0 October 3, 2005 04:26


All times are GMT -4. The time now is 15:19.