[chegdan]

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

[marupio]

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

[chegdan]

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

[chegdan]

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:

 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";

		}
	}
    }

This is dropped into the solver on the level of the PCG solver in the lduMatrix solvers. When run in parallel, it should report which interfaces are processor interfaces, but only report it once since we are using "Info". If we used "cout" then it would report it from every processor. I can find the interfaces that are processors from the interface list of type lduInterfaceFieldPtrsList&. 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).

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

            const unallocLabelList& pa = lduAddr().patchAddr(patchI);
            const scalarField& pCoeffs = coupleBouCoeffs[patchI];

that is used in

            forAll(pa, face)
            {
                sumAPtr[pa[face]] -= pCoeffs[face];
            }

Any help is much appreciated...and any discussion would also be helpful.

Dan

[RoE]

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:

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

[chegdan]

Roland,

I ended up changing an OpenFOAM core file in the lduInterface.H by adding the virtual functions:

	    //- 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