[Carlen]

Hi Foamers,

I have recently successfully combined two solvers which solve different fields on two separate meshes(one fine and one coarse). Now I need to couple the fields of these two meshes every timestep. How can I achieve this? I read from this thread http://www.cfd-online.com/Forums/ope...-parallel.html #15 that I need to use a portion of the mapFields code. However, I have two questions, where do I put this code in my solver (eg.interFoam.C)? The second question is how can I read

HashTable<word> patchMap;
myDict.lookup("patchMap") >> patchMap;

wordList cuttingPatches;
myDict.lookup("cuttingPatches") >>  cuttingPatches;

How do I define parameters in myDict?
Also, my two meshes are of the same geometry, do you think the method I mentioned above is appropriate?

Thank you in advance
Carlen

[Carlen]

Originally Posted by Carlen

Hi Foamers,

I have recently successfully combined two solvers which solve different fields on two separate meshes(one fine and one coarse). Now I need to couple the fields of these two meshes every timestep. How can I achieve this? I read from this thread http://www.cfd-online.com/Forums/ope...-parallel.html #15 that I need to use a portion of the mapFields code. However, I have two questions, where do I put this code in my solver (eg.interFoam.C)? The second question is how can I read

HashTable<word> patchMap;
myDict.lookup("patchMap") >> patchMap;

wordList cuttingPatches;
myDict.lookup("cuttingPatches") >>  cuttingPatches;

How do I define parameters in myDict?
Also, my two meshes are of the same geometry, do you think the method I mentioned above is appropriate?

Thank you in advance
Carlen

I managed to solver this problem by looking at the cavity case. Now my code for runTime mapFields looks like this:

IOdictionary carlenMapFieldsDict
          (
             IOobject
             (
              "carlenMapFieldsDict",
              runTime.constant(),
              mesh,
              IOobject::MUST_READ_IF_MODIFIED,
              IOobject::NO_WRITE     
            )
          );

          volVectorField U_L(U, mesh);
          volVectorField U_E(U_P, Primary);
          
          volScalarField p_L(p, mesh);
          volScalarField p_E(p_P, Primary);

          volScalarField T_L(T, mesh);
          volScalarField T_E(T_P, Primary);

          HashTable<word> patchMap;
          wordList cuttingPatches;
          carlenMapFieldsDict.lookup("patchMap") >> patchMap;
          carlenMapFieldsDict.lookup("cuttingPatches") >>  cuttingPatches;

          //set up an interpolation method
          meshToMesh::interpolationMethod mapMethod = meshToMesh::imCellVolumeWeight;
          
          // create the meshtomesh interpolator object
          Info << "Interpolating U,p,T Fields" << endl;
          meshToMesh InterpolatorObj
          (
            mesh,
            Primary,
            mapMethod,
            patchMap,
            addProcessorPatches(Primary, cuttingPatches)
          );

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

          InterpolatorObj.mapSrcToTgt
          (
            p,
            eqOp<scalar>(),
            p_P
          );

          InterpolatorObj.mapSrcToTgt
          (
            T,
            eqOp<scalar>(),
            T_P
          );

          Info << "Mapping complete" << endl;

The above code compiles and runs without problem. However, the fields are not mapped between two grids. I narrowed the problem to the way in which I define fields in two meshes. Cna anyone please tell me how to define the fields I want to map?

Thank you very much

carlen

[Sylv]

What do you mean by "the fields are not mapped between two grid"? The solver crashes or the interpolation is wrong/incomplete?

Sylv

[Carlen]

Originally Posted by Sylv

What do you mean by "the fields are not mapped between two grid"? The solver crashes or the interpolation is wrong/incomplete?

Sylv

Hi Marcel,
The interpolation is simply not carried out.

carlen

[Sylv]

Originally Posted by Carlen

Hi Marcel,
The interpolation is simply not carried out.

carlen

So the solver just go through the interpolation without crashing and the target field remained unchanged, right? Are the mesh identical? Is one grid coarser? If yes, which one (target or source)?

Be more detailed in your answers. It is tough to understand the problem with just "not carried out" or "doesn't work"....

Sylv

[Carlen]

Originally Posted by Sylv

So the solver just go through the interpolation without crashing and the target field remained unchanged, right? Are the mesh identical? Is one grid coarser? If yes, which one (target or source)?

Be more detailed in your answers. It is tough to understand the problem with just "not carried out" or "doesn't work"....

Sylv

