[sippanspojk]

Hi foamers,

I am reaching out to you because I need some help with a function object (fO) that I've written. The fO works fine in serial on my local machine but when I submit it to our cluster for parallel computing it does not work and I therefore need help adapting my code (which I believe should be fairly easy if you know what you are doing).

It should be mentioned that I started off by copying the forces fO and adapted that code. Furthermore I am using a modified solver myPimpleFoam which is just a copy of pimpleFoam that I have added a new volScalarField to.

This is my code (or at least the interesting part of it):

void Foam::functionObjects::ESWL::calcESWL()
{
    //- Reading fields
    const volScalarField& P = lookupObject<volScalarField>("p");
    volScalarField& ESWL = lookupObjectRef<volScalarField>("ESWL");
 
    
    //- Initializing variables, vectors and matrices at first time step
    if(std::abs(mesh_.time().value()-mesh_.time().deltaT().value()) < 0.01){

	//- Defining tot number of faces and accessing face normal vectors for later use
	forAll(P.boundaryField(),patchIt){
	    if(patchSet_[patchIt]){
		faceNormals.append(mesh_.boundary()[patchIt].nf());	
		forAll(P.boundaryField()[patchIt], faceIt) {
		    noFaces++;
		}
	    }
	}

	//- Constructing vectors
	meanVec = Foam::List<double>(noFaces,0.0);
	varVec = Foam::List<double>(noFaces,0.0);
	stdVec = Foam::List<double>(noFaces,0.0);   
	M2 = Foam::List<double>(noFaces,0.0);
	
	beta = Foam::List<double>(noFaces, 0.0);
	for(int i=0; i<noFaces; i++){
	    beta[i] = faceNormals[i][0];
	}

	//- Constructing matrices	
	C = Foam::List<Foam::List<double>>(noFaces,Foam::List<double>(noFaces,0.0));
	coVar = C;
	corrMat = C;
    }


    //- Start calculations after 10s to exclude initial uncertainties in statistics
    if(mesh_.time().value() > 10.01){
 
	scalar n = std::round((mesh_.time().value()-10)/mesh_.time().deltaT().value());

	//- Reading face pressure values
	facePressures = Foam::List<double>(noFaces, 0.0);
	int i = 0;
	forAll(P.boundaryField(), patchIt){
	    if(patchSet_[patchIt]){
		forAll(P.boundaryField()[patchIt],faceIt){
		    facePressures[i] = P.boundaryField()[patchIt][faceIt];
		    i ++;	   
		}
	    }
	}


	//- Calculating meanvalues and standard deviation
	Foam::List<double> delta = facePressures - meanVec;
	for(int i=0; i<noFaces; i++){
	    meanVec[i] += delta[i]/n;
	    if(n > 1){
		M2[i] += delta[i]*(facePressures[i]-meanVec[i]);
		varVec[i] = M2[i]/(n-1);
		stdVec[i] = sqrt(varVec[i]);
	    }
	}

	//- Rest of calculations when n > 1 to avoid divition by 0
	if(n > 1){

	    //- Calculating correlation matrix
	    for(int i=0; i<noFaces; i++){			    
		for(int j=0; j<noFaces; j++){
		    C[i][j] += delta[i]*(facePressures[j]-meanVec[j]);
		    coVar[i][j] = C[i][j]/(n-1);
		    if(stdVec[i] == 0 || stdVec[j] == 0){
			corrMat[i][j] = 0.0;
		    }else{
			corrMat[i][j] = coVar[i][j]/stdVec[i]/stdVec[j];
		    }
		}
	    }	

	    //- Calculating nomenator and denomenator
	    Foam::List<double> nomenator(noFaces, 0.0);
	    double denomenator = 0.0;	
	    for(int i = 0; i < noFaces; i++){
		for(int j = 0; j < noFaces; j++){
		    nomenator[i] += corrMat[i][j]*stdVec[j]*beta[j];
		    denomenator += beta[i]*stdVec[i]*corrMat[i][j]*stdVec[j]*beta[j];  
		}
	    }
	    
	    //- Writing values to field 
	    i = 0;
	    forAll(P.boundaryField(),patchIt){
		if(patchSet_[patchIt]){
		    forAll(P.boundaryField()[patchIt],faceIt){
			ESWL.boundaryFieldRef()[patchIt][faceIt] = meanVec[i] + 3.2 stdVec[i]*nomenator[i]/sqrt(denomenator);
		        i ++;
		    }
		}
	    }
	} 
    } 
}

Let me try to explain the basic idea of the fO.

It starts by reading the pressure field. It is then calculating some basic statistics (mean value, standard deviation, variance, etc.) of the pressure values for each cell of some specified patches (patchSet_). In order to do so these lists and lists of lists are stored during the whole simulation and updated for each time step, i.e. I am calculating the statistics on a running stream of data. These variables are therefore defined in the header file. The statistical values are then used for a type of extreme value analysis to calculate a new pressure value for each cell of the specified patches based on the running data set. Lastly these new pressure values are written to a custom made volScalarField "ESWL" that is added to the adapted pimpleFoam solver before the fO is exited and the solver can move on and compute the next time step.

I am imagining that I only need to add a few lines of code in order to get this to run in parallel, since I am just speaking to the solver twice - in the beginning when I read the field values and in the end when I write to the field. My idea is therefore to add a gather in the beginning to gather the pressure values computed on each separate processor and a scatter in the end in order to provide the new updated field to each processor for the next time step.

First of all - have I understood it somewhat correct, would it be enough to add gather/scatter to my code?

Secondly - how do I do this? I have seen example code online and I have read OpenFOAM source code but I have a hard time understanding how I should implement this correctly.


I would be very grateful if someone had the time to help me with this.

Thanks,
David