New linear system solvers

chegdan · April 28, 2009, 00:22

Hi All,

Suppose that I want to rewrite a linear system solver in another programming language and implement it in OpenFOAM. If I begin with PBiCG.C, how can I read all of the values from the "A" matrix and "b" vector so they can be moved to a new list in the other programming language? In other words, what type of sparse matrix storage is used in OpenFOAM? Are we talking "Modified Sparse Row" or "Compressed Spares Column" type or something else? How do lduAddressing, matrix_, coupleBouCoeffs_, coupleIntCoeffs_, and interfaces_ relate to each other? I have looked into PBiCG.C and it would be great to have some help. Suggestions?

Dan

gschaider · April 28, 2009, 13:31

Quote:

Originally Posted by chegdan

Hi All,

Suppose that I want to rewrite a linear system solver in another programming language and implement it in OpenFOAM. If I begin with PBiCG.C, how can I read all of the values from the "A" matrix and "b" vector so they can be moved to a new list in the other programming language? In other words, what type of sparse matrix storage is used in OpenFOAM? Are we talking "Modified Sparse Row" or "Compressed Spares Column" type or something else? How do lduAddressing, matrix_, coupleBouCoeffs_, coupleIntCoeffs_, and interfaces_ relate to each other? I have looked into PBiCG.C and it would be great to have some help. Suggestions?

Dan

I never fiddled around with the linear solvers, so I can't give you much help. But there seems to be some information on addressing in the header-file of lduAddressing. Other than that I can only recommend clicking through the Doxygen.

Bernhard

joern · April 28, 2009, 14:55

Uff, i analysed the linear Solvers last year. I try to explain, how the linear systems are organised.

First, the class fvMatrix is not a Matrix, it is a whole linear system. The class lduMatrix is a whole system, too.

I think the interesting part for you are the arrays:
upper(), lower(), diag(), source()

diag and source should be clear. upper and lower contain the U and L parts of the lduMatrix. the elements of these arrays get adressed by upperAddr(), lowerAddr() and ownerStartAddr().
upperAddr contains the rows of Upper and the colums of lower, lowerAddr the colums of upper and the rows of lower.

So the arrays are sortet for A*x multiplication.
If you want to get Transpose[A] you have to use
upper()[losortAddr()], lower()[losortAddr()], upperAddr()[losortAddr()], lowerAddr()[losortAddr()] and losortStartAddr() for adressing.

all the arrays, i mentioned, are acessible over member-functions of the lduMatrix.

hope that helps.

chegdan · May 7, 2009, 13:48

Thanks for your replies.
Bernhard: I did see the header file and it was a little bit confusing, but I think I'm getting it. Thanks for your pointers.

Joern: I will look more into what you've posted. I recently have been handed another short-term project and have to leave my efforts for a week or two more. Thanks for your help though, its greatly appreciated.

Dan

henrik · May 8, 2009, 04:09

Dear Dan,

I am not sure what you are planning to do, but it might interesting for you to know that there is a set of enhanced AMG solvers in OpenFOAM-extend. More information can be found in this paper.

http://powerlab.fsb.hr/ped/kturbo/Op...pmAIAA2007.pdf

Best Regards,

Henrik

chegdan · May 8, 2009, 12:53

Henrik,

I'll take a look at the paper. I have already downloaded the -dev version and compiled successfully. However, I have not looked at the new AMG solvers. Thanks for pointing that out. I'm actually looking at the possibility of using CUDA as a linear system solver and needed to understand how the matrices are addressed so they can be passed to the GPU. This is not my main project, only a little fun. I have seen several threads pertaining to porting openFoam to run on a GPU...but I think its a bit overkill to port more than the linear system solvers. Thanks again.

Dan

sorcerer · July 22, 2009, 08:14

Quote:

Originally Posted by chegdan

