CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

How to compute the integral in the parallel case

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

Like Tree1Likes
  • 1 Post By olesen

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   December 18, 2020, 14:25
Default How to compute the integral in the parallel case
  #1
New Member
 
Saddam Hijazi
Join Date: May 2016
Posts: 12
Rep Power: 7
SaddamH is on a distinguished road
Hi everybody, I have a question regarding the computation of integrals using fvc::domainIntegrate.



How can I compute the value of integrals over the different domains defined by the parallel decomposition ?


For example, I run a case in parallel using 6 processors and I want to compute the following quantity, where Ufield is volVectorField for velocity.



double A =fvc::domainIntegrate(Ufield & fvc::div(linearInterpolate(Ufield) & Ufield.mesh().Sf(), Ufield)).value();


when I then print the value of A using the following command

std::cout << "A is " << A << " from " << "processor" + name(Pstream::myProcNo()) << endl;


I get always - 0.156 from processor 0 - printed six times. This means that the integral is computed only over the first decomposed part of the domain. Ideally, I should be able to sum the values of A given by the computation of the integral across the six subdomains but I am not able to figure out how to do this.


Could you help me in solving this issue ?


Thank you very much,
Saddam
SaddamH is offline   Reply With Quote

Old   December 18, 2020, 15:22
Default
  #2
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,294
Rep Power: 31
olesen will become famous soon enougholesen will become famous soon enough
Quote:
Originally Posted by SaddamH View Post
Hi everybody, I have a question regarding the computation of integrals using fvc::domainIntegrate.

How can I compute the value of integrals over the different domains defined by the parallel decomposition ?

For example, I run a case in parallel using 6 processors and I want to compute the following quantity, where Ufield is volVectorField for velocity.


double A =fvc::domainIntegrate(Ufield & fvc::div(linearInterpolate(Ufield) & Ufield.mesh().Sf(), Ufield)).value();


when I then print the value of A using the following command

std::cout << "A is " << A << " from " << "processor" + name(Pstream::myProcNo()) << endl;


I get always - 0.156 from processor 0 - printed six times. This means that the integral is computed only over the first decomposed part of the domain. Ideally, I should be able to sum the values of A given by the computation of the integral across the six subdomains but I am not able to figure out how to do this.


Could you help me in solving this issue ?


Thank you very much,
Saddam

Use Pout to see information on a per-process basis. Eg,
Code:
Pout<< "Hello from " << Pstream::myProcNo() << nl;
There are different ways to "combine" data from parallel processes. For your particular case, you are looking for a reduction. Eg,
Code:
scalar value = ...;
reduce(value, sumOp<scalar>());
Take a look through the code and you will find plenty of examples.


As a gentle hint, if you are working with OpenFOAM you will generally not use std:: cout for IO, but use OpenFOAM-specific ones. For stdout, it is "Info" which is more or less cout on the master process and /dev/null on all of the others. As soon as you start trying to print out a List or Field of values, you will notice that there are no operators for std ostream, but there are for the OpenFOAM Ostream.
Also do not use "double" and "int", in OpenFOAM they are called "scalar" and "label". It may not seem important at the moment, but if you want code that will also compile with single or mixed-precision and/or 64-bit labelling (for large meshes) it is better to use these OpenFOAM typedefs.
crubio.abujas likes this.
olesen is offline   Reply With Quote

Old   December 18, 2020, 17:37
Default
  #3
New Member
 
Saddam Hijazi
Join Date: May 2016
Posts: 12
Rep Power: 7
SaddamH is on a distinguished road
Quote:
Originally Posted by olesen View Post
Use Pout to see information on a per-process basis. Eg,
Code:
Pout<< "Hello from " << Pstream::myProcNo() << nl;
There are different ways to "combine" data from parallel processes. For your particular case, you are looking for a reduction. Eg,
Code:
scalar value = ...;
reduce(value, sumOp<scalar>());
Take a look through the code and you will find plenty of examples.


As a gentle hint, if you are working with OpenFOAM you will generally not use std:: cout for IO, but use OpenFOAM-specific ones. For stdout, it is "Info" which is more or less cout on the master process and /dev/null on all of the others. As soon as you start trying to print out a List or Field of values, you will notice that there are no operators for std ostream, but there are for the OpenFOAM Ostream.
Also do not use "double" and "int", in OpenFOAM they are called "scalar" and "label". It may not seem important at the moment, but if you want code that will also compile with single or mixed-precision and/or 64-bit labelling (for large meshes) it is better to use these OpenFOAM typedefs.

Thank you very much Mark for your answer and helpful advice. I would like to write the line of codes that I modified according to what you have recommended:


scalar value = fvc::domainIntegrate(L_U_SUPmodes[i] & fvc::div(
linearInterpolate(L_U_SUPmodes[j]) & L_U_SUPmodes[j].mesh().Sf(),
L_U_SUPmodes[k])).value();


reduce(value, sumOp<scalar>());



the last two lines are inside three nested loops which loop over i, j and k, and L_U_SUPmodes is a list of velocity vector fields.


I printed the value of the scalar "value" by doing this:


Pout<< "value from " << Pstream::myProcNo() << " is " << value << nl;



I found it strange that the six values printed (I used six processors for the parallel run) are all the same and equal to 0.157144 so basically the reduce command just gives that value multiplied by the number of processors.



It is worth mention that the correct value of the integral when the code is run in serial is totally different.


Could you help me in identifying the source of this discrepancy.


Thank you very much for your help,
Saddam
SaddamH is offline   Reply With Quote

Old   December 19, 2020, 10:35
Default
  #4
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,294
Rep Power: 31
olesen will become famous soon enougholesen will become famous soon enough
Quote:
Originally Posted by SaddamH View Post
Pout<< "value from " << Pstream::myProcNo() << " is " << value << nl;

