[mturcios777]

Hello Forum-lurkers!

I was wondering if anyone had any experience with combining diferent dynamicFvMesh features togther to form new libraries. With the recent fix of the lagrangian cloud mapping, I was looking to find a way to combine mesh motion and refinement. I tried something simple first:

1. Create a new class dynamicMotionSolverRefineFvMesh. I basically took the dynamicMotionSolverFvMesh, copied it over and declared it dynamicRefineFvMesh.

2. Change the dynamicMotionSolverRefineFvMesh to do the following:

00062 bool Foam::dynamicMotionSolverRefineFvMesh::update()
00063 {
00064     fvMesh::movePoints(motionPtr_->newPoints());
00065 
00066     if (foundObject<volVectorField>("U"))
00067     {
00068         volVectorField& U =
00069             const_cast<volVectorField&>(lookupObject<volVectorField>("U"));
00070         U.correctBoundaryConditions();
00071     }
00072
00073    Foam::dynamicRefineFvMesh::update()
00074 
00075     return true;
00076 }

This compiles and can be linked to a solver in controlDict, but the solver hangs when it comes time to evolve the cloud. If I call the dynamicRefineFvMesh before the motionSolver, I get the following error, :

Selected 9 cells for refinement out of 3025.
Refined from 3025 to 3088 cells.
Selected 0 split points out of a possible 9.
#0  Foam::error::printStack(Foam::Ostream&) in "/home/gandalf/OpenFOAM/OpenFOAM-2.1.x/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#1  Foam::sigSegv::sigHandler(int) in "/home/gandalf/OpenFOAM/OpenFOAM-2.1.x/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2  Uninterpreted: 
#3  Foam::valuePointPatchField<double>::updateCoeffs() in "/home/gandalf/OpenFOAM/OpenFOAM-2.1.x/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#4  Foam::velocityComponentLaplacianFvMotionSolver::solve() in "/home/gandalf/OpenFOAM/OpenFOAM-2.1.x/platforms/linuxGccDPOpt/lib/libfvMotionSolvers.so"
#5  Foam::motionSolver::newPoints() in "/home/gandalf/OpenFOAM/OpenFOAM-2.1.x/platforms/linuxGccDPOpt/lib/libdynamicMesh.so"
#6  Foam::dynamicMotionSolverRefineFvMesh::update() in "/home/gandalf/OpenFOAM/gandalf-2.1.x/platforms/linuxGccDPOpt/lib/libdynamicMotionSolverRefineFvMesh.so"
#7  
 in "/home/gandalf/OpenFOAM/OpenFOAM-2.1.x/platforms/linuxGccDPOpt/bin/sprayDyMFoam"
#8  __libc_start_main in "/lib/libc.so.6"
#9  
 at /usr/src/packages/BUILD/glibc-2.11.3/csu/../sysdeps/i386/elf/start.S:122
Segmentation fault

Any ideas? Gracias in advance!

[Bernhard]

Are you sure it compiles, dline 73 has a ; missing, that should give an error? Did you try this with only the call to the update() of dynamicRefineFvMesh and does it work?

[mturcios777]

Thanks for looking Bernhard; the missing semicolon is present in the actual code, I just made an error in transcription. Calling the motionsolver by itself works, calling the refinement causes the program to hang at spray evolution. Perhaps I didn't fully initialize the dynamicRefineFvMesh?

[mturcios777]

I think I'm closer to the solution, just not quite sure how to implement it. Turning on some of the debug features, I see that the code is getting stuck in the tracking rescues. When Andy fixed the lagrangian spray code to include mapping, he needed clear the tets used by the tracking algorithm. Maybe the tets aren't being cleared properly and thus the tracking algorithm is getting stuck.

In other words, the approach I have should work for cases that don't use lagrangian particles. I'm going to give this a quick try without the spray and see what happens and report back.

EDIT: So I tried the same library in a case without lagrangian spray (but still with motion and refinement) and was able to do one iteration with the mesh motion before the refinement. Once the second iteration starts I get the same type of error the case with refinement before motion. Digging deeper I found that the exact line the crash occurs at is

 pointMotionU_.boundaryField().updateCoeffs();

Which is inside the solve function for all the fvMotionSolvers. If I have motion happen before refinement and I write every timestep, I get an error when writing the boundary conditions (segmentation fault at write time). It seems that updating the boundary conditions is the key. What I'm not sure this is a case setup problem or a code issue. Any help would be appreciated!

