New sixDoFRigidBody BC for laplaceFaceDecomposition, performance issues
Dear Foamers,
I managed to build a new sixDoFRigidBodyVelocity boundary type based on the original sixDoFRigidBodyDisplacement, in order to make the 6DoFs solver in 1.6extension work together with laplaceFaceDecomposition. Attached is the source code, just unpack and replace the original make files in: $FOAM_SRC/postProcessing/functionObjects/forces And in the time folder you'll need a motionU file instead of pointDisplacement. Following is a sample content of this file: /** C++ **\  =========    \\ / F ield  OpenFOAM Extend Project: Open Source CFD   \\ / O peration  Version: 1.6ext   \\ / A nd  Web: www.extendproject.de   \\/ M anipulation   \**/ FoamFile { version 2.0; format ascii; class tetPointVectorField; object motionU; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { stationaryWalls { type fixedValue; value uniform (0 0 0); } atmosphere { type fixedValue; value uniform (0 0 0); } floatingObject { type sixDoFRigidBodyVelocity; centreOfMass (0.5 0.5 0.5); momentOfInertia (0.08622222 0.08622222 0.144); mass 9.6; rhoInf 1; report on; value uniform (0 0 0); } It is working perfectly according to my tests, but there's still a little glitch somewhere in the code though. While performing setFields, you'll get a warning message but the field can be set without problem: /**\  =========    \\ / F ield  OpenFOAM Extend Project: Open source CFD   \\ / O peration  Version: 1.6ext   \\ / A nd  Web: www.extendproject.de   \\/ M anipulation   \**/ Build : 1.6ext Exec : setFields time 0 Date : Mar 27 2013 Time : 15:14:25 Host : ubuntu PID : 43873 Case : /home/yasquall/OpenFOAM/OpenFOAM1.6ext/tutorials/multiphase/interDyMFoam/ras/floatingObject nProcs : 1 SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Create time > FOAM Warning : From function dlLibraryTable::open(const fileName& functionLibName) in file db/dlLibraryTable/dlLibraryTable.C at line 86 could not load /home/yasquall/OpenFOAM/OpenFOAM1.6ext/lib/linux64GccDPOpt/libforces.so: undefined symbol: _ZN4Foam25FixedValuePointPatchFieldINS_17tetPolyPa tchFieldENS_12tetPointMeshENS_22tetPolyPatchFaceDe compENS_12tetFemMatrixENS_6VectorIdEEE8typeNameE Create mesh for time = 0 Reading setFieldsDict Setting field default values Setting volScalarField alpha1 Setting field region values Adding cells with center within box (100 100 100) (100 100 0.5368) Setting volScalarField alpha1 Adding cells with center within box (0.5 0.6 100) (100 100 0.8) Setting volScalarField alpha1 End Would be much appreciate if someone could help to have of look of the code attached and get rid of the bug. Many thanks in advance and good luck in Foaming! Xin L. 
sorry, forgot the attachement, here it is.
1 Attachment(s)
Here's the source pack:
Attachment 20180 After unpack and replace compile the code as follows: cd $FOAM_SRC/postProcessing/functionObjects/forces wmake libso and you're done. 
please use this one instead...

