# coupleBouCoeffs, coupleIntCoeffs, and interfaces oh my

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

 July 1, 2011, 19:10 coupleBouCoeffs, coupleIntCoeffs, and interfaces oh my #1 Senior Member     Daniel P. Combest Join Date: Mar 2009 Location: St. Louis, USA Posts: 538 Rep Power: 18 Hello All, I am trying to understand the parallel (bi)conjugate gradient solvers in OF and have a good grasp of the lduMatrix format. I have programmed some serial solvers myself, moved the contents of lduMatrix to my solver, and returned the result to OF. I would really like to make them work in parallel in OF. I know that OF uses the domain decomposition method to parallelize, so I would like to pull out the interface (coupling) matrices. I'm assuming that they are in the coupleBouCoeffs and coupleIntCoeffs objects and I know they are used in the Amul methods of the lduMatrix. Has anyone been able to extract these interface coefficients (coupling matrix)? Any help would be much appreciated. Dan

 July 2, 2011, 22:19 #2 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 397 Rep Power: 12 Hi Dan, I've been trying to figure this out myself. Here's what I've got so far: OpenFOAM's lduAddressing falls apart when out-of-core multiplication is required (e.g. when off-diagonal terms arise from boundary conditions). The solution is the coupling matrix, and it has to accommodate: - processor - cyclic - ggi - etc. The cyclic and ggi patch mean that a coordinate transformation may be required across the interface. The processor patch means that the associated coupling is between two different lduMatrices. I'm fuzzy on this bit, but theoretically, the effect of the coupling matrix is to insert extra columns in the main lduMatrix so that the off-diagonal terms are actually on the diagonal. For instance, if cell 5 couples with cell 72, we create an extra column 5' just after 72. Practically, I'm not sure how this is achieved. As you are probably aware, there are three components to the coupling matrix: -an "lduInterfaceFieldPtrsList" object; -internalCoeffs; and -boundaryCoeffs. I believe the lduInterfaceFieldPtrsList's job is to: -know which cells are coupled -know which processors they reside on -perform coordinate transformations The internal and boundary coeffs are associated with the 'diagonal' and source term requirements of the boundary... but I'm not entirely sure how. Lastly, to access them, use: pEqn.boundaryCoeffs()[patchNumber][faceNumber] pEqn.internalCoeffs()[patchNumber][faceNumber] Hope you can fill in the blanks... -Dave

July 21, 2011, 16:49
#3
Senior Member

Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 538
Rep Power: 18
There were a few offline discussions that I wanted to put up here for others:

Dave,
Fantastic! I will include some of the things that I've found out as well. For instance, there are several presentations by Dr. Jasak that discuss this topic very briefly (say essentially the same thing):
[1] http://wtec2.com/ScalableSoftware/do...sak-slides.pdf
[2] http://www.prace-ri.eu/IMG/pdf/Wikki.pdf
From these presentations, some important info that confirms some stuff we already know (or can find via the message board):
* In terms of code organization, each sub-domain creates its own numbering space: locally, equation numbering always starts with zero and one cannot rely on global numbering: it breaks parallel efficiency
* Processor interfaces are updated separately, involving communications
* Uses a zero halo approach, i.e. non-overlapping decomposed domains and communication concerning boundary values on patches.
Also, the domain decomposition method in this case (the case of parallelization of the linear system solvers) is not Neumann-Neumann or Dirichlet-Dirichlet (DD) etc. since those are sequential methods to solve coupled domains. For example, the DD (solves Dirichlet bc on one domain, transfers info to second domain, transmits result back to first domain via Neumann bc, then process repeated until no change between domains) method is however used in some of the conjugate solvers (cht*Foam? solvers ) that do not use the block matrix methods like those used in conjugateHeatFoam.

The way I have understood what is happens in parallelization of conjugate gradient (from the book http://www.amazon.com/Computational-...9732671&sr=8-1 ) via domain decomposition is that we first split the global domain into some parts with processor interfaces (figure 1). each of the subdomains have an lduMatrix that is stored and associated with each processor. Lets take 6 processors and have 6 domains. Now we are left with a global Ax=b system that is solved on 6 subdomains that interact with eachother. Imagine again that we assemble a global matrix (figure 2) where each of the lduMatrices for each subdomain are stored in Aii. Each of the coupling domains are stored in Aij (with j!= I) and a coupling matrix is only stored for each domain that interfaces another (i.e. A12 is stored since domains 1 and 2 interact with eachother. However A16 is not necessary since they do not exchange information to eachother.) Suppose we want to get the residual vector (r1) for the first processor, we then need to do r1 = b1 - (A11*x1 + A12*x2+A14*x4). Processor 1 will need the x vectors from processor 2 and 4 since they are coupled to processor 1. Figure 3 shows the information that processor 1 will need to calculate the residual vector for processor 1. the key in getting another solver to work is grabbing these Aij values and doing the bookkeeping. Hope this helps...my arm hurts from typing this one today.
Dan
Attached Files
 parallelIdea.pdf (31.5 KB, 69 views)

 July 21, 2011, 17:09 #4 Senior Member     Daniel P. Combest Join Date: Mar 2009 Location: St. Louis, USA Posts: 538 Rep Power: 18 Now in a little more follow up on the topic. I was able to figure out how to get the interfaces to identify themselves as processor interfaces using this code: Code: ``` register label nInterfaces = interfaces().size(); Info<<"the number of interfaces "<

October 4, 2011, 08:07
#5
New Member

Roland Engberg
Join Date: Jan 2011
Posts: 14
Rep Power: 6
Quote:
 Originally Posted by chegdan [...] However, I would like to call a method specific to a processorLduInterfaceField, especially the myProcNo() and neighbProcNo() methods ( I know that I can get myProcNo() from Pstream::myProcNo()...but i need the neighboring processor that the interface is communicating to and from). [...]
Hello Dan,

thank you for your helpful post! Did you figure out how to access myProcNo() and neighbProcNo()? It is in the boundary file so there has to be a command to retrieve this information. When I try something like the following, I get a compilation error:

Code:
```const fvPatchList& patches = psi.mesh().boundary();

forAll(patches, patchI)
{
const fvPatch& p = patches[patchI];
const polyPatch& pp = p.patch();

if (typeid(pp) == typeid(processorPolyPatch))
{
Pout << pp.neighbProcNo();

[...]```
Best wishes,
Roland

 October 4, 2011, 21:47 #6 Senior Member     Daniel P. Combest Join Date: Mar 2009 Location: St. Louis, USA Posts: 538 Rep Power: 18 Roland, I ended up changing an OpenFOAM core file in the lduInterface.H by adding the virtual functions: Code: ``` //- Return the processor number of the interface //dcombest July 26 2011 virtual int myProcNo() const { return (-1); } //- Return neigbour processor number of the interface //dcombest July 26 2011 virtual int neighbProcNo() const { return (-1); }``` I had it return -1 since it was a little bit faster than querying if the interface was a "processor" interface i.e. string compare vs. integer compare. The virtual functions will return the correct value of the myProc and naighborProc if it is in parallel. Hope this helps. Dan __________________ Dan Find me on twitter @dancombest and LinkedIn

 Tags interfaces, ldumatrix, linear solver

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

All times are GMT -4. The time now is 23:37.