CFD Online Discussion Forums - [snappyHexMesh] SnappyHexMeshing SnakeRiverCanyon tutorial

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Meshing & Mesh Conversion (https://www.cfd-online.com/Forums/openfoam-meshing/)

- - [snappyHexMesh] SnappyHexMeshing SnakeRiverCanyon tutorial (https://www.cfd-online.com/Forums/openfoam-meshing/76544-snappyhexmeshing-snakerivercanyon-tutorial.html)

SnappyHexMeshing SnakeRiverCanyon tutorial

Hi,

I´ve been trying to use snappyHexMesh on the SnakeRiverCanyon tutorial using the instructions posted by AlanR. My snappy file runs and creates output files but no ´snap´ of the mesh appears to have taken place.

I have inserted my snappyHexMeshDict below and would really appreciate any pointers or hints what I might be doing wrong.

Thanks for your help.

Thom

Code:

FoamFile

{

    version     2.0;

    format      ascii;

    class       dictionary;

    object      snappyHexMeshDict;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



// Which of the steps to run

castellatedMesh true;

snap            true;

addLayers       true;





// Geometry. Definition of all surfaces. All surfaces are of class

// searchableSurface.

// Surfaces are used

// - to specify refinement for any mesh cell intersecting it

// - to specify refinement for any mesh cell inside/outside/near

// - to 'snap' the mesh boundary to the surface

geometry

{

    AcrossRiver.stl

    {

        type triSurfaceMesh;

        name AcrossRiver;

    }



    refinementBox

    {

        type searchableBox;

        min (659531  4.7513e+06   1028);

        max (662381  4.75454e+06  1200);

    }

};







// Settings for the castellatedMesh generation.

castellatedMeshControls

{



    // Refinement parameters

    // ~~~~~~~~~~~~~~~~~~~~~



    // If local number of cells is >= maxLocalCells on any processor

    // switches from from refinement followed by balancing

    // (current method) to (weighted) balancing before refinement.

    maxLocalCells 1000000;



    // Overall cell limit (approximately). Refinement will stop immediately

    // upon reaching this number so a refinement level might not complete.

    // Note that this is the number of cells before removing the part which

    // is not 'visible' from the keepPoint. The final number of cells might

    // actually be a lot less.

    maxGlobalCells 2000000;



    // The surface refinement loop might spend lots of iterations refining just a

    // few cells. This setting will cause refinement to stop if <= minimumRefine

    // are selected for refinement. Note: it will at least do one iteration

    // (unless the number of cells to refine is 0)

    minRefinementCells 10;



    // Allow a certain level of imbalance during refining

    // (since balancing is quite expensive)

    // Expressed as fraction of perfect balance (= overall number of cells /

    // nProcs). 0=balance always.

    maxLoadUnbalance 0.10;





    // Number of buffer layers between different levels.

    // 1 means normal 2:1 refinement restriction, larger means slower

    // refinement.

    nCellsBetweenLevels 3;







    // Explicit feature edge refinement

    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



    // Specifies a level for any cell intersected by its edges.

    // This is a featureEdgeMesh, read from constant/triSurface for now.

    features

    (

        //{

        //    file "someLine.eMesh";

        //    level 2;

        //}

    );







    // Surface based refinement

    // ~~~~~~~~~~~~~~~~~~~~~~~~



    // Specifies two levels for every surface. The first is the minimum level,

    // every cell intersecting a surface gets refined up to the minimum level.

    // The second level is the maximum level. Cells that 'see' multiple

    // intersections where the intersections make an

    // angle > resolveFeatureAngle get refined up to the maximum level.



    refinementSurfaces

    {

        AcrossRiver

        {

            // Surface-wise min and max refinement level

            level (2 2);

        }

    }



    // Resolve sharp angles

    resolveFeatureAngle 30;





    // Region-wise refinement

    // ~~~~~~~~~~~~~~~~~~~~~~



    // Specifies refinement level for cells in relation to a surface. One of

    // three modes

    // - distance. 'levels' specifies per distance to the surface the

    //   wanted refinement level. The distances need to be specified in

    //   descending order.

    // - inside. 'levels' is only one entry and only the level is used. All

    //   cells inside the surface get refined up to the level. The surface

    //   needs to be closed for this to be possible.

    // - outside. Same but cells outside.



    refinementRegions

    {

        refinementBox

        {

            mode inside;

            levels ((1E15 4));

        }

    }





    // Mesh selection

    // ~~~~~~~~~~~~~~



    // After refinement patches get added for all refinementSurfaces and

    // all cells intersecting the surfaces get put into these patches. The

    // section reachable from the locationInMesh is kept.

    // NOTE: This point should never be on a face, always inside a cell, even

    // after refinement.

    locationInMesh (659535 4.7513e+06 1328);

}







// Settings for the snapping.

snapControls

{

    //- Number of patch smoothing iterations before finding correspondence

    //  to surface

    nSmoothPatch 3;



    //- Relative distance for points to be attracted by surface feature point

    //  or edge. True distance is this factor times local

    //  maximum edge length.

    tolerance 4.0;



    //- Number of mesh displacement relaxation iterations.

    nSolveIter 30;



    //- Maximum number of snapping relaxation iterations. Should stop

    //  before upon reaching a correct mesh.

    nRelaxIter 5;

}







// Settings for the layer addition.

addLayersControls

{

    // Are the thickness parameters below relative to the undistorted

    // size of the refined cell outside layer (true) or absolute sizes (false).

    relativeSizes true;



    // Per final patch (so not geometry!) the layer information

    layers

    {

        AcrossRiver_patch0

        {

            nSurfaceLayers 1;

        }

    }



    // Expansion factor for layer mesh

    expansionRatio 1.0;



    //- Wanted thickness of final added cell layer. If multiple layers

    //  is the thickness of the layer furthest away from the wall.

    //  See relativeSizes parameter.

    finalLayerThickness 0.3;



    //- Minimum thickness of cell layer. If for any reason layer

    //  cannot be above minThickness do not add layer.

    //  Relative to undistorted size of cell outside layer.

    minThickness 0.1;



    //- If points get not extruded do nGrow layers of connected faces that are

    //  also not grown. This helps convergence of the layer addition process

    //  close to features.

    nGrow 1;





    // Advanced settings



    //- When not to extrude surface. 0 is flat surface, 90 is when two faces

    //  make straight angle.

    featureAngle 30;



    //- Maximum number of snapping relaxation iterations. Should stop

    //  before upon reaching a correct mesh.

    nRelaxIter 3;



    // Number of smoothing iterations of surface normals

    nSmoothSurfaceNormals 1;



    // Number of smoothing iterations of interior mesh movement direction

    nSmoothNormals 3;



    // Smooth layer thickness over surface patches

    nSmoothThickness 10;



    // Stop layer growth on highly warped cells

    maxFaceThicknessRatio 0.5;



    // Reduce layer growth where ratio thickness to medial

    // distance is large

    maxThicknessToMedialRatio 0.3;



    // Angle used to pick up medial axis points

    minMedianAxisAngle 130;



    // Create buffer region for new layer terminations

    nBufferCellsNoExtrude 0;



    // Overall max number of layer addition iterations

    nLayerIter 50;

}







// Generic mesh quality settings. At any undoable phase these determine

// where to undo.

meshQualityControls

{

    //- Maximum non-orthogonality allowed. Set to 180 to disable.

    maxNonOrtho 65;



    //- Max skewness allowed. Set to <0 to disable.

    maxBoundarySkewness 20;

    maxInternalSkewness 4;



    //- Max concaveness allowed. Is angle (in degrees) below which concavity

    //  is allowed. 0 is straight face, <0 would be convex face.

    //  Set to 180 to disable.

    maxConcave 80;



    //- Minimum projected area v.s. actual area. Set to -1 to disable.

    minFlatness 0.5;



    //- Minimum pyramid volume. Is absolute volume of cell pyramid.

    //  Set to a sensible fraction of the smallest cell volume expected.

    //  Set to very negative number (e.g. -1E30) to disable.

    minVol 1e-13;



    //- Minimum face area. Set to <0 to disable.

    minArea -1;



    //- Minimum face twist. Set to <-1 to disable. dot product of face normal

    //- and face centre triangles normal

    minTwist 0.02;



    //- minimum normalised cell determinant

    //- 1 = hex, <= 0 = folded or flattened illegal cell

    minDeterminant 0.001;



    //- minFaceWeight (0 -> 0.5)

    minFaceWeight 0.02;



    //- minVolRatio (0 -> 1)

    minVolRatio 0.01;



    //must be >0 for Fluent compatibility

    minTriangleTwist -1;





    // Advanced



    //- Number of error distribution iterations

    nSmoothScale 4;

    //- amount to scale back displacement at error points

    errorReduction 0.75;

}





// Advanced



// Flags for optional output

// 0 : only write final meshes

// 1 : write intermediate meshes

// 2 : write volScalarField with cellLevel for postprocessing

// 4 : write current intersections as .obj files

debug 0;





// Merge tolerance. Is fraction of overall bounding box of initial mesh.

// Note: the write tolerance needs to be higher than this.

mergeTolerance 1E-6;





// ************************************************************************* //

Hi, Sorry problem solved. Turns out when I viewed my STL image and blockmesh together they did not overlap as I had taken my grid from the blockmesh file from the dynamicmesh tutorial.

Thomas

Hi everyone,

I've now managed to SHM the snake canyon tutorial stl file and another terrain file which I have created. However my problem is that each terrain created appears is very blocky (see image)s and whenever I checkMesh the output the mesh fails in 5 places.

Is solving this problem just an issue of increasing maxLocalCells and MaxGlobal cells or are there any other settings I should be adjusting as I am not confident at which settings I should adjust to fix this issue.

Can anyone suggest any solutions no matter how trivial!

Thanks,
Thom

Hi Thom.

Remember that snappyHexMesh generates new time-directories, so to check your mesh, remember to check the '3'-directory (assuming you are using a time step of unity).

If you are using foamToVTK and paraview do this:
foamToVTK -latestTime;paraview

and if things look good

cp -r ./3/polyMesh/* ./constant/polyMesh

/Mads

Thanks Mads, I've never used VTK before. Its alot quicker to load and view!

However my meshes are still of very poor quality. Will my meshes/errors improve with greater maxLocalCells and MaxGlobal cells resolution? Sorry I cannot answer these questions myself as the current machine I am using tells me it has run out of virtual memory when I increase these properties!

Thanks again guys,
Thom

Thom,

Yeah, I am using the VTK-approach all the time.
You are pointing out a major issue of the current snappyHexMesh: memory use is very high for even moderate meshe sizes. I believe it is partly due to the fact that a lot of cells are spent on areas which are cut away later on - but I have no knowledge of the implementation. SHM is a very clever tool and hopefully it will be improved even more in coming versions.

Anyway, in principle, SHM is capable of generating high quality meshes. Try to apply the following approach:

1. generate a VERY coarse mesh with blockMesh
2. use very small refinement levels in snappyHexMeshDict (say 2 2)
3. check with blockMesh;snappyHexMesh;checkmesh -latestTime

in point 3 you check your mesh count and quality and go back to item 1-2 and refine. Another possibility is to do a refineWallLayer which refines your near-STL-wall-layer. I can only do one refinement though and then the cells collapse, but you might try that. Also play with layers to get enough cells near walls.

/Mads

Hi everybody,

Mads you were quite right about the time directory. Although paraFoam was displaying time step number 2 or 3 it was still displaying the mesh from timestep 1.

I have now been able to get the snakecanyon tutorial to mesh correctly but with my own example it appears to be meshing fine in the x axis but is not creating a smooth mesh on the y (see pictures).

Can anyone speculate why this may be happening? I've been trying to sort this for a few days and any help would be most welcome!

Thanks
Thom

I want to make a 3d geometry like a tube using blockMesh utility. but it seems like there is no internal mesh in the geometry. is it possible to use the blockmesh to create 3D meshes?

Running SHM in parallel for SnakeRiverCanyon but having problems

I've been working with the snappyHexMesh version of the SnakeRiverCanyon tutorial. First of all, thanks very much to all the contributors to this thread. It has been an enormous help. I have followed Thom's original setup as closely as possible, the only difference being that I had to run SHM in parallel (4 tasks) because it ran out of memory running single process. I think I'm close, but have run into some warnings and then a segmentation fault trying to reassemble the mesh using reconstructPar.

(Before getting to this stage it was necessary to copy files from the constant directory into the processor directories using this script:

Code:

#!/bin/csh

foreach p ( processor* )

  echo $p

  foreach n ( 0 1 2 3 )

    cp $p/constant/polyMesh/pointProcAddressing $p/$n/polyMesh

    cp $p/constant/polyMesh/faceProcAddressing $p/$n/polyMesh

    cp $p/constant/polyMesh/boundaryProcAddressing $p/$n/polyMesh

    cp $p/constant/polyMesh/cellProcAddressing $p/$n/polyMesh

  end

end

)
The quoted text below is what I am seeing when I run reconstructPar. Obviously the segfault is the show-stopper, but I'm also wondering about the set of 4 warnings that precede it. Are these to be expected or might they be part of the problem too? I would very much appreciate any suggestions.

Also wondered: does anyone know of documentation or explanatory text for the SnakeRiverCanyon tutorial? There appear to be just the files themselves at the link http://repo.or.cz/w/OpenFOAM-1.6.x.git/tree/HEAD:/tutorials/mesh/moveDynamicMesh/SnakeRiverCanyon .

Thanks!

Quote:

> reconstructPar -time 3
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.7.x |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 1.7.x-4ff0ce28fdbc
Exec : reconstructPar -time 3
Date : Aug 04 2010
Time : 15:57:43
Host : rm402
PID : 12136
Case : /home/jmichal/SnakeRiverCanyon
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Time = 3

--> FOAM Warning :
From function polyMesh::readUpdateState polyMesh::readUpdate()
in file meshes/polyMesh/polyMeshIO.C at line 204
Number of patches has changed. This may have unexpected consequences. Proceed with care.
--> FOAM Warning :
From function polyMesh::readUpdateState polyMesh::readUpdate()
in file meshes/polyMesh/polyMeshIO.C at line 204
Number of patches has changed. This may have unexpected consequences. Proceed with care.
--> FOAM Warning :
From function polyMesh::readUpdateState polyMesh::readUpdate()
in file meshes/polyMesh/polyMeshIO.C at line 204
Number of patches has changed. This may have unexpected consequences. Proceed with care.
--> FOAM Warning :
From function polyMesh::readUpdateState polyMesh::readUpdate()
in file meshes/polyMesh/polyMeshIO.C at line 204
Number of patches has changed. This may have unexpected consequences. Proceed with care.
--> FOAM Warning :
From function reconstructPar
in file reconstructPar.C at line 159
readUpdate for the reconstructed mesh:0
readUpdate for the processor meshes :3
These should be equal or your addressing might be incorrect. Please check your time directories for any mesh directories.
Reconstructing FV fields

Reconstructing volScalarFields

cellLevel
#0 Foam::error:p rintStack(Foam::Ostream&) in "/home/jmichal/OpenFOAM/OpenFOAM-1.7.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::sigSegv::sigSegvHandler(int) in "/home/jmichal/OpenFOAM/OpenFOAM-1.7.x/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 ?? in "/lib64/libc.so.6"
#3 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::fvFieldReconstructor::reconstructFvVolumeFie ld<double>(Foam::IOobject const&) in "/home/jmichal/OpenFOAM/OpenFOAM-1.7.x/applications/bin/linux64GccDPOpt/reconstructPar"
#4 void Foam::fvFieldReconstructor::reconstructFvVolumeFie lds<double>(Foam::IOobjectList const&, Foam::HashSet<Foam::word, Foam::string::hash> const&) in "/home/jmichal/OpenFOAM/OpenFOAM-1.7.x/applications/bin/linux64GccDPOpt/reconstructPar"
#5 main in "/home/jmichal/OpenFOAM/OpenFOAM-1.7.x/applications/bin/linux64GccDPOpt/reconstructPar"
#6 __libc_start_main in "/lib64/libc.so.6"
#7 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/home/jmichal/OpenFOAM/OpenFOAM-1.7.x/applications/bin/linux64GccDPOpt/reconstructPar"
Segmentation fault (core dumped)

Solution to Thom's Problem

Thom,

Did you solve your problems with snappy? I had exactly the same results as you posted June 9th. It turned out that starting from the motorBike tutorial was not the right approach. The problem was system directory files (fvSolution and fvSchemes) and/or constant directory files (dynamicMeshDict and transportProperties). Try the setup in the attached example (you can run the example first like a tutorial). If you already figured this out, the solution is posted for others.

Alan

How to run snappyHexMesh in parallel?

Hi michalakes,

could you post how you managed to run snappyHexMesh in parallel?

I do the following:
1) blockMesh
2) decomposePar my case (with metis and/or simple)
3) mpirun -np 8 snappyHexMesh -parallel

Error (snippet):
----------------------------------------------
[7] --> FOAM FATAL ERROR:
[7] You have selected decomposition method decompositionMethod which is not parallel aware.
Please select one that is (hierarchical, parMetis)
[7]
[7] From function snappyHexMesh
[7] in file snappyHexMesh.C at line 341.
[7]
FOAM parallel run exiting
----------------------------------------------

Using hierarchical it starts running correctly (don't know how to use parMetis though, as there is no line in the decomposeParDict saying parMetis, and no coefficients section; and metis doesn't work either), but this crashes when I reconstruct it because of a missing file:

--> FOAM FATAL IO ERROR:
cannot open file

file: ...casedir/snappy_partest/processor0/2/polyMesh/pointProcAddressing at line 0.

From function regIOobject::readStream()
in file db/regIOobject/regIOobjectRead.C at line 61.

FOAM exiting

Edit: Just saw your script in your posting, will try that now. Anyway, the selection of methods still is not clear to me. What did you use?

2nd Edit: Get the same segmentation fault as you. Did you solve the problem?

Thank you and best regards

Pascal.

Hi Pascal,
try using reconstructParMesh utility to rebuild your mesh. I used only the hierarchical method to run snappyHexMesh in parallel, anyway if you look at the parMetisDecomp.C source code you'll see something like this:
----
// Check for user supplied weights and decomp options
if (decompositionDict_.found("metisCoeffs"))
{
const dictionary& metisCoeffs =
decompositionDict_.subDict("metisCoeffs");
....
---

Therefore this method will look for the metisCoeff subDict.

Kind regards

Norman

Once I give the moveDynamicMesh, what should I do, then??

Hi:

Just to let you know that I tried AlanR suggestion to create a mesh for a river using an stl file for the bathymetry and it worked perfect for me. Attached you will find some pictures for the mesh using SHM and AlanR case files. Just make sure that the stl file covers more area than the hex mesh, otherwise, it will not work properly.

Thanks AlanR.

Christian

Christian,

Your screenshots are exactly what I am shooting for. However, the .stl file in AlanR's example looks quite different than your example...

Also, could you elaborate a bit on the "oversized" .stl comment regarding the mesh creation?

Thanks!

Jeff

I ran across this tutorial website which may have explained why the stl file should be larger than the mesh block.

Essentially, the stl forms an internal surface inside the block where the mesh is to be refined.

If the block is larger than the stl surface, the mesh generation routine will have to "make a guess" as to what the internal characteristics of the mesh should be like, and this can lead to some unpleasant results...

Quote:

Originally Posted by AlanR (Post 272216)

Alan,

Regarding to your stl. I upload my case using moveDynamicMesh. To those who need it. Thanks for you stl.

opps, file is too large. If anyone is interested into this just Email me.