CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   MeshToMesh parallel (https://www.cfd-online.com/Forums/openfoam-programming-development/110265-meshtomesh-parallel.html)

ganeshv December 6, 2012 10:24

MeshToMesh parallel
 
hi !

I need to perform a mapping from a fine mesh to a coarse mesh every few time steps. The coarse mesh can be a sub-domain of the full fine mesh. This needs to run in parallel. My search for such a utility took me here

http://www.cfd-online.com/Forums/ope...time-step.html

But it's not very conclusive. I thought I could modify the source code of mapFields to do something in this regard. The idea is to read in every processor sub-mesh of the coarse mesh on all the processors. Initiate the desired fields to zero. Map the fields from the fine mesh to the coarse mesh and sum them all up on one processor and write it out. This would essentially become some sort of parallel mapFields utility. I can't get all processors to read in a processor mesh like so


Code:

        Time runTimeCoarseMesh
          (
          Time::controlDictName,
          word("."),
          word(".")/fileName(word("processor")+name(procI))
          );

        fvMesh coarseMesh
          (
          IOobject
          (
            "coarseSolution",
            runTimeCoarseMesh.timeName(),
            runTimeCoarseMesh,
            IOobject::MUST_READ
            )
          );


carowjp February 17, 2013 00:47

meshToMesh in parallel is incorrect
 
1 Attachment(s)
Hello,

Instead of needing to map from a course to a fine mesh I am mapping from a rotated to non-rotated mesh:

Code:

meshToMesh mtmCon(mesh, meshR);
mtmCon.interpolate(alphaR, alpha, meshToMesh::MAP); // map consistent

quaternion Rotation(axis,theta);
tensor R = Rotation.R();
meshPoints = transform(R, meshPoints);
meshR.movePoints(meshPoints);              // rotate the mesh

meshToMesh mtmInc(meshR, mesh, patchMap, cuttingPatches);
mtmInc.interpolate(alpha, alphaR, meshToMesh::INTERPOLATE); // interpolate inconsistent

The cutting patches are read from a mapFieldsDict.

In any case using meshToMesh interpolation in parallel gives incorrect results as shown in the attached image.

The image is a composite of three frames for a case where a scalar field is mapped. The left image is the original, the middle image is rotated and mapped with meshToMesh interpolation in serial (correct), and the right image mapped with meshToMesh interpolation in parallel with the mesh decomposed for 4 processors (2 x 2).

In the parallel case the mapping only takes place on processor0, so the field is 'clipped' by the processor patches.

Can anyone please shed any light on how to address this issue?

thank you,

James

carowjp February 24, 2013 20:36

meshToMesh in parallel is incorrect
 
Sorry for the bump but I am really stumped!

Can anyone please help?

regards,

James

carowjp February 26, 2013 23:56

meshToMesh in parallel maps on processor0 only
 
Hello again,

The meshToMesh interpolation seems to only map on processor0 when run in parallel. I have studied the mapFields utility and attempted to implement the parallel specific portions of code into my application.

With mapFields, the decomposed source and target meshes are read from file:

Code:

   
if (parallelSource && parallelTarget)
{
        //after reading source and target "decomposeParDict" files & determining number of subdomains
       
        List<boundBox> bbsTarget(nProcsTarget);
        List<bool> bbsTargetSet(nProcsTarget, false);

      // loop over processors for source mesh
        for (int procISource=0; procISource<nProcsSource; procISource++)
        {
            Info<< nl << "Source processor " << procISource << endl;

            Time runTimeSource
            (
                Time::controlDictName,
                rootDirSource,
                caseDirSource/fileName(word("processor") + name(procISource))
            );

            #include "setTimeIndex.H"

            fvMesh meshSource
            (
                IOobject
                (
                    sourceRegion,
                    runTimeSource.timeName(),
                    runTimeSource
                )
            );

            Info<< "mesh size: " << meshSource.nCells() << endl;

            boundBox bbSource(meshSource.bounds());

            // loop over processors for target mesh
            for (int procITarget=0; procITarget<nProcsTarget; procITarget++)
            {
                if
                (
                    !bbsTargetSet[procITarget]
                  || (
                      bbsTargetSet[procITarget]
                  && bbsTarget[procITarget].overlaps(bbSource)
                    )
                )
                {
                    Info<< nl << "Target processor " << procITarget << endl;

                    Time runTimeTarget
                    (
                        Time::controlDictName,
                        rootDirTarget,
                        caseDirTarget/fileName(word("processor")
                      + name(procITarget))
                    );

                    fvMesh meshTarget
                    (
                        IOobject
                        (
                            targetRegion,
                            runTimeTarget.timeName(),
                            runTimeTarget
                        )
                    );

                    Info<< "mesh size: " << meshTarget.nCells() << endl;

                    bbsTarget[procITarget] = meshTarget.bounds();
                    bbsTargetSet[procITarget] = true;

                    if (bbsTarget[procITarget].overlaps(bbSource))
                    {
                        if (consistent)
                        {
                            mapConsistentSubMesh
                            (
                                meshSource,
                                meshTarget,
                                mapOrder
                            );
                        }
                        else
                        {
                            mapSubMesh
                            (
                                meshSource,
                                meshTarget,
                                patchMap,
                                addProcessorPatches(meshTarget, cuttingPatches),
                                mapOrder
                            );
                        }
                    }
                }
            }
        }
    }
}

