CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Questions on dynamicTopoFvMesh (https://www.cfd-online.com/Forums/openfoam-solving/111926-questions-dynamictopofvmesh.html)

danvica January 18, 2013 03:02

Questions on dynamicTopoFvMesh
 
Hi,
I'm trying to understand how to use dynamicTopoFvMesh library to simulate a free moving (6DoF) object into a fluid.

I'm following the circCylinder3d tutorial just to understand how to setup my case that will use pimpleDyMFoam and I've got lots of questions...
just a couple to start:

1. What is tetFemSolution ?

2. In dynamicMeshDict there are position constrains:

Code:

//- Constrain surface mesh-motion on a specified cylinder
    cylindricalConstraints
    {
        //- Specify options per slip patch
        sideWall
        {
            axisPoint      (0.0 0.0 0.0);
            axisVector    (0.0 0.0 1.0);
            radius          1.0;
        }
    }

    //- Specify fixedValue patches for the motionSolver
    fixedValuePatches
    {
        topWall
        {
            type          angularOscillatingDisplacement;
            amplitude    -0.0125;
            //type          oscillatingDisplacement;
            //amplitude    (0 0 -0.01);
            axis          (1 0 0);
            origin        (0 0 3);
            angle0        0.0;
            omega        0.15;
            value        uniform (0 0 0);
        }
    }

but what is then motionU file for ?

Code:

boundaryField
{
    bottomWall
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }

    sideWall
    {
        type            slip;
    }

    topWall
    {
        type            oscillatingFixedValue;
        refValue        uniform (0 0 0);
        amplitude      uniform (0 0 0.5);
        frequency      0.1;
        //type            fixedValue;
        //value          uniform (0 0 0.1);
    }

    outlet
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
}

3. The only type of motion constrain I found (in mesquiteOptions) are:
// angularOscillatingDisplacement
// angularOscillatingVelocity
// oscillatingDisplacement
// oscillatingVelocity
// surfaceDisplacement
// surfaceSlipDisplacement

Is there something similar to the following type ?
Code:

oggetto
    {
        type            sixDoFRigidBodyDisplacement;
        centreOfMass    (0 0 0.2);
        momentOfInertia (0.000001 0.000001 0.000001);
        mass            0.01;
        rhoName        rhoInf;
        rhoInf          1000;  // needed only for solvers solving for kinematic pressure
        report          on;
        value          uniform (0 0 0);
    }

I'm sorry if these question seems I little confused but I definitely need a starting point.

Thanks.

deepsterblue January 18, 2013 09:14

1 Attachment(s)
Daniele,

Sorry if the whole thing was confusing to you. I'll briefly explain:

1. The tetFemSolution file is for the tetDecompositionMotionSolvers in OpenFOAM-1.6-ext. If you don't use that library, then this file has no effect, since the mesquiteMotionSolver does not use it any way.

2. The motionU file is a way for motionSolvers to input boundary conditions, and again, is not used by the mesquiteMotionSolver, and is input using dynamicMeshDict instead. However, it appears that the fixedValuePatches approach is too simplistic, and using more complicated patch fields like sixDoFRigidBodyDisplacement may not work this way (this particular patch writes the initial position to disk, which is required for solution restarts). I've attached a diff patch which I've tested only briefly with the Port-2.1.x branch of the github repository - for some reason, it doesn't work with OF-1.6-ext. This patch basically points the mesquiteMotionSolver to use the pointDisplacement field (which the fvMotionSolver uses) to calculate 6-DoF properties.