[wyldckat]

Hi Marco,

If I understand you correctly, you have a dynamic mesh that is refined and then moved, and/or vice-versa.

And from the last comment/result you've got, it looks to me like the "pointMotionU_" field has to be reconstructed/remapped before the coefficients can be updated! Otherwise the field is thinking of a mesh that no longer exists...

Unfortunately I'm not very familiar with this part of OpenFOAM (as many other parts), so my question back to you is this: how are all other fields being properly kept between each stage (move/refine)?


Another detail to keep in mind is this (probably you already know this): refining a mesh in the wrong place may lead to an unmovable mesh I hope your test case has enough space for both refining and moving!

As a last note, I took a quick look at how interDyMFoam refined the "damBreakWithObstacle" case and it looked like the mesh was gradually refined and then unrefined as needed... which could provide with some additional ideas for your implementation

Best regards,
Bruno

[mturcios777]

Hey Bruno,

From what I can tell, mapPolyMesh generates the mapping that is used in a host of other calls to map the fields. The mapping itself is done in $FOAM_SRC/OpenFOAM/fields/DimensionedFields/DimensionedField/MapDimensionedFields.H by the function MapDimensionedFields. The pointScalarField pointMotionU_ is not being mapped, even though the volScalarFields are all mapped as well as some surfaceScalarFields. I don't know if pointMotionU_ is not being added to the object registry, or if it is impossible to map a pointScalarField like this.

There shouldn't be any issues for this mesh, its a simple box with a single moving boundary (huge abstraction to what I will eventually need) with a very large cell size to start (this dies on the first timestep when the mesh boundaries have barely moved). I'll upload the test case and some libraries that allow you to edit some of the functionality without "contaminating" the base install.

I'm quickly finding out that dynamic meshes are a whole new level of voodoo in OF...

[mturcios777]

Originally Posted by mturcios777

Hey Bruno,
[...]
I don't know if pointMotionU_ is not being added to the object registry, or if it is impossible to map a pointScalarField like this.

Doing some more digging and with the help of debug messages and doxygen, the key is in the MapFields function of fvMesh, which calls the MapGeometricFields<template args>(mapper) and MapDimensionedFields<template args>(mapper).:

void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
{
    // Create a mapper
    const fvMeshMapper mapper(*this, meshMap);

    // Map all the volFields in the objectRegistry
    MapGeometricFields<scalar, fvPatchField, fvMeshMapper, volMesh>
    (mapper);
    MapGeometricFields<vector, fvPatchField, fvMeshMapper, volMesh>
    (mapper);
    MapGeometricFields<sphericalTensor, fvPatchField, fvMeshMapper, volMesh>
    (mapper);
    MapGeometricFields<symmTensor, fvPatchField, fvMeshMapper, volMesh>
    (mapper);
    MapGeometricFields<tensor, fvPatchField, fvMeshMapper, volMesh>
    (mapper);

    // Map all the surfaceFields in the objectRegistry
    MapGeometricFields<scalar, fvsPatchField, fvMeshMapper, surfaceMesh>
    (mapper);
    MapGeometricFields<vector, fvsPatchField, fvMeshMapper, surfaceMesh>
    (mapper);
    MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
    (mapper);
    MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
    (mapper);
    MapGeometricFields<tensor, fvsPatchField, fvMeshMapper, surfaceMesh>
    (mapper);

    // Map all the dimensionedFields in the objectRegistry
    MapDimensionedFields<scalar, fvMeshMapper, volMesh>(mapper);
    MapDimensionedFields<vector, fvMeshMapper, volMesh>(mapper);
    MapDimensionedFields<sphericalTensor, fvMeshMapper, volMesh>(mapper);
    MapDimensionedFields<symmTensor, fvMeshMapper, volMesh>(mapper);
    MapDimensionedFields<tensor, fvMeshMapper, volMesh>(mapper);

    // Map all the clouds in the objectRegistry
    mapClouds(*this, meshMap);


    const labelList& cellMap = meshMap.cellMap();

    // Map the old volume. Just map to new cell labels.
    if (V0Ptr_)
    {
        scalarField& V0 = *V0Ptr_;

        scalarField savedV0(V0);
        V0.setSize(nCells());

        forAll(V0, i)
        {
            if (cellMap[i] > -1)
            {
                V0[i] = savedV0[cellMap[i]];
            }
            else
            {
                V0[i] = 0.0;
            }
        }
    }

    // Map the old-old volume. Just map to new cell labels.
    if (V00Ptr_)
    {
        scalarField& V00 = *V00Ptr_;

        scalarField savedV00(V00);
        V00.setSize(nCells());

        forAll(V00, i)
        {
            if (cellMap[i] > -1)
            {
                V00[i] = savedV00[cellMap[i]];
            }
            else
            {
                V00[i] = 0.0;
            }
        }
    }
}