Can anyone please explain how to do this at runtime (without reading from file) by addressing the source and target meshes on each processor?

This was discussed in the thread: Map fields every time step but the issue was not fully resolved.

thanks alot,

James

Lieven April 17, 2013 10:41

Hi Jim,

Did you manage to find a solution for this issue?

Cheers,

L

carowjp April 23, 2013 19:58

meshToMesh in parallel
 
Lieven,

Sorry for the delayed reply I was on a short vacation. :cool:

Unfortunately I have not been able to resolve this issue...but I haven't given up.

best regards,

Jim

pixarzhang April 19, 2014 04:01

meshToMesh:interpolate
 
i am doing something with meshToMesh.interpolate now and encountered a problem. I used "meshToMeshInterp.interpolate(fieldTarget,fieldSour ce)"


,but when i wmake ,the terminal says

"meshandinit.H:71:65: error: no matching function for call to ‘Foam::meshToMesh::interpolate(Foam::volScalarFiel d&, Foam::volScalarField&)’
meshandinit.H:71:65: note: candidates are:
...."

it seems like there must be 4 candidates in my case.something like
{
GeometricField<Type, fvPatchField, volMesh>&,
const GeometricField<Type, fvPatchField, volMesh>&,
order,
const CombineOp& cop = eqOp<Type>()
}
when i looked into meshToMesh.H,it also has the CombineOp& cop

I try to add "meshToMesh::INTERPOLATE "in it, but I do not know "const CombineOp& cop = eqOp<Type>()"can set as what?





i am not very familiar with the source files. I also looked into the source files of mapField but do not find out what the cop is. do any body know where the CombineOp& cop comes from.


by the way,it seems that most people only used 2 candidates in interpolate like

"meshToMeshInterp.interpolate(fieldTarget,fieldSou rce)", but how can't i succeed.below is how i define the variables

fvMesh meshTarget
(
IOobject
(
fvMesh::defaultRegion,
runTimeTarget.timeName(),//runTimeTarget.constant()
runTimeTarget
)
);
meshToMesh meshToMeshInterp(mesh, meshTarget, patchMap_, cuttingPatches_);

volScalarField fieldSource
(
p,
mesh
);

volScalarField fieldTarget
(
p,
meshTarget
);

meshToMeshInterp.interpolate(fieldTarget,fieldSour ce);


any help would be appreciated.

pixarzhang April 19, 2014 07:07

Quote:

Originally Posted by pixarzhang (Post 486998)
i am doing something with meshToMesh.interpolate now and encountered a problem. I used "meshToMeshInterp.interpolate(fieldTarget,fieldSour ce)"


,but when i wmake ,the terminal says

"meshandinit.H:71:65: error: no matching function for call to ‘Foam::meshToMesh::interpolate(Foam::volScalarFiel d&, Foam::volScalarField&)’
meshandinit.H:71:65: note: candidates are:
...."

it seems like there must be 4 candidates in my case.something like
{
GeometricField<Type, fvPatchField, volMesh>&,
const GeometricField<Type, fvPatchField, volMesh>&,
order,
const CombineOp& cop = eqOp<Type>()
}
when i looked into meshToMesh.H,it also has the CombineOp& cop

I try to add "meshToMesh::INTERPOLATE "in it, but I do not know "const CombineOp& cop = eqOp<Type>()"can set as what?





i am not very familiar with the source files. I also looked into the source files of mapField but do not find out what the cop is. do any body know where the CombineOp& cop comes from.


by the way,it seems that most people only used 2 candidates in interpolate like

"meshToMeshInterp.interpolate(fieldTarget,fieldSou rce)", but how can't i succeed.below is how i define the variables

fvMesh meshTarget
(
IOobject
(
fvMesh::defaultRegion,
runTimeTarget.timeName(),//runTimeTarget.constant()
runTimeTarget
)
);
meshToMesh meshToMeshInterp(mesh, meshTarget, patchMap_, cuttingPatches_);

volScalarField fieldSource
(
p,
mesh
);

volScalarField fieldTarget
(
p,
meshTarget
);

meshToMeshInterp.interpolate(fieldTarget,fieldSour ce);


any help would be appreciated.