3. Remember to add a libs entry in your controlDict that loads libforces.so (located at $FOAM_SRC/postProcessing/functionObjects/forces, if you're curious), so that the 6DoF stuff is picked up, and those patches become available for run time selection.

Hope this helps.

danvica January 18, 2013 11:09

Hi Sandeep,
many thanks for the clarifications !

I'll try the patch and let you know.

danvica January 18, 2013 12:59

It give me a strange error:

Code:

/*---------------------------------------------------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  2.1                                  |
|  \\  /    A nd          | Web:      www.OpenFOAM.org                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
/*  Windows 32 and 64 bit porting by blueCAPE: http://www.bluecape.com.pt  *\
|  Based on Windows porting (2.0.x v4) by Symscape: http://www.symscape.com  |
\*---------------------------------------------------------------------------*/
Build  : 2.1-fb1015c3f1b5
Exec  : C:\PROGRA~2\BLUECF~1.1\OpenFOAM-2.1\platforms\linuxmingw-w64DPOpt\bin\pimpleDyMFoam.exe
Date  : Jan 18 2013
Time  : 18:50:58
Host  : "UFFTECNICO7"
PID    : 5668
Case  : f:/taps/cfd/valvmov-topo
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

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

Create mesh for time = 0

Selecting dynamicFvMesh dynamicTopoFvMesh
Selecting metric Knupp
Selecting motion solver: mesquiteMotionSolver
Selecting quality metric: AspectRatioGamma
Selecting objective function: LPtoP
Selecting optimization algorithm: FeasibleNewton
Outer termination criterion (tcOuter) was not found. Using default values.
Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type laminar
Reading field rAU if present

No field sources present


PIMPLE: no residual control data found. Calculations will employ 2 corrector loops


Starting time loop

Courant Number mean: 0.001961305481 max: 0.207096507
deltaT = 0.001111111111
Time = 0.00111111


************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001149442 1.201508832 1.208379622 3.0127907370.1286772598

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001124721 1.201507976 1.208379044 3.0144733790.1286798298

 ~~~ Mesh Quality Statistics ~~~
 Min: 0.4792098079
 Max: 0.9992508882
 Mean: 0.8898675635
 Cells: 72220
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 Topo modifier time: 0.648 s
 Bisections :: Total: 0, Surface: 0
 Collapses  :: Total: 0, Surface: 0
 Swaps      :: Total: 441, Surface: 0


--> FOAM FATAL ERROR:
attempted to assign to a const reference to constant object

    From function Foam::tmp<T>::operator=(const tmp<T>&)
    in file C:/PROGRA~2/BLUECF~1.1/OpenFOAM-2.1/src/OpenFOAM/lnInclude/tmpI.H at line 270.

FOAM aborting



Backtrace:
        ZN10StackTraceC1Ev [0x6c5815de+62]
                module: C:\PROGRA~2\BLUECF~1.1\OpenFOAM-2.1\platforms\linuxmingw-w64DPOpt\lib\libstack_trace.dll
        (No symbol) [0x44ed1a0]
        RtlCaptureContext [0x775eb6f0+0]
                module: C:\Windows\system32\kernel32.dll
        (No symbol) [0x76cc10]
        (No symbol) [0x590270]
        (No symbol) [0x5904a0]
        (No symbol) [0x1b]
        (No symbol) [0x7602a8]
        RtlAllocateHeap [0x77cd3488+232]
                module: C:\Windows\SYSTEM32\ntdll.dll

This application has requested the Runtime to terminate it in an unusual way.
Please contact the application's support team for more information.

Actually I've got this error even before the patch.

Here you can find the case: https://www.box.com/s/9pg0lcuhu20th8cnccfh

Thanks for any help.

deepsterblue January 18, 2013 14:27

1 Attachment(s)
Anyway you could give me a more meaningful stack trace? I you compile with debug symbols, that would be helpful.

I think I might have seen this error before. Perhaps the attached diff would help. Not sure if this is the same problem though.

wyldckat January 18, 2013 17:07

Greetings to all!

@Sandeep: Your latest patch "topoMapperTemplates.patch" fixed the problem!

Sorry about the nearly meaningless stack trace. I'm the guy responsible for this OpenFOAM port named blueCFD, and this particular version Daniele is using is still a preview version, for which the stack trace is still very flaky.

Although it's not necessary but as future reference, on Linux (Ubuntu 12.04 x86_64), the respective stack trace was this (although not using debug symbols):
Code:

#0  Foam::error::printStack(Foam::Ostream&) in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::error::abort() in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2  void Foam::topoMapper::storeGradients<Foam::Vector<double>, Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> >(Foam::HashTable<Foam::autoPtr<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> >, Foam::word, Foam::string::hash>&) const in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#3  Foam::topoMapper::storeGradients() const in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#4  Foam::topoMapper::storeMeshInformation() const in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#5  Foam::dynamicTopoFvMesh::resetMesh() in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#6  Foam::dynamicTopoFvMesh::update() in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#7 
 in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
#8  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#9 
 in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
Aborted (core dumped)

Best regards,
Bruno

danvica January 19, 2013 03:01

1 Attachment(s)
Thanks to both Bruno and Sandeep!

But I still have problems:

I cleared all the position BC from dynamicMeshDict and I used the following pointDisplacement dict:

Code:


dimensions      [0 1 0 0 0 0 0];

internalField  uniform (0 0 0);

boundaryField
{
    wall
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
    inlet
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
    outlet
    {
        type            fixedValue;
        value          uniform (0 0 0);
    }
    oggetto
    {
        type            sixDoFRigidBodyDisplacement;
        centreOfMass    (0 0 0.2);
        momentOfInertia (0.000001 0.000001 0.000001);
        mass            0.01;
        rhoName        rhoInf;
        rhoInf          1000;  // needed only for solvers solving for kinematic pressure
        report          on;
        value          uniform (0 0 0);
    }
}

dynamicMeshDict is:

Code:


//- Select the type of dynamicFvMesh
dynamicFvMesh          dynamicTopoFvMesh;

//- Select the type of motionSolver
motionSolverLibs ("libmesquiteMotionSolver.so" "libfvMotionSolvers.so");
solver                  mesquiteMotionSolver;

mesquiteOptions
{
    //- Optimization metric
    optMetric              AspectRatioGamma;

    //- Objective function
    objFunction            LPtoP;

    //- Optimization algorithm
    optAlgorithm            FeasibleNewton;

 
    tcInner
    {
        absGradL2            1e-4;
        cpuTime              0.5;
    }

    //- Power value for the LPtoP objective function
    pValue                  2;
    power                  2;

    //- Specify a tolerance for the surface-smoothing CG solver
    tolerance              1e-2;

    //- Specify number of CG sweeps for surface-smoothing
    nSweeps                2;

    //- Specify slip patches for the motionSolver
    slipPatches
    {
     
    }

    //- Constrain surface mesh-motion on a specified cylinder
    cylindricalConstraints
    {
     
    }

    //- Specify fixedValue patches for the motionSolver
    fixedValuePatches
    {
       
       
   
   
    }
   

    //- Specify interval for surface smoothing
    surfInterval            1;
}

//- Options for dynamicTopoFvMesh
dynamicTopoFvMesh
{
   
    allOptionsMandatory no;

   
    //debug              2;

    //- Specify the number of threads
    threads            1;

    //- Specify re-meshing interval
    //- Negative value implies no re-meshing
    interval            1;

    //- Specify whether the length-scale field
    //- should be dumped to disk
    dumpLengthScale    false;

    //- sliverThreshold specifies the
    //- quality criteria for sliver removal.
    sliverThreshold    0.35;

    //- Should the tool attempt to remove slivers
    //- that fall below the sliverThreshold value?
    removeSlivers      false;

    //- Skip mapping step. Useful while using
    //- this tool as a pre-processor
    // skipMapping        true;

    // Toggle edgeRefinement on/off
    edgeRefinement      off;

 
    refinementOptions
    {
        collapseRatio  0.5;
        bisectionRatio  1.5;
        growthFactor    1.0;

        fixedLengthScalePatches
        {
            oggetto    0.2;
        }

        //- Avoid refinement on certain patches, if necessary
        noModificationPatches
        {
            inlet;
            outlet;
            wall;
            oggetto;
        }

        //- Set floating length-scale values on certain patches
        freeLengthScalePatches
        {
            inlet;
            outlet;
            wall;
        }

     
    }

    //- Tetrahedral mesh quality metric
    tetMetric          Knupp;

    //- Avoid 2-2 swapping on certain patches
    noSwapPatches
    {}
}

This setup results in a not moving object ("oggetto" in the case I posted before).

If I then modified dynamicMeshDict adding the same position BC:

Code:

fixedValuePatches
    {
        oggetto
    {
        type            sixDoFRigidBodyDisplacement;
        centreOfMass    (0 0 0.2);
        momentOfInertia (0.000001 0.000001 0.000001);
        mass            0.01;
        rhoName        rhoInf;
        rhoInf          1000;  // needed only for solvers solving for kinematic pressure
        report          on;
        value          uniform (0 0 0);
    }
    }

Then the solver displays some information regarding centre of mass / linear velocity and angular velocity but the centre of mass is fixed. See the output below:

Code:


Create time

Create mesh for time = 0

Selecting dynamicFvMesh dynamicTopoFvMesh
Selecting metric Knupp
Selecting motion solver: mesquiteMotionSolver
Selecting quality metric: AspectRatioGamma
Selecting objective function: LPtoP
Selecting optimization algorithm: FeasibleNewton
Outer termination criterion (tcOuter) was not found. Using default values.
Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type laminar
Reading field rAU if present

No field sources present


PIMPLE: no residual control data found. Calculations will employ 2 corrector loops


Starting time loop

Courant Number mean: 0.001961305481 max: 0.207096507
deltaT = 0.001
Time = 0.001

Centre of mass: (0 0 0.2)
Linear velocity: (0 0 0)
Angular velocity: (0 0 0)


************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001149442 1.201508832 1.208379622 3.0127907370.1286772598

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001124721 1.201507976 1.208379044 3.0144733790.1286798298

 ~~~ Mesh Quality Statistics ~~~
 Min: 0.4792098079
 Max: 0.9992508882
 Mean: 0.8898675635
 Cells: 72220
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 Topo modifier time: 0.441044 s
 Bisections :: Total: 0, Surface: 0
 Collapses  :: Total: 0, Surface: 0
 Swaps      :: Total: 441, Surface: 0
 Mapping time: 0.051005 s
 Reordering time: 0.135013 s
 ~~~ No flux correction ~~~

 ~~~ Mesh Quality Statistics ~~~
 Min: 0.4792098079
 Max: 0.9992508882
 Mean: 0.8888376015
 Cells: 72226
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

GAMG:  Solving for pcorr, Initial residual = 1, Final residual = 0.009784238691, No Iterations 14
GAMG:  Solving for pcorr, Initial residual = 0.1318015667, Final residual = 0.005524330453, No Iterations 1
GAMG:  Solving for pcorr, Initial residual = 0.0274078549, Final residual = 0.001671230422, No Iterations 1
time step continuity errors : sum local = 0.003519238986, global = 2.951992158e-005, cumulative = 2.951992158e-005
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.2099944886, Final residual = 0.001067176115, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.2176037194, Final residual = 0.001119524463, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.9144495281, Final residual = 0.006921714113, No Iterations 1
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.006135893841, No Iterations 16
GAMG:  Solving for p, Initial residual = 0.2245647768, Final residual = 0.001579813968, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.04407841501, Final residual = 0.0004149844197, No Iterations 2
time step continuity errors : sum local = 0.0005491589084, global = 1.44487661e-005, cumulative = 4.396868768e-005
PIMPLE: iteration 2
DILUPBiCG:  Solving for Ux, Initial residual = 0.4474193232, Final residual = 9.817080748e-007, No Iterations 64
DILUPBiCG:  Solving for Uy, Initial residual = 0.4486384426, Final residual = 7.49347033e-007, No Iterations 61
DILUPBiCG:  Solving for Uz, Initial residual = 0.6700159044, Final residual = 8.998574123e-007, No Iterations 61
GAMG:  Solving for p, Initial residual = 0.2749979373, Final residual = 0.001482151148, No Iterations 14
GAMG:  Solving for p, Initial residual = 0.1748427853, Final residual = 0.001382902354, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.03707025143, Final residual = 4.711430883e-008, No Iterations 22
time step continuity errors : sum local = 5.367153172e-007, global = -5.911960027e-008, cumulative = 4.390956808e-005
ExecutionTime = 8.398 s  ClockTime = 9 s

Courant Number mean: 1.488762178 max: 300.8061065
deltaT = 1.639344262e-005
Time = 0.00101639

Centre of mass: (0 0 0.2)
Linear velocity: (4.896733259e-005 2.424161936e-005 0.006302055408)
Angular velocity: (2.485263244e-005 -0.0001594212928 -4.137107884e-005)


************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72068 free elements of 72226 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001124721 1.203706166  1.21061955 3.0144733790.1291942728

BTW: I left the solver running for a while and even with a not-zero velocity the center of mass remained fixed. Why ?

Last question / problem:

I made the tests with:

edgeRefinement off;

I think I don't need it so far. However when I enable the refinement it returns an error (see enclosed file).


Any idea or comment ?

danvica January 19, 2013 05:23

A small update:

@Sandeep: in the first patch you sent me there was an additional flag called "usePointDisplacement".

Code:

if (optionsDict.found("usePointDisplacement"))
+    {
+        boundaryConditions_.reset
+        (
+            new pointVectorField
+            (
+                IOobject
+                (
+                    "pointDisplacement",
+                    Mesh_.time().timeName(),
+                    Mesh_,
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                pointMesh::New(Mesh_)
+            )
+        );
+    }

From what I can understand it seems that you have to enable it in order to use the pointDisplacement file.

I tried but it returns more or less the same error:

Code:

/*---------------------------------------------------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  2.1                                  |
|  \\  /    A nd          | Web:      www.OpenFOAM.org                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
/*  Windows 32 and 64 bit porting by blueCAPE: http://www.bluecape.com.pt  *\
|  Based on Windows porting (2.0.x v4) by Symscape: http://www.symscape.com  |
\*---------------------------------------------------------------------------*/
Build  : 2.1-fb1015c3f1b5
Exec  : C:\PROGRA~2\BLUECF~1.1\OpenFOAM-2.1\platforms\linuxmingw-w64DPOpt\bin\pimpleDyMFoam.exe
Date  : Jan 19 2013
Time  : 11:16:40
Host  : "UFFTECNICO7"
PID    : 5212
Case  : f:/taps/cfd/valvmov-topo
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

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

Create mesh for time = 0

Selecting dynamicFvMesh dynamicTopoFvMesh
Selecting metric Knupp
Selecting motion solver: mesquiteMotionSolver
Selecting quality metric: AspectRatioGamma
Selecting objective function: LPtoP
Selecting optimization algorithm: FeasibleNewton
Outer termination criterion (tcOuter) was not found. Using default values.
Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type laminar
Reading field rAU if present

No field sources present


PIMPLE: no residual control data found. Calculations will employ 2 corrector loops


Starting time loop

Courant Number mean: 1.961305481e-005 max: 0.00207096507
deltaT = 1.19047619e-005
Time = 1.19048e-005

Centre of mass: (0 0 0.2)
Linear velocity: (0 0 0)
Angular velocity: (0 0 0)
Solving for point motion:  Initial residual: 1 Final residual: 0.007879279805 No Iterations: 12
Solving for point motion:  Initial residual: 1 Final residual: 0.005761454167 No Iterations: 18

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.000929751 1.202198694  1.20924222 3.0127907370.1303266922

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001152746 1.202272074 1.209234895 3.0212553940.1295796677

 ~~~ Mesh Quality Statistics ~~~
 Min: 0.4784923963
 Max: 0.9992322401
 Mean: 0.8895440621
 Cells: 72220
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 Topo modifier time: 0.443 s
 Bisections :: Total: 0, Surface: 0
 Collapses  :: Total: 0, Surface: 0
 Swaps      :: Total: 472, Surface: 0
 Mapping time: 0.058 s
 Reordering time: 0.139 s
 ~~~ No flux correction ~~~

 ~~~ Mesh Quality Statistics ~~~
 Min: 0.4784923963
 Max: 0.9992322401
 Mean: 0.8884998041
 Cells: 72226
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

GAMG:  Solving for pcorr, Initial residual = 1, Final residual = 0.008812231801, No Iterations 14
GAMG:  Solving for pcorr, Initial residual = 0.131029343, Final residual = 0.00548001763, No Iterations 1
GAMG:  Solving for pcorr, Initial residual = 0.02724234412, Final residual = 0.001672274306, No Iterations 1
time step continuity errors : sum local = 4.373115074e-005, global = 4.508582388e-007, cumulative = 4.508582388e-007
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.8981178464, Final residual = 0.02090304535, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.8276717749, Final residual = 0.01968424129, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.7550246921, Final residual = 0.01660081128, No Iterations 1
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.007054332566, No Iterations 14
GAMG:  Solving for p, Initial residual = 0.1219790556, Final residual = 0.000873474047, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.02592907917, Final residual = 0.0002061318523, No Iterations 2
time step continuity errors : sum local = 5.656009101e-006, global = 2.517583376e-007, cumulative = 7.026165764e-007
PIMPLE: iteration 2
DILUPBiCG:  Solving for Ux, Initial residual = 0.2348991934, Final residual = 8.590299483e-007, No Iterations 4
DILUPBiCG:  Solving for Uy, Initial residual = 0.233874458, Final residual = 8.983213027e-008, No Iterations 5
DILUPBiCG:  Solving for Uz, Initial residual = 0.089038337, Final residual = 5.353018298e-008, No Iterations 5
GAMG:  Solving for p, Initial residual = 0.0406207147, Final residual = 0.000269779655, No Iterations 8
GAMG:  Solving for p, Initial residual = 0.06275711829, Final residual = 0.0005373520539, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.01329922478, Final residual = 5.808675466e-008, No Iterations 17
time step continuity errors : sum local = 1.18446448e-009, global = -2.537820206e-010, cumulative = 7.023627944e-007
ExecutionTime = 6.299 s  ClockTime = 6 s

Courant Number mean: 0.1408376737 max: 2.181599535
deltaT = 1.411564626e-005
Time = 2.60204e-005



Backtrace:
        ZN10StackTraceC1Ev [0x6c5815de+62]
                module: C:\PROGRA~2\BLUECF~1.1\OpenFOAM-2.1\platforms\linuxmingw-w64DPOpt\lib\libstack_trace.dll
        (No symbol) [0x4148d40]
        RtlCaptureContext [0x775eb6f0+0]
                module: C:\Windows\system32\kernel32.dll
        (No symbol) [0x22ff30]
        (No symbol) [0x68e0000]
        (No symbol) [0x3342c10]
        (No symbol) [0x5d12f0]
        (No symbol) [0x5e0080]
        (No symbol) [0x48923b0]
        (No symbol) [0x6bee310]
        (No symbol) [0x48923c0]
        (No symbol) [0x19200010000f]
        (No symbol) [0x53002b002b0033]
        (No symbol) [0x206002b002b]
        (No symbol) [0x3f1]
        (No symbol) [0x100000300]
        (No symbol) [0x7e]

Any idea ?

wyldckat January 19, 2013 14:05

Hi!

First note: I don't know how to help Daniele any further on how to set-up this case for having the loose object float away with the flow. I'm only able so far to diagnose why things were crashing.


Second note: Just to inform readers that I've already sent an email to Daniele informing about the problems that were happening with this case in particular, at least about the error itself.
In summary:
  1. The case with these latest settings did not crash on Linux.
  2. Doing a full check mesh:
    Code:

    checkMesh -allTopology -allGeometry
    revealed that there were around 70 something cells with very small or null determinant.
  3. The conclusion coincides with an older diagnoses I've done in the past: on Linux it is possible to stretch a bit more the boundaries of memory access, such as accessing a negative location in an array (e.g array[-1]), while on Windows this simply leads to a crash.
  4. Full debug build for Windows didn't provide much more information, although it might have been because I mistakenly left the "symbol strip" option turned on... either way, the current implementation of stack tracing on blueCFD is still a bit flaky :(
    edit: nope, without stripping symbols it still doesn't give a decent stack trace... going to have to revamp the system...
Best regards,
Bruno

danvica January 20, 2013 10:35

3 Attachment(s)
Hi,
first I want to thank Bruno for the support.

Then, just to complete the thread, I enclose some pictures of the mesh.

It's a simple two diameter pipe, with a sphere inside of the bigger diameter part.

In order to mesh it using just tet elements I used Netgen. Everything was fine for checkmesh (skew and not-orto) but I hadn't check the small-determinant cells.

I don't know why normally these faulty cells are not so important...or maybe they are and they are the hidden cause of lots of crashes... or maybe it's just me that I'm stick to Windows.

However the main problem is that I don't know how to refine the mesh in order to avoid them.

You can see their position in the last picture (white coloured). I neither don't know what caused them. The geometry is quite simple and I cannot understand how Netgen can create such cells in those (simple) position.


I googled a little but there's not so many information avaible on the topic.

BTW I tried:
Netgen itself: mesh resolution seems to have no correlation with the number of null-determinant cells.
GMsh, Engrid: Windows version are too instables.
MeshLab: Still studying it. Help needed
Salome: I wanted to locally refine the edges but it crash immediately.


Anyway I'll leave this thread,for now, and I'll create another in the right meshing forum to try to solve this problem. I'll come back asap.

But if you have any idea post it wherever you want !! :)

deepsterblue January 21, 2013 09:01

Daniele,

If I vaguely recall, I remember having trouble with objects falling under gravity using pimpleDyMFoam as well, and that using interDyMFoam (with a bogus alpha field of all fluid) seemed to fix the problem. Again, this is a little weird since gravity is supposed to be picked up from the boundary condition in pimpleDyMFoam, and there's some funky settings where a registered gravity field is sometimes used instead. You may want to also check if pressure is being dealt with in a consistent manner for both situations.

Your 'cells with small determinants' seems okay to me from a visual observation, and if checkMesh gives you a decent non-orthogonality, you should be good on that front. I've run meshes with far worse and still succeeded.

danvica January 21, 2013 09:47

Thanks for your support Sandeep.

I will definitely try interDyMFoam but so far the main problem seems caused by the 0-det cells.

Bruno made a deep analysis of the issue that, I want to confirm, happens just on this Windows porting.

In few words, it brings to an out-of-range array access. Linux manages better than Windows this kind of problems.

I'm still looking for a way to get rid of this cells but, as you said, normally OF can manage worse meshes without problems: I don't know if I'm going to find the right tool.

Thanks.

wyldckat January 28, 2013 17:02

Greetings to all!

I've finally managed to figure out what's wrong! And I don't have enough experience to fix this :(.
It's not a problem with cell determinant being null or very small...
It's the first patch on this thread that's broken!

The problem is this part of the patch:
Code:

@@ -3583,6 +3604,14 @@ void mesquiteMotionSolver::applyFixedValuePatches()
    // Fetch the sub-dictionary
    const dictionary& optionsDict = subDict("mesquiteOptions");
 
+    if (boundaryConditions_.valid())
+    {
+        //boundaryConditions_().boundaryField().updateCoeffs();
+        boundaryConditions_().correctBoundaryConditions();
+
+        refPoints_ += boundaryConditions_().internalField();
+    }
+    else

The call to "boundaryConditions_().correctBoundaryConditions() " leads to a crash on Linux as well and gives this stack trace:
Code:

#0  Foam::error::printStack(Foam::Ostream&) in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::sigSegv::sigHandler(int) in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2  in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::valuePointPatchField<Foam::Vector<double> >::evaluate(Foam::UPstream::commsTypes) in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#4  Foam::GeometricField<Foam::Vector<double>, Foam::pointPatchField, Foam::pointMesh>::GeometricBoundaryField::evaluate() in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libfiniteVolume.so"
#5  Foam::mesquiteMotionSolver::applyFixedValuePatches() in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libmesquiteMotionSolver.so"
#6  Foam::mesquiteMotionSolver::solve() in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libmesquiteMotionSolver.so"
#7  Foam::motionSolver::newPoints() in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicMesh.so"
#8  Foam::dynamicTopoFvMesh::update() in "/home/user/OpenFOAM/user-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#9 
 in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11 
 in "/home/user/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
Segmentation fault (core dumped)

The case ready to be tested is available here: https://www.dropbox.com/s/bs3ehwn6xq...topo_v2.tar.gz
Simply run:
Code:

pimpleDyMFoam
edit: I'm using the (almost) latest version of OpenFOAM 2.1.x. As far as I can understand, the problem is due to the "boundaryConditions_" variable not being updated along with the remaining mesh, which leads to this weird crash...

Best regards,
Bruno

deepsterblue January 28, 2013 21:42

Bruno,

Is this after a topology change? If so, could you check if the patch field is being mapped after topo-changes? If not, perhaps the pointMesh mapper is being incorrectly called (or maybe not at all).

Thanks

wyldckat January 29, 2013 04:25

Hi Sandeep,

I think there are topology changes and some mapping seems to occur... but not for the second iteration, if I understand it correctly. I've used the following debug flags to turn on some more reports on what OpenFOAM does (placed inside the file "~/.OpenFOAM/2.1.x/controlDict"):
Code:

DebugSwitches
{
    primitiveMesh      1;
    mesquiteMotionSolver 1;
    polyMesh        1;
}

The respective output from running pimpleDyMFoam:
Code:

Create time

Create mesh for time = 0

Selecting dynamicFvMesh dynamicTopoFvMesh
void polyMesh::initMesh() : initialising primitiveMesh
primitiveMesh::clearGeom() : clearing geometric data
primitiveMesh::clearAddressing() : clearing topology
primitiveMesh::reset : mesh reset to nInternalPoints:-1 nPoints:14143 nEdges:-1 nInternalFaces:140933 nFaces:147947 nCells:72220
primitiveMesh::calcCells() : calculating cells
Selecting metric Knupp
primitiveMesh::calcEdges(const bool) : calculating edges, pointEdges and optionally faceEdges
primitiveMesh::faceEdges() : calculating faceEdges
primitiveMesh::edgeFaces() : calculating edgeFaces
Selecting motion solver: mesquiteMotionSolver
Reading options for mesquiteMotionSolver
Selecting quality metric: AspectRatioGamma
Selecting objective function: LPtoP
Selecting optimization algorithm: FeasibleNewton
Outer termination criterion (tcOuter) was not found. Using default values.
Reading field p

Reading field U

Reading/calculating face flux field phi

primitiveMesh::calcFaceCentresAndAreas() : Calculating face centres and face areas
primitiveMesh::calcFaceCentresAndAreas() : Finished calculating face centres and face areas
polyMesh::globalData() const : Constructing parallelData from processor topology
This needs the patch faces to be correctly matched
primitiveMesh::calcCellCentresAndVols() : Calculating cell centres and cell volumes
primitiveMesh::calcCellCentresAndVols() : Finished calculating cell centres and cell volumes
Selecting incompressible transport model Newtonian
Selecting turbulence model type laminar
Reading field rAU if present

No field sources present


PIMPLE: no residual control data found. Calculations will employ 2 corrector loops


Starting time loop

Courant Number mean: 1.961305481e-008 max: 2.07096507e-006
deltaT = 0.001111111111
Time = 0.00111111

nTet      = 72220
nHex      = 0
nPyr      = 0
nPrism    = 0
nPolyhedra = 0
Applying fixed-value patches, if any
Centre of mass: (0 0 0.2)
Linear velocity: (0 0 0)
Angular velocity: (0 0 0)

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001149442 1.201508832 1.208379622 3.0127907370.1286772598

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72062 free elements of 72220 total elements.
  There were no inverted elements detected.
  No entities had undefined values for any computed metric.

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.001124721 1.201507976 1.208379044 3.0144733790.1286798298

 ~~~ Mesh Quality Statistics ~~~
 Min: 0.4792098079
 Max: 0.9992508882
 Mean: 0.8898675635
 Cells: 72220
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 Topo modifier time: 0.365021 s
 Bisections :: Total: 0, Surface: 0
 Collapses  :: Total: 0, Surface: 0
 Swaps      :: Total: 441, Surface: 0
primitiveMesh::calcCellCells() : calculating cellCells
 Mapping time: 0.042002 s
 Reordering time: 0.103006 s
void polyMesh::clearAddressing() : clearing topology
primitiveMesh::clearAddressing() : clearing topology
void polyMesh::setInstance(const fileName& inst) : Resetting file instance to "0.00111111"
void polyMesh::initMesh() : initialising primitiveMesh
primitiveMesh::clearGeom() : clearing geometric data
primitiveMesh::clearAddressing() : clearing topology
primitiveMesh::reset : mesh reset to nInternalPoints:-1 nPoints:14143 nEdges:-1 nInternalFaces:140945 nFaces:147959 nCells:72226
void polyMesh::clearGeom() : clearing geometric data
primitiveMesh::clearGeom() : clearing geometric data
void polyMesh::clearAddressing() : clearing topology
primitiveMesh::clearAddressing() : clearing topology
void polyMesh::updateMesh(const mapPolyMesh&) : updating addressing and (optional) pointMesh/pointFields
void polyMesh::setInstance(const fileName& inst) : Resetting file instance to "0.00111111"
Mapping volTensorField grad(U)
primitiveMesh::calcCellCentresAndVols() : Calculating cell centres and cell volumes
primitiveMesh::calcFaceCentresAndAreas() : Calculating face centres and face areas
primitiveMesh::calcFaceCentresAndAreas() : Finished calculating face centres and face areas
primitiveMesh::calcCellCentresAndVols() : Finished calculating cell centres and cell volumes
tmp<scalarField> polyMesh::movePoints(const pointField&) :  Moving points for time 0.001111111111 index 1
bool primitiveMesh::checkMeshMotion(const pointField& newPoints, const bool report) const: checking mesh motion
Min volume = 3.690509017e-010.  Total volume = 0.001987458117.  Cell volumes OK.
Min area = 8.235135209e-007.  Face areas OK.
Pyramid volumes OK.
Non-orthogonality check OK.
Mesh motion check OK.
primitiveMesh::clearGeom() : clearing geometric data
primitiveMesh::calcCellCentresAndVols() : Calculating cell centres and cell volumes
primitiveMesh::calcFaceCentresAndAreas() : Calculating face centres and face areas
primitiveMesh::calcFaceCentresAndAreas() : Finished calculating face centres and face areas
primitiveMesh::calcCellCentresAndVols() : Finished calculating cell centres and cell volumes
 ~~~ No flux correction ~~~
tmp<scalarField> polyMesh::movePoints(const pointField&) :  Moving points for time 0.001111111111 index 1
bool primitiveMesh::checkMeshMotion(const pointField& newPoints, const bool report) const: checking mesh motion
Min volume = 3.698076689e-010.  Total volume = 0.001987458117.  Cell volumes OK.
Min area = 8.244142416e-007.  Face areas OK.
Pyramid volumes OK.
Non-orthogonality check OK.
Mesh motion check OK.
primitiveMesh::clearGeom() : clearing geometric data
primitiveMesh::calcCellCentresAndVols() : Calculating cell centres and cell volumes
primitiveMesh::calcFaceCentresAndAreas() : Calculating face centres and face areas
primitiveMesh::calcFaceCentresAndAreas() : Finished calculating face centres and face areas
primitiveMesh::calcCellCentresAndVols() : Finished calculating cell centres and cell volumes
Clearing out mesquiteMotionSolver for topo-changes

 ~~~ Mesh Quality Statistics ~~~
 Min: 0.4792098079
 Max: 0.9992508882
 Mean: 0.8888376015
 Cells: 72226
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

polyMesh::globalData() const : Constructing parallelData from processor topology
This needs the patch faces to be correctly matched
GAMG:  Solving for pcorr, Initial residual = 1, Final residual = 0.009784238691, No Iterations 14
GAMG:  Solving for pcorr, Initial residual = 0.1318015667, Final residual = 0.005524330453, No Iterations 1
GAMG:  Solving for pcorr, Initial residual = 0.0274078549, Final residual = 0.001671230422, No Iterations 1
time step continuity errors : sum local = 3.91026554e-008, global = 3.279991287e-010, cumulative = 3.279991287e-010
PIMPLE: iteration 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.001776048218, Final residual = 9.263663799e-011, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.001404779735, Final residual = 7.728408671e-011, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.00987428103, Final residual = 8.50084699e-010, No Iterations 1
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.009588103798, No Iterations 14
GAMG:  Solving for p, Initial residual = 0.1302883199, Final residual = 0.0009036470934, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.02737724854, Final residual = 0.0002262485978, No Iterations 2
time step continuity errors : sum local = 5.294873174e-009, global = 1.106183578e-010, cumulative = 4.386174865e-010
PIMPLE: iteration 2
DILUPBiCG:  Solving for Ux, Initial residual = 0.1834820032, Final residual = 5.636594204e-007, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.1830492352, Final residual = 5.109273914e-007, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.05322444189, Final residual = 1.831625782e-007, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.01625682761, Final residual = 0.0001395819713, No Iterations 10
GAMG:  Solving for p, Initial residual = 0.07015457008, Final residual = 0.0004928326396, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.01514164752, Final residual = 7.103417641e-008, No Iterations 14
time step continuity errors : sum local = 1.648956177e-012, global = 4.143241703e-013, cumulative = 4.390318106e-010
ExecutionTime = 4.949 s  ClockTime = 5 s

Courant Number mean: 0.02715254901 max: 0.4297419588
deltaT = 0.00126984127
Time = 0.00238095

primitiveMesh::calcCells() : calculating cells
nTet      = 72226
nHex      = 0
nPyr      = 0
nPrism    = 0
nPolyhedra = 0
Applying fixed-value patches, if any

Best regards,
Bruno

wyldckat February 3, 2013 13:53

1 Attachment(s)
Greetings to all!

@Sandeep: any chance you can (help us) figure out what's missing here?

Because after roughly 5h of studying dynamic meshing, studying the "mesquiteMotionSolver" and doing several trial-and-errors... I haven't managed to fix this problem.

Attached is a patch with the code where I stopped trying to fix this.

In essence, my findings are as follows:
  • "mesquiteMotionSolver" doesn't do things the same way that OpenFOAM 2.1.x does things. From what I could understand:
    • 2.1.x relies on already existing code for moving/mapping points and fields.
    • "mesquiteMotionSolver" does things directly.
  • Since "mesquiteMotionSolver" does things directly, "boundaryConditions_()" is just junk when "mesquiteMotionSolver::update()" is completed (well, somewhere in the middle...).
  • Given this situation, the only hack I can figure out that might work is to backup the last displacement on the boundaries and reconstruct the field when necessary.
  • Attached is an attempt at this, but is fails miserably, probably because using the method "boundaryConditions_().boundaryField().set(" is a very bad idea...
  • The only other way I can imagine is to write out the boundaries onto file and then recreate them, the same way that the fixed patches are made.
Best regards,
Bruno

deepsterblue February 3, 2013 23:02

I think I might know why things are not being mapped. The boundaryConditions_ member is instantiated using a pointMesh::New function, which may not be registering itself correctly for all subsequent pointField mapping.

I say this because the only field that is currently being mapped (from your logs) is the volTensorField grad(U). It would be helpful to ensure that the pointMesh exists during the mapping process (in MapGeometricFields.H), so that all pointFields associated with it are mapped accordingly as well.

Another possibility is to have the boundaryConditions_ pointVectorField mapped during the mesquiteMotionSolver::update(const mapPolyMesh& mpm) call. I think a pointMapper will need to be instantiated for this, but I can't remember off the top of my head - perhaps you could take a cue from other parts of OpenFOAM, like the fvMotionSolver::updateMesh() functions.

wyldckat February 4, 2013 05:19

Hi Sandeep,

Many thanks! I'll look into this further.

Best regards,
Bruno

wyldckat February 18, 2013 18:09

1 Attachment(s)
Hi Sandeep,

Unfortunately, I'm unable to figure out how to do the point mapping. All of my tests seem to indicate that "boundaryConditions_()" no longer has enough information when it is working inside "void mesquiteMotionSolver::update(const mapPolyMesh& mpm)".

Attached is a hack that seems to work, although it has to resort to the file system for storing the "pointDisplacement" field. It's applied on top of your two patches, if I'm not mistaken.
It's a (very) sub-optimal solution, but it gets the job done... or at least I hope so.

edit: after several iterations, it looks like the sphere is starting to get deformed... I can only guess that it's because there are points that are being erased/added, which are not mapped out with this hack... :(

Best regards,
Bruno

deepsterblue February 19, 2013 09:25

Bruno,

I'm willing to bet that the points on the sphere correspond to un-mapped points. Can you add a noModificationPatches / noSwapPatches entry in dynamicMeshDict with your "sphere" patch, and see if the deformations still happen?

If you look at the Foam:: polyMesh::updateMesh(const mapPolyMesh& mpm) function in polyMeshUpdate.C, you'll see this bit:

Code:

    if (thisDb().foundObject<pointMesh>(pointMesh::typeName))
    {
        const_cast<pointMesh&>
        (
            thisDb().lookupObject<pointMesh>
            (
                pointMesh::typeName
            )
        ).updateMesh(mpm);
    }

As Mattijs put it, this is indeed a bit of a hack, and a proper callback mechanism needs to be put into place for MeshObject (like the way it's done in 1.6-ext). You may need to add this to the mesquiteMotionSolver::updateMesh member function as well, and see if the boundary conditions get mapped.

wyldckat February 19, 2013 10:21

Hi Sandeep,

Many thanks for the feedback! I should have posted the several tests I've made so far, because I did try several different ways of mapping the "pointDisplacement" field and respective mesh.

That one in particular I think didn't work, although it's possible that I didn't manage to define the proper usage of each object and respective methods.

Later today I'll post the attempts I've made at mapping the field.

In the mean time, the problem is also hat I did further (re)searching on this topic, which lead me to thinking that I would need several more complex hacks that I found in OpenFOAM's source code, specifically the ones in:
  • "fvMeshDistribute.C":
    Code:

    00544    // Store boundary fields (we only do this for surfaceFields)
    00545    PtrList<FieldField<fvsPatchField, scalar> > sFlds;
    00546    saveBoundaryFields<scalar, surfaceMesh>(sFlds);
    ...
    00555
    00556    // Change the mesh (no inflation). Note: parallel comms allowed.
    00557    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
    00558
    00559    // Update fields. No inflation, parallel sync.
    00560    mesh_.updateMesh(map);
    00561
    00562    // Map patch fields using stored boundary fields. Note: assumes order
    00563    // of fields has not changed in object registry!
    00564    mapBoundaryFields<scalar, surfaceMesh>(map, sFlds);
    ...

    More specifically, I would need to borrow code from "fvMeshDistributeTemplates.C" and "fvMeshDistribute.C", for specifically mapping the point mesh objects, instead of volume and surface meshes.
  • Then there is:
    Code:

    void Foam::displacementFvMotionSolver::updateMesh(const mapPolyMesh& mpm)
    But seems even more complex and only details about point movement...
  • And I know there is a similar hack to the one you've showed, but made for surface meshes... but I can't find it.
Like I've said, later today I'll post the attempts I've made so far.


Best regards,
Bruno

wyldckat February 19, 2013 17:30

2 Attachment(s)
Hi Sandeep,

OK, I've got mixed results...

First for the partially-good news:
While using the hack to use the file system and by adding the sphere (named "oggetto") to the "noModificationPatches" and "noSwapPatches", it does allow to keep the sphere's surface in perfect condition while it moves.
The downside is that it crashes when the cells in the way start to get too tight and eventually collapse into negative volumes.

The error message:
Code:

Courant Number mean: 0.04541765902 max: 4.191785362
deltaT = 3.333333333e-05
Time = 0.0202

sixDoFRigidBodyMotion constraints converged in 1 iterations
Constraint force: (0.6537136637 0.02252467437 0)
Constraint moment: (0 0 0)
Centre of mass: (-1.842749434e-08 1.08451596e-08 0.2001177404)
Linear velocity: (-0.0004352669538 8.108849543e-05 0.04280597798)
Angular velocity: (1.573833728 -1.569414796 -0.03449086776)

************** QualityAssessor(free only) Summary **************

  Evaluating quality for 72216 free elements of 72375 total elements.
  THERE ARE 1 INVERTED ELEMENTS.
  (Elements invalid at 1 of 72216 sample locations.)

  1 OF 72216 ENTITIES EVALUATED TO AN UNDEFINED VALUE FOR AspectRatioGamma

          metric    minimum    average        rms    maximum    std.dev.
AspectRatioGamma 1.002851198  15.115086 3721.202667    1000000 3721.171969

 Warning: Minimum cell quality is: -0.07980588317 at cell: 68575

 ~~~ Mesh Quality Statistics ~~~
 Min: -0.07980588317
 Max: 0.9981037065
 Mean: 0.8637096283
 Cells: 72375
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



--> FOAM FATAL ERROR:
Encountered negative cell-quality!
Edge: 89295: (3253 3261)
vertexHull: 5(227 9861 10215 9856 3297)
Minimum Quality: -0.07980588317

    From function scalar dynamicTopoFvMesh::computeMinQuality(const label eIndex, labelList& hullVertices) const
    in file edgeSwap.C at line 2246.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::error::abort() in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2  Foam::dynamicTopoFvMesh::computeMinQuality(int, Foam::List<int>&) const in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#3  Foam::dynamicTopoFvMesh::swap3DEdges(void*) in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#4  Foam::dynamicTopoFvMesh::threadedTopoModifier() in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#5  Foam::dynamicTopoFvMesh::update() in "/home/bmss/OpenFOAM/bmss-2.1.x/platforms/linux64GccDPOpt/lib/libdynamicTopoFvMesh.so"
#6 
 in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8 
 in "/home/bmss/OpenFOAM/OpenFOAM-2.1.x/platforms/linux64GccDPOpt/bin/pimpleDyMFoam"
Aborted (core dumped)

The "dynamicMeshDict":
Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM Extend Project: Open Source CFD        |
|  \\    /  O peration    | Version:  1.6-ext                              |
|  \\  /    A nd          | Web:      www.extend-project.de                |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      dictionary;
    location    "constant";
    object      dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

//- Select the type of dynamicFvMesh
dynamicFvMesh          dynamicTopoFvMesh;

//- Select the type of motionSolver
motionSolverLibs ("libmesquiteMotionSolver.so" "libfvMotionSolvers.so");
solver                  mesquiteMotionSolver;

mesquiteOptions
{
    usePointDisplacement yes;
 
    //- Optimization metric
    optMetric              AspectRatioGamma;

    //- Objective function
    objFunction            LPtoP;

    //- Optimization algorithm
    optAlgorithm            FeasibleNewton;

    //- Termination criteria sub-dictionaries
    //- (takes default values if not specified)
    //-  Specifying an empty sub-dictionary
    //-  terminates with available options
    tcInner
    {
        absGradL2            1e-4;
        cpuTime              0.5;
        // lele: iterationLimit 10;
    }

    // tcOuter
    // {}

    //- For composite functions, two objectives need to be specified
    // firstFunction          LPtoP;
    // secondFunction          LInf;

    //- For scaled functions, scale and objective needs to be specified
    // scaleFunction        PMeanP;
    // scale                1.5;

    //- Power value for the LPtoP objective function
    pValue                  2;
    power                  2;

    //- Specify a tolerance for the surface-smoothing CG solver
    tolerance              1e-2;

    //- Specify number of CG sweeps for surface-smoothing
    nSweeps                2;

    //- Specify slip patches for the motionSolver
    slipPatches
    {
        //wall;
        //oggetto;
    }

    //- Constrain surface mesh-motion on a specified cylinder
    cylindricalConstraints
    {
        //- Specify options per slip patch
    //    sideWall
    //    {
    //        axisPoint      (0.0 0.0 0.0);
    //        axisVector    (0.0 0.0 1.0);
    //        radius          1.0;
    //    }
    }

    //- Specify fixedValue patches for the motionSolver
    fixedValuePatches
    {
      oggetto
      {
          type            sixDoFRigidBodyDisplacement;
          centreOfMass    (0 0 0.2);
          momentOfInertia (0.000001 0.000001 0.000001);
          mass            0.01;
          rhoName        rhoInf;
          rhoInf          1000;  // needed only for solvers solving for kinematic pressure
          report          on;
          value          uniform (0 0 0);
          constraints
          {
            maxIterations 500000000;
            fixedLine1
            {
              sixDoFRigidBodyMotionConstraint fixedLine;
              tolerance 1e-6;
              relaxationFactor 0.7;
              fixedLineCoeffs
              {
                refPoint (0 0 0.2);
                direction (0 0 1);
              }
            }
          }
      }
    }
   
    //Other fixedValuePatch types can be found here:
    //src/dynamicMesh/meshMotion/fvMotionSolver/pointPatchFields/derived
   
    // angularOscillatingDisplacement
    // angularOscillatingVelocity
    // oscillatingDisplacement
    // oscillatingVelocity
    // surfaceDisplacement
    // surfaceSlipDisplacement

    //- Specify interval for surface smoothing
    surfInterval            1;
}

//- Options for dynamicTopoFvMesh
dynamicTopoFvMesh
{
    //- Should all options be made mandatory?
    //- Useful for first-time use.
    allOptionsMandatory no;

    //- Set run-time debug level [0-5]
    //debug              2;

    //- Specify the number of threads
    threads            1;

    //- Specify re-meshing interval
    //- Negative value implies no re-meshing
    interval            1;

    //- Specify whether the length-scale field
    //- should be dumped to disk
    dumpLengthScale    false;

    //- sliverThreshold specifies the
    //- quality criteria for sliver removal.
    sliverThreshold    0.35;

    //- Should the tool attempt to remove slivers
    //- that fall below the sliverThreshold value?
    removeSlivers      true;

    //- Skip mapping step. Useful while using
    //- this tool as a pre-processor
    // skipMapping        true;

    // Toggle edgeRefinement on/off
    edgeRefinement      no;

    //- Options for edge-bisection/collapse.
    //-  The edgeRefinement flag must be set for
    //-  the following options to have effect
    refinementOptions
    {
        collapseRatio  0.5;
        bisectionRatio  1.5;
        growthFactor    1.0;

        //- By default, existing boundary edge-lengths
        //- are used for length-scales.
        //- Length-scale can be fixed for certain patches.
        fixedLengthScalePatches
        {
            oggetto    0.2;
        }

        //- Avoid refinement on certain patches, if necessary
        noModificationPatches
        {
            inlet;
            outlet;
            wall;
            oggetto;
        }

        FreeSurfacePatch
        {
            oggetto;
        }
       
       
        //- Set floating length-scale values on certain patches
        freeLengthScalePatches
        {
            inlet;
            outlet;
            wall;
        }

        //- Limit lengthScales to specified values, if necessary
        // minLengthScale  0.1;
        // maxLengthScale  0.3;

        //- Field-based refinement options
        // fieldRefinement  gamma;
        // fieldLengthScale 0.005;
        // lowerRefineLevel 0.001;
        // upperRefineLevel 0.999;
        // maxRefineLevel  4;
        // meanScale        0.015;
    }

    //- If the number of modifications are to be limited, set this option
    // maxModifications  1000;

    //- Load custom libraries for metrics
    // tetMetricLibs      ("libtetMetrics.so");

    //- Tetrahedral mesh quality metric
    tetMetric          Knupp;

    //- Avoid 2-2 swapping on certain patches
    noSwapPatches
    {
      oggetto;
    }
}

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

Attached are the two images that show the initial position and the final position nearly before the crash, due to a collapsed cell.
The bad news:
I still haven't been able to make the proper update to work. I've tried to use the following code, but without success (this is the patch for the code):
Code:

----------------- mesquiteMotionSolver/mesquiteMotionSolver.C -----------------
index 468be60..a482416 100644
@@ -4738,6 +4738,18 @@ void mesquiteMotionSolver::update(const mapPolyMesh& mpm)
        Info<< "Clearing out mesquiteMotionSolver for topo-changes" << endl;
    }
 
+    // pointMesh
+    if (Mesh_.thisDb().foundObject<pointMesh>(pointMesh::typeName))
+    {
+        const_cast<pointMesh&>
+        (
+            Mesh_.thisDb().lookupObject<pointMesh>
+            (
+                pointMesh::typeName
+            )
+        ).updateMesh(mpm);
+    }
+
    motionSolver::updateMesh(mpm);
 
    if (surfaceSmoothing_)

The crash occurs just the same, namely whenever the "boundaryConditions_()" are accessed after the mesh has been changed.

Other attempts I did in the past were, both in "mesquiteMotionSolver::update(...)":
  • Force points to move:
    Code:

        const_cast<pointMesh&>(boundaryConditions_().mesh())
            .movePoints(Mesh_.points());

  • Force an update of the mesh:
    Code:

        const_cast<pointMesh&>(boundaryConditions_().mesh())
            .updateMesh(mpm);

Both crashed while trying to make the move, presumably due to a non-existing mesh to be moved...

Best regards,
Bruno

deepsterblue February 20, 2013 09:36

Bruno,

You could have the smoother work a little harder at each time step, and see if the simulation goes further. Specify a higher cpuTime value for tcInner and see if that helps. Also, you have edgeRefinement set to 'no', which means that cells cannot be collapsed, and that probably explains why you see negative volumes.

You can play around with the growth factor to coarsen the mesh away from the 'oggetto' patch (not too agressive, otherwise you'll get too many collapses, which you can restrict using the maxModifications entry).

I'll look into the mapping when I have some time on my hands.

Sandeep

wyldckat February 20, 2013 17:53

Hi Sandeep,

Many thanks for the hints!

Unfortunately, by removing the cells that were in the way, the saved "pointDisplacement" of the internal field has a different number of cells, which leads the solver stopping to complain about the different array sizes.

Nonetheless, this does give me (again) the idea for adopting a hack similar to the strategy used for "fixedValuePatches"...

Best regards,
Bruno

deepsterblue February 21, 2013 10:16

Bruno,

Can you give me a test case that I can use? Also, you'll have to list out the details of how to reproduce the problem.

Sandeep

wyldckat February 21, 2013 12:23

Hi Sandeep,

Many thanks! Although I'll only be able to provide a clean/smaller case and more explicit instructions in about 5 hours from now :(


In the mean time, the case is the one mentioned in post #13, namely: https://www.dropbox.com/s/bs3ehwn6xq...topo_v2.tar.gz

Unpack and simply run pimpleDyMFoam.
It should crash as soon as "boundaryConditions_().correctBoundaryConditions() " is called the second time around, if I'm not mistaken.
moveDynamicMesh is not enough, due to the need for forces to work, which depend on "U" and "p"...


This is assuming that the whole "dynamicTopoFvMesh" code only has the first two source code patches you provided on this thread, for using with OpenFOAM 2.1.x.

Unless you prefer the hack I tried to use, relying on the file system, for which the "pointDisplacement" file was placed in the "constant" folder; the source code patch is on post #19.

Best regards,
Bruno

deepsterblue February 21, 2013 16:52

Bruno,

I've committed a fix and your case seems to be running with point fields mapped. I haven't tried extensively with edgeRefinement on, so you might want to test it out.

Unfortunately, since your length scales vary quite drastically across the 'wall' patch, the algorithm has a hard time trying to figure out what the optimal length scale is, and so you tend to get excessive bisections / collapses. You have a length scale for 'oggetto' specified as 0.2, while the paraview ruler tells me it's 0.002 - so you may want to check that.

Until I have a more general solution, one workaround is to limit the number at each time-step using the maxModifications entry, or perhaps you could split up the patch into several smaller ones with a uniform length scale for each.

Thanks for the bug-report.

wyldckat February 21, 2013 18:47

Hi Sandeep,

You're very welcome about the bug report! And many thanks for the fix and for the suggestions!

And no wonder I couldn't figure it out... I was only looking at the source code for "mesquiteMotionSolver" :rolleyes:

Best regards,
Bruno

joegi.geo February 26, 2013 11:38

patch motion with dynamicTopoFvMesh
 
Hi Sandeep,

I wonder if there is a difference between the implementation of patch motion between dynamicTopoFvMesh and dynamicMotionSolverFvMesh?. I am trying to do a simple oscillating body simulation by using both libraries and the output of the patch motion is totally different (about an order of magnitude in amplitude).

Using dynamicMotionSolverFvMesh i set the moving patch in pointDisplacement as follows:

Code:

    wing
    {
        type            oscillatingDisplacement;
        value          uniform ( 0 0 0 );
        amplitude      ( 0 0.2 0 );
        omega          6.28318;
    }

which gives me the desired output.

Using dynamicTopoFvMesh/mesquite I set the moving patch in dynamicMeshDict exactly in the same way, by just adding

Code:

    wing
    {
        type            oscillatingDisplacement;
        value          uniform ( 0 0 0 );
        amplitude      ( 0 0.2 0 );
        omega          6.28318;
    }

to dynamicMeshDict in the mesquiteOptions entry,and the results are totally different. I also tried by using the pointDisplacement dictionary and the output didn't change.

Any clue?,

jg

deepsterblue February 26, 2013 11:48

I see that both your entries are the same, but how different is "totally different"? What was your output? If you post your entire dynamicMeshDict file, that would be helpful.

joegi.geo February 26, 2013 12:06

patch motion with dynamicTopoFvMesh
 
Hi,

Ok, this is the entry where i set the patch motion


Code:

solver                  mesquiteMotionSolver;

dynamicFvMesh dynamicTopoFvMesh;

dynamicFvMeshLibs        ("libdynamicTopoFvMesh.so");

motionSolverLibs  ("libfvMotionSolvers.so"  "libmesquiteMotionSolver.so");

mesquiteOptions
{
    // Optimization metric
    optMetric              AspectRatioGamma;

    // Objective function
    objFunction            LPtoP;

    // Optimization algorithm
    optAlgorithm            FeasibleNewton;

    // Termination criteria sub-dictionary (takes default values if not specified)
    // Specifying an empty sub-dictionary terminates with available options
    tcInner
    {
      absGradL2            1e-4;
      cpuTime              1.25;
    }
    // tcOuter
    // {}

    // For composite functions, two objectives need to be specified
    firstFunction          LPtoP;
    secondFunction          LInf;

    // For scaled functions, scale and objective needs to be specified
    // scaleFunction        PMeanP;
    // scale                1.5;

    // Power value for the LPtoP objective function
    pValue                  2;
    power                  2;

    // Specify a tolerance for the CG solver
    tolerance              1e-2;

    // Specify number of CG sweeps
    nSweeps                1;

    // Specify slip patches for the motionSolver
    slipPatches
    {

    }

    // Set run-time debug level
    debug                  0;

    // Specify interval for surface smoothing
    surfInterval            1;



    //- Specify fixedValue patches for the motionSolver
    fixedValuePatches
    {

        ball
        {
        type            oscillatingDisplacement;
            amplitude    (0.1 0 0);
            omega        6.28316;
            value        uniform (0 0 0);
    }

    }

}

If dynamicTopoFvMesh uses the same patch motion class as dynamicMotionSolverFvMesh (which I think it does), I should get the same result for both dynamic meshing libraries, am i right?. In my case the output of both classes are totally different.

Even if I set the entry usePointDisplacement yes, and use the pointDisplacement dictionary (same one as the one used with dynamicMotionSolverFvMesh), the patch motion of both libraries are different.

deepsterblue February 26, 2013 12:12

And what's the output?

joegi.geo February 26, 2013 12:29

The simulation runs fine (I am using moveDynamicMesh), the smoothing works fine, the remeshing works fine. But, nevertheless the fact that I am using the same patch motion for both dynamic meshing classes (libdynamicTopoFvMesh and dynamicMotionSolverFvMesh), I think the motion I am getting with libdynamicTopoFvMesh is not the right one, the amplitude is like one order of magnitude higher compare to the output of dynamicMotionSolverFvMesh.

deepsterblue February 26, 2013 12:43

Ah... I know what the problem is:

It appears that the oscillatingDisplacement patchField provides an absolute value of displacement (from the mean), while the applyFixedValuePatches() function in mesquiteMotionSolver expects an incremental displacement vector for each time-step.

You see the acceleration because the displacement is continuously updated as:

refPoints_ += dPointField; in mesquiteMotionSolver.C

You could write a new boundary condition similar to oscillatingDisplacement which provides displacement increments, or alternatively, I also provide a timeVaryingDisplacement patchField in the repository that uses an interpolation table from file, which might be easier. I would have to look into a fix for the absolute vs. incremental confusion.

Can you try changing the "+=" operator in applyFixedValuePatches to simply "=" and see if that makes a difference?

joegi.geo February 26, 2013 13:12

It crash immediately.

Let's see if I can write the bc, in anycase I am not in a hurry, so probably you will come with a solution first.

Does this affect the behavior of the 6DOF solver? it is safe to use?


Have a nice day,

joel

deepsterblue February 26, 2013 13:14

Can you post a test case with dynamicMeshDict entries for both solver types?

joegi.geo February 27, 2013 03:26

1 Attachment(s)
Hi Sandeep,

You can use the ballTranslation tutorial from the last openfoam workshop. Use the attached files to run the case using dynamicMotionSolverFvMesh.
I tested everything and runs fine.

Running this two cases you clearly see the difference in the patch motion between the two solvers.

Have a nice day,

joel

deepsterblue February 27, 2013 10:54

joel,

Can you try the angularOscillating variant with zero angular velocity and see if there's a difference between both solvers?

Thanks

joegi.geo March 1, 2013 05:44

Hi Sandeep,

Same problem with angularOscillatingDisplacement.


Joel

joegi.geo March 13, 2013 06:24

Hi Sandeep,

Any update regarding this issue?


All times are GMT -4. The time now is 11:07.