ENH: dynamicRefineFvMesh: initial commit. Not working. See README.txt

README.txt

+Notes from merging Integration-TUD-addOns and plus.
+
+0)
+- Changed
+dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+to created internal faces out-of-nothing.
+- Added the mapNewInternalFaces to dynamicRefineFvMesh
+- unset FOAM_SETNAN, FOAM_SIGFPE and run
+
+    testCases/testAMRandLoadBalancing/damBreakWithObstacle
+
+- However when writing V0 we noticed a nan in the V0 field.
+  (this doesn't get written anymore but it might indicate a bug)
+
+
+1) Tried moving additional mapping of surface fields to
+
+    virtual dynamicRefineFvMesh::mapFields(const mapPolyMesh&);
+
+so it would get called from the fvMesh::updateMesh. However
+this mapping (explicitly) gets done using unadapted addressing
+and it causes the cell addressing to be wrong:
+
+
+[1] --> FOAM FATAL ERROR:
+[1] index 826015320 out of range 0 ... 71079
+[1]
+[1]     From function void Foam::UList<T>::checkIndex(Foam::label) const [with T = Foam::cell; Foam::label = int]
+[1]     in file /home/preston2/mattijs/OpenFOAM/work/OpenFOAM-plus.integration-TUD/src/OpenFOAM/lnInclude/UListI.H at line 106.
+
+
+From dynamicRefineFvMesh::mapNewInternalFaces : (I think)
+
+                    cell faceOwner = this->cells()[owner[facei]];
+
+So this adaptation should be done in a separate pass
+
+
+2) Mapping new internal faces:
+- currently done by averaging the values of 'properly' mapped
+  faces from owner and neighbour.
+- this would not work if a cell gets split into 3x3x3
+- instead this should be done by some geometric interpolation
+  from point values?
+- have selection mechanism instead or use oriented flag on surfaceFields
+
+3) Current averaging is valid only for internal faces & only scalar/vector

src/dynamicFvMesh/Make/files

@@ -5,6 +5,7 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
 dynamicMultiMotionSolverFvMesh/dynamicMultiMotionSolverFvMesh.C
 dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
 dynamicRefineFvMesh/dynamicRefineFvMesh.C
+dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
 dynamicMotionSolverListFvMesh/dynamicMotionSolverListFvMesh.C
 
 simplifiedDynamicFvMesh/simplifiedDynamicFvMeshes.C

src/dynamicFvMesh/Make/options

 EXE_INC = \
+    -g -DFULLDEBUG -O0 -g \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
+    -ldecompositionMethods \
     -lmeshTools \
     -ldynamicMesh \
     -lfiniteVolume

src/dynamicFvMesh/dynamicRefineBalancedFvMesh/balanceParDict

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.1.1                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      balanceParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 2;
+
+method          ptscotch;
+
+constraints
+{
+    refinementHistory
+    {
+        type    refinementHistory;
+    }
+}
+// ************************************************************************* //