I tried it out finally!!
the fourh candidate can set as " Foam::eqOp<scalar>() ",and the .C file can be compiled successfully.

despite i don't know where the class comes from ,i searched the src fold,but can't find any file named *eqOp* or related.

wish someone could tell me where to know what Foam::eqOp<scalar>() is .And all the questions before are also not resolved.

wyldckat April 19, 2014 07:53

Quick answer: "Foam::eqOp<Type>()" is a special class with a single operator-method implemented, meant for helping with the process of combining data from several parallel processes.
By assigning a specific variable type in "Type", you specify which kind of data is being compared.

pixarzhang April 20, 2014 12:42

Thank you for your explaination

but I am still not very clear.
I found a statement in meshToMeshinterplate.c which is "
const CombineOp& cop
...
cop(toF[celli],fromf[adr[celli]]) ",
I guess it give the value of fromf[adr[celli]] to toF[celli]?

do you know where the source file of the class is located, or where I can find it,I want to see more details if could.

wyldckat April 20, 2014 13:41

Hi pixarzhang,

Which OpenFOAM version are you using? Because I'm looking at OpenFOAM 2.3.x and the code has changed considerably from your description. Therefore knowing which version of OpenFOAM you're using would make it a lot easier to explain this.

Best regards,
Bruno

pixarzhang April 21, 2014 03:58

HI,wyldckat,
my OpenFOAM version is 2.2.0,I haven't updated it.
thanks again.

wyldckat April 26, 2014 07:35

Hi pixarzhang,

I've finally managed to look into this. So what happens is this:
  1. In OpenFOAM 2.1 it uses this line of code:
    Code:

    toF[celli] = fromVf[adr[celli]];
  2. In OpenFOAM 2.2, it uses this line of code:
    Code:

    cop(toF[celli], fromVf[adr[celli]]);
  3. The first one would simply state that the value "toF[celli]" on the current cell mesh "celli", should be the value that was on "fromVf[adr[celli]]" for the old cell address "adr[celli]".
  4. The limitation of the first line of code is that it requires for both cell values (old and new locations) have to be in the same processor (aka sub-domain).
  5. The second implementation, done in OpenFOAM 2.2, allows for the value to "fromVf[adr[celli]]" to come from another processor (aka sub-domain), which is why it is using "cop" as a combine operator, instead of standard equal operator.
This implementation apparently worked, but on OpenFOAM 2.3 they changed everything completely for this class "meshToMesh", apparently using a more robust and parallel way of doing its objectives.

Best regards,
Bruno

pixarzhang April 27, 2014 03:25

thanks,that helps a lot

Sylv April 29, 2014 05:49

To complete the answer of wyldckat, the meshtomesh class in 2.3.X is the only one (to my knowledge), which can permform the interpolation from a decomposted source mesh (srcMesh) to a decomposted target mesh (tgtMesh).

In my own solver, I interpolate the U field from a coarse source mesh (srcMesh) to a fine target mesh (tgtMesh). The tgtMesh represent only a small portion of the srcMesh. Here is the code, which does the job:
Code:


// a simple function
wordList addProcessorPatches
(
    const fvMesh& meshTarget,
    const wordList& cuttingPatches
)
{
    // Add the processor patches to the cutting list
    HashSet<word> cuttingPatchTable;
    forAll(cuttingPatches, i)
    {
        cuttingPatchTable.insert(cuttingPatches[i]);
    }

    const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();

    forAll(pbm, patchI)
    {
        if (isA<processorPolyPatch>(pbm[patchI]))
        {
            const word& patchName = pbm[patchI].name();
            cuttingPatchTable.insert(patchName);
        }
    }

    return cuttingPatchTable.toc();
}

// generate the source mesh "srcMesh", the target mesh "tgtMesh",
// the source velocity field "U"  and the target velocity field
// "Utgt" (fill it with vector::zero or some dummy value...)

A very long piece of code here. You should know how to do it....

// get patchMap (I read the hash table from a dict...)
HashTable<word> patchMap;
myDict.lookup("patchMap") >> patchMap;

// get the cutting patch. This the target patches which cut the source mesh
// I read this list from a dictionnary
wordList cuttingPatches;
myDict.lookup("cuttingPatches") >>  cuttingPatches;

// As I want to be able to run in parallel, the extra processor patches of the target
// mesh are also cutting patch, therefore, I need to update "cuttingPatch"
wordList allCuttingPatches;
allCuttingPatches = addProcessorPatches(tgtMesh, cuttingPatches)

//set a map method
meshToMesh::interpolationMethod mapMethod = meshToMesh::imCellVolumeWeight;

// create the meshtomesh interpolator object
meshToMesh interpolatorObj
(
            srcMesh,
            tgtMesh,
            mapMethod,
            patchMap,
            allCuttingPatches
);