Volume, surface and dimensionedFields are all being mapped, but I don't see anything being done for point fields.

I'm not sure how to implement this. What would be the proper template parameters to handle pointScalarFields? I'm thinking something like

MapDimensionedFields<scalar, fvMeshMapper, X>(mapper)

where X is the magic argument

[wyldckat]

Hi Marco,

"pointMesh" seems to be the answer... better yet, it already has a dedicated class and methods for this: http://foam.sourceforge.net/docs/cpp/a06760_source.html -> "src/OpenFOAM/meshes/pointMesh/pointMesh.C"

00035 void Foam::pointMesh::mapFields(const mapPolyMesh& mpm)
00036 {
00037     // Create a mapper
00038     const pointMeshMapper m(*this, mpm);
00039 
00040     MapGeometricFields<scalar, pointPatchField, pointMeshMapper, pointMesh>(m);
00041     MapGeometricFields<vector, pointPatchField, pointMeshMapper, pointMesh>(m);
00042     MapGeometricFields
00043     <
00044         sphericalTensor,
00045         pointPatchField,
00046         pointMeshMapper,
00047         pointMesh
00048     >(m);
00049     MapGeometricFields<symmTensor, pointPatchField, pointMeshMapper, pointMesh>
00050     (m);
00051     MapGeometricFields<tensor, pointPatchField, pointMeshMapper, pointMesh>(m);
00052 }

Although when compared to the list you've posted, this is rather limited... but then again, there isn't much that points are meant to do in OpenFOAM, besides securing the mesh

Among the methods listed in that file, you'll also find "movePoints" and "updateMesh". It almost seems like "pointMotionU_" should be able to update itself, if you can use it properly...

[mturcios777]

Thanks Bruno! That function does look promising, so long as I can figure out how to access it through the dynamicFvMesh (as I only want to do this when there is a topology change).

If I continue to get stumped I'll upload what I have so the group mind can have a look.

EDIT: Some progress. Since I need access to the mapPolyMesh objects generated at refinement, I copied over the dynamicFvRefineMesh classs and added a motionSolver to it in a similar manner to dynamicFvMotionSolver. In the refine() and unrefine() member function, I added the following code right after the meshCutter and polyTopoChange have refined/unrefined the mesh:

pointMesh motionUpdate(*this);
motionUpdate.movePoints(points());
motionUpdate.updateMesh(map);

This compiles and runs, and the mapping function for the pointMotionU fields gets called. However, it is not actually mapped due to the following in MapGeometricField.H:

00118         if (&field.mesh() == &mapper.mesh())
00119         {
00120             if (polyMesh::debug)
00121             {
00122                 Info<< "Mapping " << field.typeName << ' ' << field.name()
00123                     << endl;
00124             }
00125 
00126             // Map the internal field
00127             MapInternalField<Type, MeshMapper, GeoMesh>()
00128             (
00129                 field.internalField(),
00130                 mapper
00131             );
00132 
00133             // Map the patch fields
00134             typename GeometricField<Type, PatchField, GeoMesh>
00135             ::GeometricBoundaryField& bfield = field.boundaryField();
00136             forAll(bfield, patchi)
00137             {
00138                 // Cannot check sizes for patch fields because of
00139                 // empty fields in FV and because point fields get their size
00140                 // from the patch which has already been resized
00141                 //
00142 
00143                 bfield[patchi].autoMap(mapper.boundaryMap()[patchi]);
00144             }
00145 
00146             field.instance() = field.time().timeName();
00147         }
00148         else if (polyMesh::debug)
00149         {
00150             Info<< "Not mapping " << field.typeName << ' ' << field.name()
00151                 << " since originating mesh differs from that of mapper."
00152                 << endl;
00153         }

Instead of mapping I get the "Not mapping" message (because I have debug options on). I think I need to use another object to create the pointMesh, or create it at a different point.

[wyldckat]

Hi Marco,

Since you didn't share the code and test case, I can't tinker with this myself

