[itsme_kit]

Hi OpenFOAMers,

My problem is similar to these two threads:

How to get cell indices and manipulate field data

runTime mapFields from fineMesh to coarseMesh

I want to do data exchange for internal field values between two identical meshes with different mesh density.

Each mesh has own solver. For example, solver1 gets results of velocity and sends them to solver2 for temperature calculation. And temperature is returned to solver1 for velocity calculation, the loop continues like this.

Those two threads make use of the code of mapfields and implement it in solver.

I'm not sure if I need to merge those two solvers as a new solver before using mapfields code in the solver.

Where to the put the code in the solver?

Similar work in OpenFOAM has been found in journal paper, it states that "The computations on the two meshes are fully
coupled, meaning that data (temperature from macro- to meso-,
and solid fraction from meso- to macroscale) is exchanged in both
directions at every time step via a linear interpolation framework"

I'm not sure how the author implement the data exchange, limited information provided in the paper.

Thanks in advance.

[ch1]

Hello,
with "two identical meshes with different mesh density" you mean that they have the same geometry but different mesh refinement, right? Which version of OpenFOAM do you use?

In my thread, which you refer to, I used foam-extend 4.0. So for that I can give you some hints.

To use a second mesh I do the following (code can be inserted after e.g. #include "createTime.H" ):

        
    fvMesh emMesh
    (
        Foam::IOobject
        (
            word("emRegion"),
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

    fvMesh mesh //keep default variable name to avoid compilation errors with the default fluid solvers
    (
        Foam::IOobject
        (
            word("fluidRegion"),
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

This implies that the case folder has got roughly the following structure.

../caseName/
|
|-- 0/
| |-- emRegion/ --> initial conditions T, U, p,...
| |-- fluidRegion/ --> initial conditions T, U, p,...
|
|-- const/
| |-- emRegion/polyMesh/ --> emMesh
| |-- fluidRegion/polyMesh/ --> mesh (fluidMesh)
|-- system/
| |-- emRegion/ --> fvSolution, fvSchemes
| |-- fluidRegion/ --> fvSolution, fvSchemes

In your case, the "emMesh" is where you want to solve the temperature equation and the default "mesh" is where the fluid equations are solved.

In createFields.H of your solver, you should assign the fields to the corresponding meshes.
A field in the e.g. "emMesh" would be like

    volScalarField A
    (
        IOobject
        (
            "A",
            runTime.timeName(),
            emMesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        emMesh
    );

The standard fluid fields eventually may remain unchanged like

    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

However, if a field happens to appear in both meshes and you want to map that field between the meshes, I recommend to do the following:

    volScalarField Tfluid
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    volScalarField Tem
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            emMesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        emMesh
    );

So it will be clear from the code, to which mesh the variable is assigned.


The last thing you need to know is how to properly map + interpolate the field between the meshes at run time.

This code can be executed after #include "createTime.H" as well:

// read the mapFieldsDict located in /case/system
IOdictionary mapFieldsDict
(
	IOobject
	(
		"mapFieldsDict",
		runTime.system(),
		runTime,
		IOobject::MUST_READ,
		IOobject::NO_WRITE,
		false
	)
);

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

Just create an empty mapFieldsDict file and put it into the system directory:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      mapFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// List of pairs of source/target patches for mapping
patchMap
(
);

// List of target patches cutting the source domain (these need to be
// handled specially e.g. interpolated from internal values)
cuttingPatches
(
);

And here is a code example of map + interpolation, which will be executed at runtime:

// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh(); 
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();     

// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
    meshSource,
    meshTarget,
    patchMap,
    cuttingPatches
);

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

Info<< nl
    << "Mapping fields for time " << meshSource.time().timeName()
    << nl << endl;

// Interpolate field
meshToMeshInterp.interpolate
(
    Tfluid, //target field
    Tem, //source field
    meshToMesh::INTERPOLATE
);

I guess that's everything you need to get started.
Note, that I encountered some "undesired extrapolation bug", which I mentioned previously here in this thread and which I could not resolve yet. But this problem should not occur in your case, because you mentioned that your geometry is consistent.

Best regards,
Christian

[frobaux]

Wow,

Thank you for that complete explanation, I wish I had it when I started to write mine!
Everything is in there!

I just want to add two things.

First, I use

meshToMeshInterp.mapSrcToTgt(srcField,  plusEqOp<scalar>(), tgtField)

instead of "meshToMeshInterp.interpolate(..)" (Probably different version, I am on V17.12+)
Second, I am pretty sure that it is not what we call fully coupled.

Correct me if I am wrong, but fully coupled means that only one matrice is inverted and describe the entire problem. like Ax=B with x beeing ALL unknowns (both T in the coarse mesh and solid fraction in the finer). On the contrary, a normal coupling would have a loop solving for one then for the other, until convergence is reached.

But I would have no clue on how to implement a fully coupling in OpenFOAM.

[ch1]

Originally Posted by frobaux

Wow,

Thank you for that complete explanation, I wish I had it when I started to write mine!
Everything is in there!

I just want to add two things.

First, I use

meshToMeshInterp.mapSrcToTgt(srcField,  plusEqOp<scalar>(), tgtField)

instead of "meshToMeshInterp.interpolate(..)" (Probably different version, I am on V17.12+)
Second, I am pretty sure that it is not what we call fully coupled.

Correct me if I am wrong, but fully coupled means that only one matrice is inverted and describe the entire problem. like Ax=B with x beeing ALL unknowns (both T in the coarse mesh and solid fraction in the finer). On the contrary, a normal coupling would have a loop solving for one then for the other, until convergence is reached.

But I would have no clue on how to implement a fully coupling in OpenFOAM.

Yes, this command might be different due to different OF versions. I am using foam-extend 4.0 and if you want to implement fully-coupled solvers you can find it there too. The library is called fvBlockMatrix.H. See also pUCoupledFOAM in foam-extend to get started.

I don't know what itsme_kit is doing with the field variables or how the equations look like. So maybe what was meant with "fully coupled" is that the coupling is bidirectional but if the coupling is nonlinear you cannot solve it as a fully coupled matrix, so a segregated approach is needed.

[itsme_kit]

Hi Christian

Thanks for your detailed explanation.

I'm back to CFD recently.

I will look into your suggestion immediately and let you know what I got.

Thanks a million again!

Originally Posted by ch1

Yes, this command might be different due to different OF versions. I am using foam-extend 4.0 and if you want to implement fully-coupled solvers you can find it there too. The library is called fvBlockMatrix.H. See also pUCoupledFOAM in foam-extend to get started.

I don't know what itsme_kit is doing with the field variables or how the equations look like. So maybe what was meant with "fully coupled" is that the coupling is bidirectional but if the coupling is nonlinear you cannot solve it as a fully coupled matrix, so a segregated approach is needed.

[itsme_kit]

Hi Christian

I completely followed your procedure, it is fully working after several trials.

Thanks very much.

There is one potential shortcoming in this procedure I wish to discuss:
Solver2 has to wait for execution until solver1 completes its job.

I wonder if we could execute these two solvers simultaneously, for example, solver1 solves and maps fields into solver2 and solver2 solves fields at each time step.

Looking forward to hearing from you.

Best Regards,

Ke

Originally Posted by ch1

Yes, this command might be different due to different OF versions. I am using foam-extend 4.0 and if you want to implement fully-coupled solvers you can find it there too. The library is called fvBlockMatrix.H. See also pUCoupledFOAM in foam-extend to get started.

I don't know what itsme_kit is doing with the field variables or how the equations look like. So maybe what was meant with "fully coupled" is that the coupling is bidirectional but if the coupling is nonlinear you cannot solve it as a fully coupled matrix, so a segregated approach is needed.

[ch1]

Hey Ke,

sorry for my late response. So if I understand correctly you mean an implementation like this:

while (runTime.run())
{

eqnSolverA();

mapFieldsFromAtoB();

eqnSolverB();

}

Actually, this should work without a problem if you integrate the code of your two solvers, which at the moment are separate applications if I got you correctly, into a single application, which will be your new combined solver.

Maybe you can share your code one more time, so I can better understand the problem and what you are trying to achieve.

Best regards,
Christian

[itsme_kit]

Hi Christian

I have completed the mapping between two solvers in one code as you suggested.

However, my case needs a little bit more complexity, for example, solverA runs at timestep T=0.0001, and solverB runs at timestep T=0.01, so the mapping should be implemented at every T=0.01.

I really appreciate any idea of doing mapping in this way.

Many thanks

Originally Posted by ch1

Hey Ke,

sorry for my late response. So if I understand correctly you mean an implementation like this:

while (runTime.run())
{

eqnSolverA();

mapFieldsFromAtoB();

eqnSolverB();

}

Actually, this should work without a problem if you integrate the code of your two solvers, which at the moment are separate applications if I got you correctly, into a single application, which will be your new combined solver.

Maybe you can share your code one more time, so I can better understand the problem and what you are trying to achieve.

Best regards,
Christian

[randolph]

Originally Posted by ch1

Hello,
with "two identical meshes with different mesh density" you mean that they have the same geometry but different mesh refinement, right? Which version of OpenFOAM do you use?

In my thread, which you refer to, I used foam-extend 4.0. So for that I can give you some hints.

To use a second mesh I do the following (code can be inserted after e.g. #include "createTime.H" ):

        
    fvMesh emMesh
    (
        Foam::IOobject
        (
            word("emRegion"),
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

    fvMesh mesh //keep default variable name to avoid compilation errors with the default fluid solvers
    (
        Foam::IOobject
        (
            word("fluidRegion"),
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

This implies that the case folder has got roughly the following structure.

../caseName/
|
|-- 0/
| |-- emRegion/ --> initial conditions T, U, p,...
| |-- fluidRegion/ --> initial conditions T, U, p,...
|
|-- const/
| |-- emRegion/polyMesh/ --> emMesh
| |-- fluidRegion/polyMesh/ --> mesh (fluidMesh)
|-- system/
| |-- emRegion/ --> fvSolution, fvSchemes
| |-- fluidRegion/ --> fvSolution, fvSchemes

In your case, the "emMesh" is where you want to solve the temperature equation and the default "mesh" is where the fluid equations are solved.

In createFields.H of your solver, you should assign the fields to the corresponding meshes.
A field in the e.g. "emMesh" would be like

    volScalarField A
    (
        IOobject
        (
            "A",
            runTime.timeName(),
            emMesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        emMesh
    );

The standard fluid fields eventually may remain unchanged like

    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

However, if a field happens to appear in both meshes and you want to map that field between the meshes, I recommend to do the following:

    volScalarField Tfluid
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    volScalarField Tem
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            emMesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        emMesh
    );

So it will be clear from the code, to which mesh the variable is assigned.


The last thing you need to know is how to properly map + interpolate the field between the meshes at run time.

This code can be executed after #include "createTime.H" as well:

// read the mapFieldsDict located in /case/system
IOdictionary mapFieldsDict
(
	IOobject
	(
		"mapFieldsDict",
		runTime.system(),
		runTime,
		IOobject::MUST_READ,
		IOobject::NO_WRITE,
		false
	)
);

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

Just create an empty mapFieldsDict file and put it into the system directory:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      mapFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// List of pairs of source/target patches for mapping
patchMap
(
);

// List of target patches cutting the source domain (these need to be
// handled specially e.g. interpolated from internal values)
cuttingPatches
(
);

And here is a code example of map + interpolation, which will be executed at runtime:

// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh(); 
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();     

// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
    meshSource,
    meshTarget,
    patchMap,
    cuttingPatches
);

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

Info<< nl
    << "Mapping fields for time " << meshSource.time().timeName()
    << nl << endl;

// Interpolate field
meshToMeshInterp.interpolate
(
    Tfluid, //target field
    Tem, //source field
    meshToMesh::INTERPOLATE
);

I guess that's everything you need to get started.
Note, that I encountered some "undesired extrapolation bug", which I mentioned previously here in this thread and which I could not resolve yet. But this problem should not occur in your case, because you mentioned that your geometry is consistent.

Best regards,
Christian

Christian,

I follow your approach (trying to define a volScalarField on the second mesh) but an error complains that the emMesh was not declared in this scope. Any suggestions?

Thanks,

volScalarField Tem
(
IOobject
(
"T",
runTime.timeName(),
emMesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
emMesh
);

[Luke99]

Originally Posted by randolph

Christian,

I follow your approach (trying to define a volScalarField on the second mesh) but an error complains that the emMesh was not declared in this scope. Any suggestions?

Thanks,

volScalarField Tem
(
IOobject
(
"T",
runTime.timeName(),
emMesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
emMesh
);

I am currently facing the same issue. Did someone find a solution to this yet? It also happens if you just want to rename "mesh" to something else.


EDIT:
I might be wrong but from what I understand it would mean that I would have to declare emMesh (or whatever you call your 2nd mesh) outside of the main function. Yet the question is where and how does this happen for mesh, that might give insights into how to fix this issue. For the record I encountered this problem while trying to add a second dynamic mesh to interfoam, but even a static mesh throws the same error.

[Luke99]

For everyone facing the same issue:
As we know it is possible to have multiple meshes in solvers like chtMultiRegionFoam i poked and modified it and came to following result:


if you create your multpile meshes in createMeshes.H(or wherever you create them) as given above or in the same manner as done in chtMultiRegionFoam just add the following line at the start of your main loop:

#define CREATE_MESH createMeshesPostProcess.H

inside createMeshesPostProcess.H or however you might call it, write:

#include "createMeshes.H"

I don`t know if this is a proper solution, but it works for my solver in OpenFOAM2212
This also works for one of the meshes being a dynamicFvMesh& insted of fvMesh&. I haven´t tested if it will work for more though, but i don't see a problem there.


Edit: I realised that runtime data transfer would be too time consuming so I will not be following up on it