Hi Marcel,
I managed to make the interpolation work from a fine mesh (mesh) to a coarse (Secondary) mesh (they are of the same geometry) using code below:

       
          IOdictionary carlenMapFieldsDict
          (
             IOobject
             (
              "carlenMapFieldsDict",
              runTime.constant(),
              mesh,
              IOobject::MUST_READ_IF_MODIFIED,
              IOobject::NO_WRITE     
            )
          );

          HashTable<word> patchMap;
          wordList cuttingPatches;
          carlenMapFieldsDict.lookup("patchMap") >> patchMap;
          carlenMapFieldsDict.lookup("cuttingPatches") >>  cuttingPatches;

          //set up an interpolation method
          meshToMesh::interpolationMethod mapMethod = meshToMesh::imCellVolumeWeight;
          
          // create the meshtomesh interpolator

          IOobject Utarget
          (
              "U",
              Secondary.time().timeName(),
              Secondary,
              IOobject::MUST_READ
          );
          volVectorField Utgt(Utarget, Secondary);

          meshToMesh Interp
          (
            mesh,
            Secondary,
            mapMethod,
            patchMap,
            addProcessorPatches(Secondary, cuttingPatches)
          );
         
          Info << "Interpolating U Field" << endl;
          Interp.mapSrcToTgt
          (
            U,
            eqOp<vector>(),
            Utgt
          ); 
          Utgt.write();

 
          Info << "Mapping complete" << endl;

However, the outcome is not ideal. You can see in the picture that the U field mapped on the coarse field is smeared. I wonder if I implement the mapping correctly.

Carlen

[Carlen]

Anyone please help me with this problem. I am stuck.
carlen

[wyldckat]

Hi Carlen,

I have this thread of yours on my to-do list, but I won't be able to look into this any time soon

Nonetheless, please provide a test case and the code ready to be used/tested, so that it takes less time for someone to step in and help.

In addition, you might want to take a look into the source code for this library/toolkit: https://github.com/smenon/dynamicTopoFvMesh - it has the ability to re-mesh the whole domain on-the-fly, although last I saw, it could only handle tetrahedral meshes.

This issue also reminds me of a similar topic: http://www.cfd-online.com/Forums/ope...pdate-fsi.html - you might be able to find some additional helpful ideas on that thread!!

Best regards,
Bruno

[Carlen]

Thank you bruno,
I solved the problem by applying the same initial boundary condition for both the fine and coarse meshes. Now the interpolations is working as intended.

[Maff]

Hi Carlen,

I have been trying to implement mapField during runTime, in a similar way to the one you used (thank you for the useful posts btw).

I stumbled upon the same problem you had: the mapped field is smeared exactly as you posted.

You wrote that you solved the problem by applying the same initial boundary condition. It doesn't seem to work for me. Could you be more specific?

[Maff]

I solved my problem by replacing

eqOp<vector>()

with

plusEqOp<vector>()

I'm not sure why, but it seems to work.

[Ebrahim]

Originally Posted by Maff

I solved my problem by replacing

eqOp<vector>()

with

plusEqOp<vector>()

I'm not sure why, but it seems to work.

Hi Maff,

I have the same issue as yours and I solved it by replacing eqOp<vector>() with plusEqOp<vector>() as well.
Did you find the reason? Actually, I do not know what is changed by this replacement!

/Ebrahim

[Maff]

Does anyone know how to save the meshToMesh addressing, so that it doesn't need to calculate it at every restart?
With some meshes it's a very time consuming process.

I found that with meshToMesh debug set to 1, it saves the srcToTgtCellAddr_, but I would need a way to write the whole result of

          meshToMesh Interp
          (
            mesh,
            Secondary,
            mapMethod,
            patchMap,
            addProcessorPatches(Secondary, cuttingPatches)
          );

so that it can read it if present.
I guess I would need to write my own meshToMesh.C.

Any idea?

[Elliptic CFD]

Originally Posted by Maff

Does anyone know how to save the meshToMesh addressing, so that it doesn't need to calculate it at every restart?
With some meshes it's a very time consuming process.

I found that with meshToMesh debug set to 1, it saves the srcToTgtCellAddr_, but I would need a way to write the whole result of

          meshToMesh Interp
          (
            mesh,
            Secondary,
            mapMethod,
            patchMap,
            addProcessorPatches(Secondary, cuttingPatches)
          );

so that it can read it if present.
I guess I would need to write my own meshToMesh.C.

Any idea?

Which of the meshToMesh libraries in OpenFOAM are you using?

I had the same problem of the cell addressing recalculating at every timestep when I used meshToMesh.H. I now use meshToMesh0.H and the interpolation between two regions is much faster, although, I am clueless on how to get it to work for parallel cases.

By the way, it seems that meshToMesh0.H is only available for the third and fourth versions of OpenFOAM.

[Elliptic CFD]

Originally Posted by Ebrahim

Hi Maff,

I have the same issue as yours and I solved it by replacing eqOp<vector>() with plusEqOp<vector>() as well.
Did you find the reason? Actually, I do not know what is changed by this replacement!

/Ebrahim

In OpenFOAM 3+, mapFields.C (which uses meshToMesh0) has the eqOp in the code while mapFieldsPar (which uses meshToMesh) has the plusEqOp in the code.

The difference may depend on which of the mesh to mesh libraries that you use.

Cheers

[hdotyao]

Hi Carlen,
I am also trying to implement mapFields into a solver. (chtMultiRegionFoam) May I ask where do you put your mapFields code in your solver.c?


Best,
Yao