
[Sponsors] 
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: 613
Rep Power: 22 
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: 436
Rep Power: 15 
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 outofcore multiplication is required (e.g. when offdiagonal 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 offdiagonal 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: 613
Rep Power: 22 
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...sakslides.pdf [2] http://www.praceri.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 subdomain 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. nonoverlapping 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 NeumannNeumann or DirichletDirichlet (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=81 ) 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 

July 21, 2011, 17:09 

#4 
Senior Member
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 613
Rep Power: 22 
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 "<<nInterfaces<<"\n"; word thisInterfaceType; forAll (interfaces(), interfaceI) { if(interfaces().set(interfaceI)){ //get the type of the interface thisInterfaceType = interfaces()[interfaceI].interfaceFieldType(); if(thisInterfaceType == "processor"){ Info<<"Interface "<<interfaceI<<" of "<<nInterfaces<<" is a "<<thisInterfaceType<<" boundary\n"; } } } Lastly, I can grab the coefficients from the coupleBouCoeffs and coupleIntCoeffs objects in a similar way that lduMatrixATmul.C did in the sumA method, but... 1) when referring to face...is this the owner or neighbor (i.e. upper or lower) values that I want to use in a coupling matrix (see previous post attachment)? 2) basically, how is the addressing working for the sumA method similar to Code:
const unallocLabelList& pa = lduAddr().patchAddr(patchI); const scalarField& pCoeffs = coupleBouCoeffs[patchI]; Code:
forAll(pa, face) { sumAPtr[pa[face]] = pCoeffs[face]; } Dan 

October 4, 2011, 08:07 

#5  
New Member
Roland Engberg
Join Date: Jan 2011
Posts: 14
Rep Power: 8 
Quote:
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(); [...] Roland 

October 4, 2011, 21:47 

#6 
Senior Member
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 613
Rep Power: 22 
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); } Dan 

Tags 
interfaces, ldumatrix, linear solver 
Thread Tools  
Display Modes  