src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicMeshDict

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicRefineBalancedFvMesh;
+
+refinementControls
+{
+    enableRefinementControl  true;
+    
+    fields // must be scalarFields
+    (
+        //alpha (min max refineLevel)
+        alpha (0.01 0.99 2) // refine cells where alpha in [0.01:0.99] with maximal 2 refinement layers
+    );
+    
+    interface // must be a scalarField (only one dictionary!)
+    (
+        alpha // refine interface (found based on snGrad of alpha > 0.1) 
+        {
+            innerRefLayers 2; // describes how many cell layers inside phase alpha are to be refined
+            outerRefLayers 5; // describes how many cell layers outside phase alpha are to be refined
+            
+            // optional settings:
+            maxRefineLevel 4; // max refinement layers; Default: maxRefinement from dynamicRefineFvMeshCoeffs is used
+            // to get slower than 2:1 refinement; add #nAddLayers between each refinement level at the interface
+            nAddLayers 1; //Default: 0 
+        }
+    );
+    
+    gradients // must be scalars
+    (
+        // arguments as in 'fields'
+        // min/max values are based on mag(fvc::grad(volScalarField)) * cellVolume
+        T    (0.01 10 1)
+    );
+    
+    curls // must be vectors
+    (
+        // arguments as in 'fields'
+        // min/max values are based on mag(fvc::curl(volVectorField))
+        U   (0.5 1 2)
+    );
+    
+    regions
+    (
+        boxToCell
+        {
+            minLevel 1;
+            
+            box (-1 0.001 0.002)(1 0.005 0.003);
+        }
+        
+    );
+}
+
+dynamicRefineFvMeshCoeffs
+{
+    // Extra entries for balancing
+    enableBalancing true;
+    allowableImbalance 0.15;
+
+    // How often to refine
+    refineInterval  10;
+    
+    // Field to be refinement on (set it to 'internalRefinementField' to use the
+    // refinementControls dictionary entries above)
+    field           internalRefinementField;
+    
+    // Refine field inbetween lower..upper
+    lowerRefineLevel 0.5; // do not change
+    upperRefineLevel 3.5; // maxRefinement+0.5
+    
+    // If value < unrefineLevel unrefine
+    unrefineLevel   -0.5; // do not change
+    
+    // Have slower than 2:1 refinement
+    nBufferLayers   4;
+    
+    // Refine cells only up to maxRefinement levels
+    maxRefinement   3;
+    
+    // Stop refinement if maxCells reached
+    maxCells        200000;
+    
+    // Flux field and corresponding velocity field. Fluxes on changed
+    // faces get recalculated by interpolating the velocity. Use 'none'
+    // on surfaceScalarFields that do not need to be reinterpolated.
+    correctFluxes
+    (
+        (phi Urel)
+        (phiAbs U)
+        (phiAbs_0 U_0)
+        (nHatf none)
+        (rho*phi none)
+        (ghf none)
+    );
+
+    // List of non-flux surface<Type>Fields to be mapped
+    // only for new internal faces (AMR refine)
+    mapSurfaceFields
+    (
+        Uf
+        Uf_0
+    );
+    
+    // Write the refinement level as a volScalarField
+    dumpLevel       true;
+}
+
+
+// ************************************************************************* //