I found it strange that the six values printed (I used six processors for the parallel run) are all the same and equal to 0.157144 so basically the reduce command just gives that value multiplied by the number of processors.
You have reached the wrong conclusion. The reduce method uses the operation that you have specified. In this case a sumOp. Thus you have summed the values from all processors, which is what you originally asked.
After the global reduction, the value MUST be the same on all processors, anything else would be a bug.
olesen is offline   Reply With Quote

Old   December 19, 2020, 11:20
Default
  #5
New Member
 
Saddam Hijazi
Join Date: May 2016
Posts: 12
Rep Power: 7
SaddamH is on a distinguished road
Quote:
Originally Posted by olesen View Post
You have reached the wrong conclusion. The reduce method uses the operation that you have specified. In this case a sumOp. Thus you have summed the values from all processors, which is what you originally asked.
After the global reduction, the value MUST be the same on all processors, anything else would be a bug.

Ok I think I misunderstood the purpose of the reduce method. I thought that the integration will be carried out by each processor separately, which would have required a summation step of the values of the integrals. This is not the case according to what I understand from you. However, I still don't understand why the value of the integral computed in serial (which is the right one) differs from the one computed in parallel.


Thank you very much for your help.
SaddamH is offline   Reply With Quote

Old   December 19, 2020, 12:10
Default
  #6
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,294
Rep Power: 31
olesen will become famous soon enougholesen will become famous soon enough
I think that at this stage you need to write out what exactly you are doing. That is to say, if you are doing the reduce inside or outside of the loop etc. This might reveal some mismatch of concepts, and explain why you have a different value from the serial case. You say that the value is different than the serial case, but are not specific. Are they different by 1-2 %, different by 10 times, different by a factor equal to the number of processors? This all makes it incredibly difficult to provide anything remotely helpful.
olesen is offline   Reply With Quote

Old   December 19, 2020, 12:33
Default
  #7
New Member
 
Saddam Hijazi
Join Date: May 2016
Posts: 12
Rep Power: 7
SaddamH is on a distinguished road
Quote:
Originally Posted by olesen View Post
I think that at this stage you need to write out what exactly you are doing. That is to say, if you are doing the reduce inside or outside of the loop etc. This might reveal some mismatch of concepts, and explain why you have a different value from the serial case. You say that the value is different than the serial case, but are not specific. Are they different by 1-2 %, different by 10 times, different by a factor equal to the number of processors? This all makes it incredibly difficult to provide anything remotely helpful.

I am sorry for not being clear. I will write the full method which basically computes a third order tensor (with the usage of c++ the library eigen) by looping over a list of volVectorFields (velocity fields) which is defined and read in the main of my application.




Eigen::Tensor<double, 3> convectiveTensorParallel(label NUmodes,
label NPmodes,
label NSUPmodes)
{
label Csize = NUmodes + NSUPmodes + liftfield.size();
Eigen::Tensor<double, 3> C_tensor;
C_tensor.resize(Csize, Csize, Csize);

for (label i = 0; i < Csize; i++)
{
for (label j = 0; j < Csize; j++)
{
for (label k = 0; k < Csize; k++)
{
scalar value = fvc::domainIntegrate(L_U_SUPmodes[i] & fvc::div(
linearInterpolate(L_U_SUPmodes[j]) & L_U_SUPmodes[j].mesh().Sf(),
L_U_SUPmodes[k])).value();

reduce(value, sumOp<scalar>());
C_tensor(i,j,k) = value;

}
}
}

//if (Pstream:arRun())
//{
// reduce(C_tensor, sumOp<Eigen::Tensor<double, 3>>());
//}


return C_tensor;
}


The important part is the one inside the nested loop, where I compute the value of the integral. The value of the integral for the first loop is 3.548312443e-05 if computed in serial.



In parallel I get the value 0.157144 for the integral and after the reduce command it becomes 0.9428640461 (I use six processors, this latter value is six times the earlier one)


Note that Csize is equal also to the size of the list L_U_SUPmodes.


Thank you very much.
SaddamH is offline   Reply With Quote

Old   December 20, 2020, 17:06
Default
  #8
Member
 
Rodrigo
Join Date: Mar 2010
Posts: 94
Rep Power: 13
guin is on a distinguished road
What about splitting your code into simpler blocks that shall make your life easier when having to debug them? For instance, following lines should just work in a less cryptic manner:

Code:
// ... (after whatever update done for Ufield)

//NOTE: useful when planning to interpolate it later...
Ufield.correctBoundaryConditions(); 

//NOTE: nested just to ensure we don't mess up with your other variable names
{
    const surfaceVectorField Ufield_f = linearInterpolate(Ufield);

    const surfaceScalarField phifield = Ufield_f & Ufield.mesh().Sf();

    const volScalarField UdivU
    (
        Ufield 
      & fvc::div(phifield,"div(phifield)")
    );

    //NOTE: "fvc::domainIntegrate()" is already parallel aware
    const dimensionedScalar A = fvc::domainIntegrate(UdivU); 

    Info<<"your integrated value = "<<A.value()<< endl;
}
guin is offline   Reply With Quote

Reply

Tags
integral computing, parallel code, parallel computation, parallel decomposition

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Sampling for graphs in parallel for AMI case C-L OpenFOAM Bugs 5 January 15, 2018 11:55
implementation of mapFields into parallel transient case simpomann OpenFOAM Pre-Processing 4 August 2, 2016 04:41
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 05:36
[OpenFOAM] parallel case with paraview Ingenieur ParaView 2 August 16, 2013 08:24
Run in parallel a 2mesh case cosimobianchini OpenFOAM Running, Solving & CFD 2 January 11, 2007 06:33


All times are GMT -4. The time now is 04:06.