This is not my main project, only a little fun. I have seen several threads pertaining to porting openFoam to run on a GPU...but I think its a bit overkill to port more than the linear system solvers. Thanks again.

Hi Dan,

I think you are right with only porting the linear system solvers.
How is the state of your project? Would you share your experience?

Greetings,
Heiko

chegdan · August 25, 2009, 16:25

I had been working on developing in CUDA for a bit, but took some time off since another project that took priority was put in front of me. I haven't finished my solver yet, but now I'm looking into OpenCL as an alternative to CUDA. OpenCL might be better option for those not using nvidia products. I will look into it in the next few months once I have more time.

Dan

NickG · March 23, 2010, 08:13

Hi Dan

How's the port going? Any luck?

I may be a bit after the event but...

Looking at the thesis http://ce.et.tudelft.nl/publicationf...063_thesis.pdf it would seem that compressed row or modified sparse row formats are better for parallelisation than the ldu format, although I wonder if some extra temporary local variables could overcome this problem. The compressed row format already has a converter written in openFoam so unless you want a lot of rewriting it may be easiest to use that, although you will lose some of the time that you'd gain in CUDA with the conversion process.

If you want to know how all this hangs together then you'll have to ask someone else as I've only just started down this path...

Cheers

Nick

johndeas · April 30, 2010, 05:41

Hi,

would it be possible to have a test case that uses RPM(3) as a solver ?

Thanks !

NickG · April 30, 2010, 08:05

Hi John

What's RPM(3)?

johndeas · April 30, 2010, 10:22

3rd Recursive Projection Method, this is what the article cited by Prof. Hrvoje Jasak is about.

Article to be found at http://powerlab.fsb.hr/ped/kturbo/Op...pmAIAA2007.pdf

April 28, 2009, 00:22	New linear system solvers	#1
chegdan Senior Member Daniel P. Combest Join Date: Mar 2009 Location: St. Louis, USA Posts: 621 Rep Power: 0	Hi All, Suppose that I want to rewrite a linear system solver in another programming language and implement it in OpenFOAM. If I begin with PBiCG.C, how can I read all of the values from the "A" matrix and "b" vector so they can be moved to a new list in the other programming language? In other words, what type of sparse matrix storage is used in OpenFOAM? Are we talking "Modified Sparse Row" or "Compressed Spares Column" type or something else? How do lduAddressing, matrix_, coupleBouCoeffs_, coupleIntCoeffs_, and interfaces_ relate to each other? I have looked into PBiCG.C and it would be great to have some help. Suggestions? Dan

April 28, 2009, 14:55		#3
joern Member Joern Bader Join Date: Mar 2009 Posts: 33 Rep Power: 17	Uff, i analysed the linear Solvers last year. I try to explain, how the linear systems are organised. First, the class fvMatrix is not a Matrix, it is a whole linear system. The class lduMatrix is a whole system, too. I think the interesting part for you are the arrays: upper(), lower(), diag(), source() diag and source should be clear. upper and lower contain the U and L parts of the lduMatrix. the elements of these arrays get adressed by upperAddr(), lowerAddr() and ownerStartAddr(). upperAddr contains the rows of Upper and the colums of lower, lowerAddr the colums of upper and the rows of lower. So the arrays are sortet for A*x multiplication. If you want to get Transpose[A] you have to use upper()[losortAddr()], lower()[losortAddr()], upperAddr()[losortAddr()], lowerAddr()[losortAddr()] and losortStartAddr() for adressing. all the arrays, i mentioned, are acessible over member-functions of the lduMatrix. hope that helps. solefire, songwukong, hua1015 and 1 others like this.