//do the interpolation between U of the source mesh and Utgt of the target mesh
interpolatorObj.mapSrcToTgt
        (
            U,
            eqOp<vector>(),
            Utgt
        );

I hope it will help

Marcel

pixarzhang May 2, 2014 08:27

it does help.thanks~

chriss85 July 16, 2014 09:48

I'm now facing a similar task.
My solver uses a complete mesh that is then separated into different regions.
In the solver, I need to map a field between the full mesh and the separate regions.

I don't care about the boundary conditions of the fields as I can set them manually in the solver to zeroGradient. Do I need a patchMap and cuttingPatches then, or can I leave these lists empty (apart from the processor patches as you've shown before)?
If I leave them empty, do I need to manually assign a zeroGradient BC, or is it kept from the original field? If it is, a call to .correctBoundaryConditions() might be required.

I'm using OF 2.3.x, is the code you posted above valid for this version and parallel runs?

Also, for performance reasons, can we cache the mesh-to-mesh addressing? In my first tries the mapping has been quite slow.

openfoam_aero December 21, 2023 13:58

Quote:

Originally Posted by Sylv (Post 488849)
To complete the answer of wyldckat, the meshtomesh class in 2.3.X is the only one (to my knowledge), which can permform the interpolation from a decomposted source mesh (srcMesh) to a decomposted target mesh (tgtMesh).

In my own solver, I interpolate the U field from a coarse source mesh (srcMesh) to a fine target mesh (tgtMesh). The tgtMesh represent only a small portion of the srcMesh. Here is the code, which does the job:
Code:


// a simple function
wordList addProcessorPatches
(
    const fvMesh& meshTarget,
    const wordList& cuttingPatches
)
{
    // Add the processor patches to the cutting list
    HashSet<word> cuttingPatchTable;
    forAll(cuttingPatches, i)
    {
        cuttingPatchTable.insert(cuttingPatches[i]);
    }

    const polyBoundaryMesh& pbm = meshTarget.boundaryMesh();

    forAll(pbm, patchI)
    {
        if (isA<processorPolyPatch>(pbm[patchI]))
        {
            const word& patchName = pbm[patchI].name();
            cuttingPatchTable.insert(patchName);
        }
    }

    return cuttingPatchTable.toc();
}

// generate the source mesh "srcMesh", the target mesh "tgtMesh",
// the source velocity field "U"  and the target velocity field
// "Utgt" (fill it with vector::zero or some dummy value...)

A very long piece of code here. You should know how to do it....

// get patchMap (I read the hash table from a dict...)
HashTable<word> patchMap;
myDict.lookup("patchMap") >> patchMap;

// get the cutting patch. This the target patches which cut the source mesh
// I read this list from a dictionnary
wordList cuttingPatches;
myDict.lookup("cuttingPatches") >>  cuttingPatches;

// As I want to be able to run in parallel, the extra processor patches of the target
// mesh are also cutting patch, therefore, I need to update "cuttingPatch"
wordList allCuttingPatches;
allCuttingPatches = addProcessorPatches(tgtMesh, cuttingPatches)

//set a map method
meshToMesh::interpolationMethod mapMethod = meshToMesh::imCellVolumeWeight;

// create the meshtomesh interpolator object
meshToMesh interpolatorObj
(
            srcMesh,
            tgtMesh,
            mapMethod,
            patchMap,
            allCuttingPatches
);

//do the interpolation between U of the source mesh and Utgt of the target mesh
interpolatorObj.mapSrcToTgt
        (
            U,
            eqOp<vector>(),
            Utgt
        );

I hope it will help

Marcel


Since it has been a few years and since this thread MASSIVELY helped me, I wish to contribute.

There are a few syntax changes that have taken place over the years. Particularly, the handling of the decomposition of the source and target meshes differently on the processors is improved thanks to the LOD and AABB methods.


I use Openfoam v1812 so please bear this in mind!

In light of that, we can change the code to

Code:

    meshToMesh InterpolatorObj
    (
    mesh,
    target,
    Foam::meshToMesh::interpolationMethod::imCellVolumeWeight,
    Foam::meshToMesh::procMapMethod::pmAABB,
    false
    );

    //do the interpolation between U of the source mesh and Utgt of the target mesh
    InterpolatorObj.mapSrcToTgt
    (
    U,
    plusEqOp<vector>(),
    U_target
    );

Here the only thing that has changed is that we add the procMapMethod pmAABB and false for interpolating all patches (I presume this is bnecause my case is not consistent so it would make sense to not interpolate all patches).

Hope this is helpful for someone who is working on this now!

I am calling this piece of code at the end of the simple loop (I am running a steady state case and I do not need the interpolations at every iteration). But as someone pointed out here, it would be nice to know how to have the cell addressing computed once and carry forward.


All times are GMT -4. The time now is 16:42.