src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 Tyler Voskuilen
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is a derivative work of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+
+Class
+    Foam::dynamicRefineBalancedFvMesh
+
+SourceFiles
+    dynamicRefineBalancedFvMesh.C
+    dynamicRefineBalancedFvMeshTemplates.C
+
+Authors
+    T.G. Voskuilen ( https://github.com/tgvoskuilen/meshBalancing )
+    Daniel Deising <deising@mma.tu-darmstadt.de>
+    Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
+    All rights reserved.
+
+Description
+    A fvMesh with built-in multi-criterion refinement
+    and run-time load balancing.
+
+    You may refer to this software as :
+    //- full bibliographic data to be provided
+
+    This code has been developed by :
+        Daniel Rettenmaier (main developer).
+    
+    Method Development and Intellectual Property :
+        T.G. Voskuilen (Purdue University)
+        Timothée Pourpoint <timothee@purdue.edu> (Purdue University
+        Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
+        Daniel Deising <deising@mma.tu-darmstadt.de>
+        Holger Marschall <marschall@csi.tu-darmstadt.de>
+        Dieter Bothe <bothe@csi.tu-darmstadt.de>
+        Cameron Tropea <ctropea@sla.tu-darmstadt.de>
+
+        School of Aeronautics and Astronautics
+        Purdue University
+        Mathematical Modeling and Analysis
+        Institute for Fluid Mechanics and Aerodynamics
+        Center of Smart Interfaces
+        Technische Universitaet Darmstadt
+       
+    If you use this software for your scientific work or your publications,
+    please don't forget to acknowledge explicitly the use of it.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicRefineBalancedFvMesh_H
+#define dynamicRefineBalancedFvMesh_H
+
+
+#include "dynamicFvMesh.H"
+#include "dynamicRefineFvMesh.H"
+#include "hexRef8.H"
+#include "PackedBoolList.H"
+#include "Switch.H"
+#include "decompositionMethod.H"
+#include "fvMeshDistribute.H"
+#include "mapDistributePolyMesh.H"
+#include "HashTable.H"
+#include "topoSetSource.H"
+#include "cellSet.H"
+#include "PtrDictionary.H"
+#include "dictionaryEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class dynamicRefineBalancedFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicRefineBalancedFvMesh
+:
+    public dynamicRefineFvMesh
+{
+   
+private:
+
+        //- (non-flux)SurfaceFields to map
+        HashTable<word> mapSurfaceFields_;
+
+        //-
+        volScalarField* internalRefinementFieldPtr_;
+        volScalarField* targetLevelPtr_;
+        volScalarField* isLevelPtr_;
+        
+        //-
+        HashTable< List<scalar> > fields_;
+        
+        //-
+        HashTable< List<scalar> > gradFields_;
+        
+        //-
+        HashTable< List<scalar> > curlFields_;
+        
+        //-
+        HashTable< dictionary > interface_;
+        
+        //-
+        PtrList<entry> refinedRegions_; 
+        
+        //-
+        Switch enableRefinementControl_;
+        
+        //-
+        void updateRefinementField();
+        
+        //-
+        void readRefinementDict();
+        
+        //-
+        List<scalar> readRefinementPoints();
+        
+        //- Disallow default bitwise copy construct
+        dynamicRefineBalancedFvMesh(const dynamicRefineBalancedFvMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dynamicRefineBalancedFvMesh&);
+
+        //-
+        label topParentID(label p);
+
+        //-
+        label nBufferLayers_;
+
+        //-
+        bool rebalance_;
+
+        //- Map non-flux surface<Type>Fields for new internal faces
+        //  (from cell splitting)
+        template <class T>
+        void mapNewInternalFaces(const labelList& faceMap);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicRefineBalancedFvMesh");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        explicit dynamicRefineBalancedFvMesh(const IOobject& io);
+
+
+    //- Destructor
+    virtual ~dynamicRefineBalancedFvMesh();
+
+
+    // Member Functions
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+
+        //- Map all fields in time using given map
+        virtual void mapFields(const mapPolyMesh& mpm);
+
+        //- Template to update all volField boundaries
+        template<class Type> void correctBoundaries();
+
+        //-
+        bool rebalance()
+        {
+            return rebalance_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "dynamicRefineBalancedFvMeshTemplates.C"
+#endif
+
+#endif
+
+// ************************************************************************* //

src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 Tyler Voskuilen
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "volMesh.H"
+#include "fvPatchField.H"
+#include "surfaceFields.H"
+
+template <class T>
+void Foam::dynamicRefineBalancedFvMesh::mapNewInternalFaces
+(
+    const labelList& faceMap
+)
+{
+    typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
+    HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
+
+    const labelUList& owner = this->faceOwner();
+    const labelUList& neighbour = this->faceNeighbour();
+    const dimensionedScalar deltaN = 1e-8 / pow(average(this->V()), 1.0 / 3.0);
+
+    forAllIter(typename HashTable<GeoField*>, sFlds, iter)
+    {
+        GeoField& sFld = *iter();
+        if (mapSurfaceFields_.found(iter.key()))
+        {   
+            if (debug)
+            {
+                InfoInFunction << iter.key()<< endl;
+            }
+
+
+            // Create flat version of field
+            Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
+            {
+                forAll(sFld.internalField(), iFace)
+                {
+                    tsFld[iFace] = sFld.internalField()[iFace];
+                }
+                forAll(sFld.boundaryField(), iPatch)
+                {
+                    const fvsPatchField<T>& pfld = sFld.boundaryField()[iPatch];
+
+                    label start = pfld.patch().start();
+                    forAll(pfld, faceI)
+                    {
+                        tsFld[faceI+start] = pfld[faceI];
+                    }
+                }
+            }
+
+
+            // Loop over all faces
+            for (label facei = 0; facei < nInternalFaces(); facei++)
+            {
+                label oldFacei = faceMap[facei];
+
+                //- map surface field on newly generated faces
+                if (oldFacei == -1)
+                {
+                    //- find owner and neighbour cell
+                    const cell& ownFaces = this->cells()[owner[facei]];
+                    const cell& neiFaces = this->cells()[neighbour[facei]];
+
+                    //- loop over all owner/neighbour cell faces
+                    //- and find already mapped ones (master-faces):
+                    T tmpValue = pTraits<T>::zero;
+                    scalar magFld = 0;
+                    label counter = 0;
+
+                    //- simple averaging of all neighbour master-faces
+                    forAll(ownFaces, iFace)
+                    {
+                        if (faceMap[ownFaces[iFace]] != -1)
+                        {
+                            tmpValue += tsFld[ownFaces[iFace]];
+                            magFld += mag(tsFld[ownFaces[iFace]]);
+                            counter++;
+                        }
+                    }
+
+                    forAll(neiFaces, iFace)
+                    {
+                        if (faceMap[neiFaces[iFace]] != -1)
+                        {
+                            tmpValue = tmpValue + tsFld[neiFaces[iFace]];
+                            magFld += mag(tsFld[neiFaces[iFace]]);
+                            counter++;
+                        }
+                    }
+
+                    if(counter > 0)
+                    {
+                        if 
+                        (   
+                            GeometricField<T, fvsPatchField, surfaceMesh>::typeName 
+                                == "surfaceScalarField"
+                        )
+                        {
+                            tmpValue /= counter;
+                        }                            
+                        else if 
+                        (    
+                            GeometricField<T, fvsPatchField, surfaceMesh>::typeName 
+                                == "surfaceVectorField"
+                        )
+                        {
+                            magFld /= counter;
+                            tmpValue *= magFld/(mag(tmpValue)+deltaN.value());
+                        }
+                        else
+                        {
+                            FatalErrorInFunction
+                                << "mapping implementation only valid for"
+                                << " scalar and vector fields! \n Field "
+                                << sFld.name() << " is of type: "
+                                << GeometricField<T, fvsPatchField, surfaceMesh>::typeName 
+                                << abort(FatalError);
+                        }
+                    }
+
+                    sFld[facei] = tmpValue;
+                }
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::dynamicRefineBalancedFvMesh::correctBoundaries()
+{
+    typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
+
+    HashTable<GeoField*> flds(this->objectRegistry::lookupClass<GeoField>());
+
+    forAllIter(typename HashTable<GeoField*>, flds, iter)
+    {
+        GeoField& fld = *iter();
+        
+        //mimic "evaluate" but only for coupled patches (processor or cyclic)
+        // and only for blocking or nonBlocking comms (no scheduled comms)
+        if
+        (
+            Pstream::defaultCommsType == Pstream::commsTypes::blocking
+         || Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
+        )
+        {
+            label nReq = Pstream::nRequests();
+
+            forAll(fld.boundaryField(), patchi)
+            {
+                if(fld.boundaryField()[patchi].coupled())
+                {
+                    fld.boundaryFieldRef()[patchi].initEvaluate
+                    (
+                        Pstream::defaultCommsType
+                    );
+                }
+            }
+
+            // Block for any outstanding requests
+            if
+            (
+                Pstream::parRun()
+             && Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
+            )
+            {
+                Pstream::waitRequests(nReq);
+            }
+
+            forAll(fld.boundaryField(), patchi)
+            {
+                if(fld.boundaryField()[patchi].coupled())
+                {
+                    fld.boundaryFieldRef()[patchi].evaluate
+                    (
+                        Pstream::defaultCommsType
+                    );
+                }
+            }
+        }
+        else
+        {
+            //Scheduled patch updates not supported
+            FatalErrorIn
+            (
+                "dynamicRefineBalancedFvMeshTemplates::correctBoundaries"
+            )   << "Unsuported communications type "
+                << Pstream::commsTypeNames[Pstream::defaultCommsType]
+                << exit(FatalError);
+        }
+    }
+}
+
+// ************************************************************************* //