My guesses are:

• The way you're trying to update the point mesh, somehow it doesn't match properly with the existing mesh pointer.
• But in the end, if "pointMotionU" doesn't want to be remapped, destroy it and build a new one

Best regards,
Bruno

[mturcios777]

Well, I was only going to upload when I was stumped or if I had a working version for the benefit of others. Unfortunately its the first one, so here you go (there isn't really much there except the addition of a mesh motion solver to the dynamicRefineFv class)

I currently have motion before refinement as that helped my determine what the issue was, but they could just as easily be switched so long as the pointMotionU field is mapped properly.

I'm curious about the destroy/create option you mentioned, as doing that is something I've pondered for some time. Would that work with a change of topology? Isn't that pretty much the same thing as mapping?

EDIT: While doing a search of the forums, I noticed this thread: http://www.cfd-online.com/Forums/ope...time-step.html

Would this be an option? Doing topology changes I would need to supply the old and new meshes, and I'm not sure how to do that...

[wyldckat]

Originally Posted by mturcios777

I'm curious about the destroy/create option you mentioned, as doing that is something I've pondered for some time. Would that work with a change of topology? Isn't that pretty much the same thing as mapping?

EDIT: While doing a search of the forums, I noticed this thread: http://www.cfd-online.com/Forums/ope...time-step.html

Would this be an option? Doing topology changes I would need to supply the old and new meshes, and I'm not sure how to do that...

By destroying the pointMotionU field, it could then be created again in the new mesh, being derived from the final mesh, not from a relation between the before and after meshes (as mapFields probably does). Although not very efficient, but optimum is enemy of good.

OK, I'm going to try and play with this neat example Let's see how much of OpenFOAM/C++ voodoo I can conjure up

[mturcios777]

Originally Posted by wyldckat

By destroying the pointMotionU field, it could then be created again in the new mesh, being derived from the final mesh, not from a relation between the before and after meshes (as mapFields probably does). Although not very efficient, but optimum is enemy of good.

OK, I'm going to try and play with this neat example Let's see how much of OpenFOAM/C++ voodoo I can conjure up

Sounds good. The concern I have is that by recreating the field we lose information, since pointMotionU is describing the velocity of all the points of the mesh, and if reset to zero might cause issues with the motionSolver or solver.

[wyldckat]

This is somewhat heavy... OK, the case you've provided requires a couple more libraries, but that's just adding too much entropy to the mix.

I've started creating a case based on the tutorial "mesh/moveDynamicMesh/simpleHarmonicMotion", but I stopped when I noticed that a field is needed on which it bases for refinement. Attached is where I stopped.

This is still going to require a call to setFields and a way to forcefully load the "T" field into memory
That and probably it also needs flux correction, even if the field is basically static.

Nonetheless, this case seems the simplest for isolating how each component works in action.

[wyldckat]

Hi Marco,

Unfortunately I'm unable to do everything I want to do... and although this is an interesting subject I would like to get to understand better, I'm completely out of time to look into this and it's probably over my head experience-wise

I noticed you kept going forward, specially with this bug report: http://www.openfoam.org/mantisbt/view.php?id=638

Good luck!
Bruno

[mturcios777]

Hey everyone!

Just to update the thread; turns out the pointer to the pointMesh was getting deleted before the mapping, and was thus not registered to be mapped in the polyMesh::updateMesh function. The latest version from git should be working now.

I have yet to test this with lagrangian sprays. Stay tuned...

EDIT: So it appears to work with lagrangian clouds. I'm going to play with this shiny new toy. Just so people don't get frustrated if their attempts don't work right off the bat, I had to modify the thermodynamics and turbulence libraries to have them work with dynamic meshes.

[nima3906m]

Dear All

Seems that no one has posted here for a long time. I am facing a problem for while which I think combining mesh motion and mesh refinement can solve but I am not sure if anyone has found a solution for it yet.

I want to simulate to rectangular boxes floating close to each other with small gap in between while the water level is really important for me I need to keep track of bodies motion based on waves hitting them. The problem occurs when I want both at the same time! as I need a really fine mesh in the area between two bodies but this leads to mesh distortion due to high movement levels in just few time steps. I have tried to restrain my bodies movement by designing some springs in between but apparently I should define such high stiffness that makes no sense. I was wondering if anyone can help how I can solve this problem. Can I combine movement of mesh with refinement in this case considering the fact that I have very little knowledge on code writing!