May 8, 2009, 12:53	I'll take a look, thanks .	#6
chegdan Senior Member Daniel P. Combest Join Date: Mar 2009 Location: St. Louis, USA Posts: 621 Rep Power: 0	Henrik, I'll take a look at the paper. I have already downloaded the -dev version and compiled successfully. However, I have not looked at the new AMG solvers. Thanks for pointing that out. I'm actually looking at the possibility of using CUDA as a linear system solver and needed to understand how the matrices are addressed so they can be passed to the GPU. This is not my main project, only a little fun. I have seen several threads pertaining to porting openFoam to run on a GPU...but I think its a bit overkill to port more than the linear system solvers. Thanks again. Dan

August 25, 2009, 16:25	OpenCL may or may not be better	#8
chegdan Senior Member Daniel P. Combest Join Date: Mar 2009 Location: St. Louis, USA Posts: 621 Rep Power: 0	I had been working on developing in CUDA for a bit, but took some time off since another project that took priority was put in front of me. I haven't finished my solver yet, but now I'm looking into OpenCL as an alternative to CUDA. OpenCL might be better option for those not using nvidia products. I will look into it in the next few months once I have more time. Dan

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Direct vs Iterative Linear Solvers for non-linears	bill	Main CFD Forum	16	November 5, 2014 07:18
CFX11 + Fortran compiler ?	Mohan	CFX	20	March 30, 2011 18:56
Linear system solver	elias	FLUENT	0	November 19, 2008 09:30
flop count for linear system solvers	no_name	Main CFD Forum	2	December 7, 2005 22:10
More advise on libraries to solve linear system	Ricardo Bonon	Main CFD Forum	0	May 18, 2000 18:37

May 7, 2009, 13:48		#4
chegdan Senior Member Daniel P. Combest Join Date: Mar 2009 Location: St. Louis, USA Posts: 621 Rep Power: 0	Thanks for your replies. Bernhard: I did see the header file and it was a little bit confusing, but I think I'm getting it. Thanks for your pointers. Joern: I will look more into what you've posted. I recently have been handed another short-term project and have to leave my efforts for a week or two more. Thanks for your help though, its greatly appreciated. Dan

May 8, 2009, 04:09		#5
henrik Senior Member Henrik Rusche Join Date: Mar 2009 Location: Wernigerode, Sachsen-Anhalt, Germany Posts: 281 Rep Power: 18	Dear Dan, I am not sure what you are planning to do, but it might interesting for you to know that there is a set of enhanced AMG solvers in OpenFOAM-extend. More information can be found in this paper. http://powerlab.fsb.hr/ped/kturbo/Op...pmAIAA2007.pdf Best Regards, Henrik

March 23, 2010, 08:13		#9
NickG Member Nick Gardiner Join Date: Apr 2009 Location: Chichester, UK Posts: 94 Rep Power: 17	Hi Dan How's the port going? Any luck? I may be a bit after the event but... Looking at the thesis http://ce.et.tudelft.nl/publicationf...063_thesis.pdf it would seem that compressed row or modified sparse row formats are better for parallelisation than the ldu format, although I wonder if some extra temporary local variables could overcome this problem. The compressed row format already has a converter written in openFoam so unless you want a lot of rewriting it may be easiest to use that, although you will lose some of the time that you'd gain in CUDA with the conversion process. If you want to know how all this hangs together then you'll have to ask someone else as I've only just started down this path... Cheers Nick

April 30, 2010, 05:41		#10
johndeas Senior Member John Deas Join Date: Mar 2009 Posts: 160 Rep Power: 17	Hi, would it be possible to have a test case that uses RPM(3) as a solver ? Thanks !

April 30, 2010, 08:05		#11
NickG Member Nick Gardiner Join Date: Apr 2009 Location: Chichester, UK Posts: 94 Rep Power: 17	Hi John What's RPM(3)?

April 30, 2010, 10:22		#12
johndeas Senior Member John Deas Join Date: Mar 2009 Posts: 160 Rep Power: 17	3rd Recursive Projection Method, this is what the article cited by Prof. Hrvoje Jasak is about. Article to be found at http://powerlab.fsb.hr/ped/kturbo/Op...pmAIAA2007.pdf