Hi Xin,
The warning in setFields is connected with the library/libraries, which you are adding in the controlDict. The warning tells you that setFields do not know how to handle one of these libraries, because setFields does not know the tetFields in OF16ext. However, since you are not setting the values in the tetFields, but merely in the volField's, then I would not be concerned with the warning as long as the simulations are running correctly. setField merely reads all the libraries defined in the controlDict, though, for testing you could remove that line in controlDict, when you execute setFields and the add again afterward. This will most probably give the same results, though without the warning:o) Good luck, Niels 
Quote:
Many thanks for the reply! Actually besides this warning message, what concerns me more is the performance of the laplaseFaceDecomposition method compared to the traditional cell centered methods. Followings are the outputs which were extracted from two different runs using the aforementioned two types of motion solvers on a same mesh (5m cells approx, on 64 CPUs): laplaseFaceDecomposition solver: Courant Number mean: 0.00494884363762 max: 0.363069639485 velocity magnitude: 8.8242247445 deltaT = 0.0226342437319 Time = 0.0862012 Correct mesh motion diffusion field. DICPCG: Solving for motionUx, Initial residual = 0.119538225675, Final residual = 7.22629688705e06, No Iterations 1000 DICPCG: Solving for motionUy, Initial residual = 0.172646722126, Final residual = 1.21688108815e05, No Iterations 1000 DICPCG: Solving for motionUz, Initial residual = 0.0634080697595, Final residual = 5.68417442456e06, No Iterations 1000 Constraint fixedLine1 error (0.000352425614415 0 4.70215607589e05) force (17661971.327 0 2356507.08661) moment (0 0 0) not converged Constraint fixedAxis1 angle 2.10734242554e08 force (0 0 0) moment (38738.9313677 630311.39261 0) converged Constraint fixedLine1 error (0.000140970245766 0 1.88086243043e05) force (7064788.53079 0 942602.834681) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (5.63880983064e05 0 7.52344972099e06) force (2825915.41232 0 377041.133837) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (2.25552393226e05 0 3.0093798884e06) force (1130366.16493 0 150816.453535) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (9.02209572903e06 0 1.20375195678e06) force (452146.465971 0 60326.5814851) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (3.60883829161e06 0 4.81500782357e07) force (180858.586388 0 24130.6325762) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (1.44353531664e06 0 1.92600312587e07) force (72343.4345553 0 9652.25301269) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (5.77414126658e07 0 7.70401253902e08) force (28937.3738221 0 3860.90122288) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (2.30965650663e07 0 3.08160501561e08) force (11574.9495289 0 1544.36048915) moment (0 0 0) not converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged Constraint fixedLine1 error (9.23862602652e08 0 1.23264189966e08) force (4629.97981154 0 617.744142247) moment (0 0 0) converged Constraint fixedAxis1 angle 0 force (0 0 0) moment (0 0 0) converged sixDoFRigidBodyMotion constraints converged in 10 iterations Constraint force: (29433532.2251 0 3927099.98159) Constraint moment: (38738.9313677 630311.39261 0) Centre of mass: (3.63452602707e08 4.61882629007e07 11.3919999948) Linear velocity: (3.84183318467e07 1.3029827429e05 2.47155315804e07) Angular velocity: (9.92921402604e09 3.39167470597e07 1.35764137427e07) Orientation: (1 9.87594129196e10 1.23281745943e08 9.87594128113e10 1 8.77995122748e11 1.23281745943e08 8.77995000996e11 1) Euler angles: (8.77995122748e11 1.23281745943e08 9.87594129196e10) Execution time for mesh.update() = 214.59 s time step continuity errors : sum local = 1.14901701909e07, global = 2.51398523649e10, cumulative = 3.09829389781e05 GAMGPCG: Solving for pcorr, Initial residual = 1, Final residual = 0.00778738257516, No Iterations 1 GAMGPCG: Solving for pcorr, Initial residual = 0.00692915558121, Final residual = 0.00013905757398, No Iterations 2 time step continuity errors : sum local = 1.09535009227e09, global = 2.52645030719e11, cumulative = 3.09829137136e05 Courant Number mean: 0.00544396271522 max: 0.397594961785 velocity magnitude: 8.78304056663 MULES: Solving for alpha1 Liquid phase volume fraction = 0.565937062613 Min(alpha1) = 5.77953282643e25 Max(alpha1) = 1.00008410779 MULES: Solving for alpha1 Liquid phase volume fraction = 0.565920905776 Min(alpha1) = 1.50126700999e22 Max(alpha1) = 1.00007889189 MULES: Solving for alpha1 Liquid phase volume fraction = 0.565904748939 Min(alpha1) = 1.34168685408e22 Max(alpha1) = 1.00007387603 GAMG: Solving for pd, Initial residual = 0.000456710576829, Final residual = 1.52726175675e06, No Iterations 2 GAMG: Solving for pd, Initial residual = 1.64329099672e06, Final residual = 1.64329099672e06, No Iterations 0 time step continuity errors : sum local = 2.44125357119e06, global = 1.51108067885e08, cumulative = 3.09980245204e05 GAMG: Solving for pd, Initial residual = 1.70322629837e06, Final residual = 1.70322629837e06, No Iterations 0 GAMG: Solving for pd, Initial residual = 1.70322629837e06, Final residual = 1.70322629837e06, No Iterations 0 time step continuity errors : sum local = 2.53029278378e06, global = 3.100707169e08, cumulative = 3.10290315921e05 GAMG: Solving for pd, Initial residual = 1.70693050507e06, Final residual = 1.70693050507e06, No Iterations 0 GAMGPCG: Solving for pd, Initial residual = 1.70693050507e06, Final residual = 7.52057072396e08, No Iterations 2 time step continuity errors : sum local = 1.11724706312e07, global = 2.02980661923e10, cumulative = 3.10292345727e05 DILUPBiCG: Solving for epsilon, Initial residual = 0.860443214355, Final residual = 7.32310869309e08, No Iterations 3 DILUPBiCG: Solving for k, Initial residual = 0.754911596303, Final residual = 9.05265195856e08, No Iterations 3 ExecutionTime = 1074.3 s ClockTime = 1091 s laplaseSBStressDisplacement: Interface Courant Number mean: 2.08817723812e08 max: 0.0185460642856 Courant Number mean: 6.87160108967e05 max: 0.349542791659 deltaT = 0.000174156445196 Time = 0.0138694 sixDoFRigidBodyMotion constraints converged in 3 iterations Constraint force: (64487998.4846 0 5539864.3848) Constraint moment: (0 0 0) Centre of mass: (1.27390699003e09 8.06065326688e07 11.3919999999) Linear velocity: (5.83587203568e09 2.7689595378e05 1.64589850408e09) Angular velocity: (1.65356983894e05 1.00019045262e05 4.49188174605e07) Orientation: (1 6.85648172833e10 6.71106184628e08 6.85651256307e10 1 4.59461478782e08 6.71106184313e08 4.59461479242e08 1) Euler angles: (4.59461478782e08 6.71106184628e08 6.85648172833e10) GAMG: Solving for cellDisplacementx, Initial residual = 0.00720968353153, Final residual = 1.91358325779e08, No Iterations 3 GAMG: Solving for cellDisplacementy, Initial residual = 0.0285346506588, Final residual = 3.64593231008e08, No Iterations 3 GAMG: Solving for cellDisplacementz, Initial residual = 0.0173456488167, Final residual = 7.17668065688e08, No Iterations 2 Execution time for mesh.update() = 3.5 s time step continuity errors : sum local = 2.27644433709e10, global = 6.2416247832e15, cumulative = 7.46683578252e12 GAMGPCG: Solving for pcorr, Initial residual = 1, Final residual = 5.03987327995e08, No Iterations 14 GAMGPCG: Solving for pcorr, Initial residual = 0.0202290342766, Final residual = 5.69562969465e08, No Iterations 8 time step continuity errors : sum local = 2.49864598147e10, global = 2.19616475645e19, cumulative = 7.46683600213e12 MULES: Solving for alpha1 Phase1 volume fraction = 0.565726291588 Min(alpha1) = 1.20738285791e16 Max(alpha1) = 1.00000000004 MULES: Solving for alpha1 Phase1 volume fraction = 0.565726291587 Min(alpha1) = 1.20301624124e16 Max(alpha1) = 1.00000000004 MULES: Solving for alpha1 Phase1 volume fraction = 0.565726291587 Min(alpha1) = 1.19866388138e16 Max(alpha1) = 1.00000000004 GAMG: Solving for p_rgh, Initial residual = 3.43466231407e06, Final residual = 1.31067977375e08, No Iterations 3 GAMG: Solving for p_rgh, Initial residual = 1.02334734839e07, Final residual = 4.13903094003e09, No Iterations 2 time step continuity errors : sum local = 5.99217004483e13, global = 6.97485426546e15, cumulative = 7.4738108564e12 GAMG: Solving for p_rgh, Initial residual = 1.53852919472e08, Final residual = 2.13417817015e09, No Iterations 1 GAMGPCG: Solving for p_rgh, Initial residual = 3.47303297889e09, Final residual = 3.47303297889e09, No Iterations 0 time step continuity errors : sum local = 5.0280953268e13, global = 5.66583750124e15, cumulative = 7.4794766939e12 smoothSolver: Solving for epsilon, Initial residual = 0.0226768400394, Final residual = 3.62783564708e11, No Iterations 2 smoothSolver: Solving for k, Initial residual = 0.0477686014334, Final residual = 8.33785770767e11, No Iterations 2 ExecutionTime = 957.89 s ClockTime = 1090 s ================================================== =========================================== As can be seen from the above, the laplaseFaceDecomposition takes much much longer clock time than its counterpart. I am not sure whether this performance is typical, or maybe there's something tricky slipped by my eyes when playing with the parameters. Would be much appreciated if anyone could help to diagnose! Here's my tetFemSolution: motionU { solver ICCG; preconditioner { preconditioner FDIC; tolerance 1e6; relTol 0; smoother FDICGaussSeidel; nPreSweeps 1; nPostSweeps 3; nBottomSweeps 3; cacheAgglomeration false; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } tolerance 1e06; relTol 0; nSweeps 3; } dynamicMeshDict: dynamicFvMesh dynamicMotionSolverFvMesh; motionSolverLibs ("libfvMotionSolver.so"); solver laplaceFaceDecomposition; diffusivity quadratic inverseDistance 1; distancePatches ( HULL ); frozenDiffusion off; twoDMotion no; and lastly, the motionU: internalField uniform (0 0 0); boundaryField { Symmetry_Y_LOW { type fixedValue; value uniform (0 0 0); } Current_Z_LOW { type fixedValue; value uniform (0 0 0); } HULL { type sixDoFRigidBodyVelocity; centreOfMass (0 0 11.392); momentOfInertia (2.6422E10 2.5584E10 4.3501E10); mass 4.2791E7; rhoInf 1; // needed only for solvers solving for kinematic pressure report on; value uniform (0 0 0); constraints { maxIterations 500; fixedLine1 { sixDoFRigidBodyMotionConstraint fixedLine; tolerance 1e7; relaxationFactor 0.7; fixedLineCoeffs { refPoint (0 0 11.392); direction (0 1 0); } } fixedAxis1 { sixDoFRigidBodyMotionConstraint fixedAxis; tolerance 1e07; relaxationFactor 0.7; fixedAxisCoeffs { refPoint (0 0 11.392); axis ( 0 0 1 ); } } } } Inlet_X_LOW { type fixedValue; value uniform (0 0 0); } Outlet_Pressure_X_HIGH { type fixedValue; value uniform (0 0 0); } Outlet_Pressure_Z_HIGH { type fixedValue; value uniform (0 0 0); } Symmetry_Y_HIGH { type fixedValue; value uniform (0 0 0); } } 
Hi Xin,
A couple of observations in that context. 1. The size of the matrix for the laplaceFaceDecomposition is considerably larger than for the standard cellmotion utilities. The former solves a set of equations in a finite element context, i.e. for each point in the computational mesh. However, as the name suggests a decomposition of the computational mesh is taking place. For the facedecomposition (cannot recall for the celldecomposition), the total number of elements in your matrix are nPoints + nFaces + nCells and then substract a bit due to boundary conditions. So in terms of numbers, let us assume that you have a cube with N^3 cells, then you will have (N+1)^3 points and 3 * N^2(N+1) faces, meaning that the laplaceFaceDecomposition solves a system with order of magnitude 5N^3 elements plus a bit in the order N^2. Obviously, this will make iteration time slower. 2. Second point is your choice of linear solvers. It is a bit unfair towards the laplaceFaceDecomposition, that you are using GAMG on the cellmotion method, whereas you apply PCG on the other. For large meshes it is well known that GAMG by far outperforms PCG, e.g. you can see that the PCG do not even converge! Unfortunately, there are no GAMG solver available for the laplaceFaceDecomposition, but it would be a significant contribution. 3. Is it only due to differences in the outputting, that you do not seem to have constraints added for the mesh motion in the case of the cellmotion method? This might cause a difference in the systems you are comparing in terms of complexity, accelerations of the object, etc. Kind regards Niels 
Hi Niels,
Thanks. As for your point 3, i confirm that a line and axis constrains are applied to both. For some reasons inside the code, which I haven't dig out why, the details for each iteration of the constrains are reported in faceDecomposition of 1.6ext, but not in cellmotion of 2.1.1. And I set the "report" option on for both solvers. For point 2, I'd like to know how much work will need to be done to write a GAMG solver for faceDecomposition? I once found following segment from a tutorial provided by Dr. Jasak and now cannot recall the source: motionU { solver amgSolver; cycle Wcycle; policy AAMG; nPreSweeps 0; nPostSweeps 2; groupSize 4; minCoarseEqns 30; nMaxLevels 100; scale on; smoother GaussSeidel; minIter 0; maxIter 500; tolerance 1e7; relTol 0.0; }; I thus guess there's already been done by someone out there. Best, Xin 
And anyone please help to elaborate more on the general rule of the usage of "libs" section, as well as functionObjectLibs command in the controlDict file. When exactly do one need to include a lib explicitly and when do not? many thanks in advance!

Quote:
Do you have idea why the different types of DyMFoams and tools recognize my BC when working with faceDecomposition by setting DFACE_DECOMP, but throw out error message like "unkown boundary type ..." when I am trying to use cellDecomposition by setting DCELL_DECOMP? For example, decomposePar dosen't allow me to use my BC compiled with DCELL_DECOMP. I am pretty sure setFields will no doubt give same error message if I try to set my motionU field as well. Regards, Xin 
Hi XIn,
Thanks for pointing my attention on the working AMG motion solver. With respect to DCELL_DECOMP, then are you changing this only in your Make/options of your boundary conditions? In case of yes, then it is because OF1.6ext are compiler _either_ with face or cell decomposition. Kind regards, Niels 
Quote:
I got the best performance by using this configure: motionU { solver amgSolver; cycle Vcycle; policy AAMG; nPreSweeps 2; nPostSweeps 2; groupSize 4; minCoarseEqns 20; nMaxLevels 100; scale on; smoother DICGaussSeidel; minIter 0; maxIter 100; tolerance 1e7; relTol 0; }; And it's performance is as follows: ================================================== ====================== Courant Number mean: 0.00537268930754 max: 0.0920963615186 velocity magnitude: 1.74018304312 deltaT = 0.00141156462585 Time = 0.102602 Correct mesh motion diffusion field. amgSolver: Solving for motionUx, Initial residual = 0.0143791197762, Final residual = 6.33224460301e08, No Iterations 13 amgSolver: Solving for motionUy, Initial residual = 0.0123148790208, Final residual = 8.59546457804e08, No Iterations 12 amgSolver: Solving for motionUz, Initial residual = 0.0389445726386, Final residual = 7.11929446358e08, No Iterations 16 Execution time for mesh.update() = 9.66 s ================================================== ====================== To make a direct comparison to the ICCG solver: ================================================== ====================== Courant Number mean: 0.00537275160328 max: 0.0920963606907 velocity magnitude: 1.74018284511 deltaT = 0.00141156462585 Time = 0.102602 Correct mesh motion diffusion field. DICPCG: Solving for motionUx, Initial residual = 0.0143801344452, Final residual = 9.53697128507e08, No Iterations 57 DICPCG: Solving for motionUy, Initial residual = 0.0123160324063, Final residual = 9.36460962979e08, No Iterations 56 DICPCG: Solving for motionUz, Initial residual = 0.038947341729, Final residual = 9.36730620876e08, No Iterations 67 Execution time for mesh.update() = 5.49 s ================================================== ====================== To be fair, they are extracted from the same point of the simulation, and strangely enough the amgSolver managed convergence in much lesser iterations, but in almost double the CPU time. Does this imply at this mesh size (70K), ICCG solver is more cost efficient? Or is there anything I've missed? 
It does indeed look like a valid conclusion. My general experience is that the mesh motion solvers are just quite slow  at least the tetFem with face decomposition.
Did you have a change to test the cell decomposition, because I have remembered how this works, since this will only decompose to the cell centre, i.e. the number of points in the matrix will be (N+1)^3 + N^3 why it is only 40% of the face decomposed mesh motion. Kind regards Niels 
Hi
I have just tried utilising the constraints for the tetMesh motion, and I have found that it accelerated my mesh update by 25%, purely because I but a constraint on the velocities in the resolved boundary layer, so the points close to the boundary moves as a solid motion (currently, only a simple 2D problem). The reason for the acceleration is that I obtain a better initial residual, so the required number of iterations are smaller. Approximately 30% of the time is still used on creating the matrix itself. You would need to do a bit of programming, since my method is mostly a hardcoded fix to resolve one specific problem. Have a nice weekend, Niels 
I just gave up on this facedecomposition method because of performance issue. I did try the celldecomp method, but it's even worse than the face one, simply because the quality of the tet mesh due to cell splitting is much worse and this leads to a ill conditioned matrix.
PS. Ngj, Please have a look of my newly opened topic on your wave2Foam. I got problem regarding the relaxationZone. many thanks. 
All times are GMT 4. The time now is 11:02. 