[Ya_Squall2010]

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.6-extension 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.6-ext |
| \\ / A nd | Web: www.extend-project.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.6-ext |
| \\ / A nd | Web: www.extend-project.de |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 1.6-ext
Exec : setFields -time 0
Date : Mar 27 2013
Time : 15:14:25
Host : ubuntu
PID : 43873
Case : /home/ya-squall/OpenFOAM/OpenFOAM-1.6-ext/tutorials/multiphase/interDyMFoam/ras/floatingObject
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

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

--> FOAM Warning :
From function dlLibraryTable:pen(const fileName& functionLibName)
in file db/dlLibraryTable/dlLibraryTable.C at line 86
could not load /home/ya-squall/OpenFOAM/OpenFOAM-1.6-ext/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.

[Ya_Squall2010]

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.

[Ya_Squall2010]

The previous one is wrong... sorry for inconvenience...

Attachment 20180

[ngj]

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 OF16-ext.

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)

Good luck,

Niels

[Ya_Squall2010]

Originally Posted by ngj

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 OF16-ext.

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)

Good luck,

Niels

Hi Niels,

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.22629688705e-06, No Iterations 1000
DICPCG: Solving for motionUy, Initial residual = 0.172646722126, Final residual = 1.21688108815e-05, No Iterations 1000
DICPCG: Solving for motionUz, Initial residual = 0.0634080697595, Final residual = 5.68417442456e-06, No Iterations 1000
Constraint fixedLine1 error (0.000352425614415 0 -4.70215607589e-05) force (-17661971.327 -0 2356507.08661) moment (0 0 0) not converged
Constraint fixedAxis1 angle 2.10734242554e-08 force (0 0 0) moment (-38738.9313677 630311.39261 0) converged
Constraint fixedLine1 error (0.000140970245766 0 -1.88086243043e-05) 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.63880983064e-05 0 -7.52344972099e-06) 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.25552393226e-05 0 -3.0093798884e-06) 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.02209572903e-06 0 -1.20375195678e-06) 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.60883829161e-06 0 -4.81500782357e-07) 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.44353531664e-06 0 -1.92600312587e-07) 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.77414126658e-07 0 -7.70401253902e-08) 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.30965650663e-07 0 -3.08160501561e-08) 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.23862602652e-08 0 -1.23264189966e-08) 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.63452602707e-08 4.61882629007e-07 11.3919999948)
Linear velocity: (3.84183318467e-07 1.3029827429e-05 -2.47155315804e-07)
Angular velocity: (9.92921402604e-09 -3.39167470597e-07 -1.35764137427e-07)
Orientation: (1 -9.87594129196e-10 1.23281745943e-08 9.87594128113e-10 1 8.77995122748e-11 -1.23281745943e-08 -8.77995000996e-11 1)
Euler angles: (-8.77995122748e-11 1.23281745943e-08 9.87594129196e-10)

Execution time for mesh.update() = 214.59 s
time step continuity errors : sum local = 1.14901701909e-07, global = -2.51398523649e-10, cumulative = 3.09829389781e-05
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.09535009227e-09, global = -2.52645030719e-11, cumulative = 3.09829137136e-05
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.77953282643e-25 Max(alpha1) = 1.00008410779
MULES: Solving for alpha1
Liquid phase volume fraction = 0.565920905776 Min(alpha1) = -1.50126700999e-22 Max(alpha1) = 1.00007889189
MULES: Solving for alpha1
Liquid phase volume fraction = 0.565904748939 Min(alpha1) = -1.34168685408e-22 Max(alpha1) = 1.00007387603
GAMG: Solving for pd, Initial residual = 0.000456710576829, Final residual = 1.52726175675e-06, No Iterations 2
GAMG: Solving for pd, Initial residual = 1.64329099672e-06, Final residual = 1.64329099672e-06, No Iterations 0
time step continuity errors : sum local = 2.44125357119e-06, global = 1.51108067885e-08, cumulative = 3.09980245204e-05
GAMG: Solving for pd, Initial residual = 1.70322629837e-06, Final residual = 1.70322629837e-06, No Iterations 0
GAMG: Solving for pd, Initial residual = 1.70322629837e-06, Final residual = 1.70322629837e-06, No Iterations 0
time step continuity errors : sum local = 2.53029278378e-06, global = 3.100707169e-08, cumulative = 3.10290315921e-05
GAMG: Solving for pd, Initial residual = 1.70693050507e-06, Final residual = 1.70693050507e-06, No Iterations 0
GAMGPCG: Solving for pd, Initial residual = 1.70693050507e-06, Final residual = 7.52057072396e-08, No Iterations 2
time step continuity errors : sum local = 1.11724706312e-07, global = 2.02980661923e-10, cumulative = 3.10292345727e-05
DILUPBiCG: Solving for epsilon, Initial residual = 0.860443214355, Final residual = 7.32310869309e-08, No Iterations 3
DILUPBiCG: Solving for k, Initial residual = 0.754911596303, Final residual = 9.05265195856e-08, No Iterations 3
ExecutionTime = 1074.3 s ClockTime = 1091 s


laplaseSBStressDisplacement:

Interface Courant Number mean: 2.08817723812e-08 max: 0.0185460642856
Courant Number mean: 6.87160108967e-05 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.27390699003e-09 8.06065326688e-07 11.3919999999)
Linear velocity: (-5.83587203568e-09 2.7689595378e-05 1.64589850408e-09)
Angular velocity: (1.65356983894e-05 1.00019045262e-05 -4.49188174605e-07)
Orientation: (1 -6.85648172833e-10 6.71106184628e-08 6.85651256307e-10 1 -4.59461478782e-08 -6.71106184313e-08 4.59461479242e-08 1)
Euler angles: (4.59461478782e-08 6.71106184628e-08 6.85648172833e-10)

GAMG: Solving for cellDisplacementx, Initial residual = 0.00720968353153, Final residual = 1.91358325779e-08, No Iterations 3
GAMG: Solving for cellDisplacementy, Initial residual = 0.0285346506588, Final residual = 3.64593231008e-08, No Iterations 3
GAMG: Solving for cellDisplacementz, Initial residual = 0.0173456488167, Final residual = 7.17668065688e-08, No Iterations 2
Execution time for mesh.update() = 3.5 s
time step continuity errors : sum local = 2.27644433709e-10, global = 6.2416247832e-15, cumulative = -7.46683578252e-12
GAMGPCG: Solving for pcorr, Initial residual = 1, Final residual = 5.03987327995e-08, No Iterations 14
GAMGPCG: Solving for pcorr, Initial residual = 0.0202290342766, Final residual = 5.69562969465e-08, No Iterations 8
time step continuity errors : sum local = 2.49864598147e-10, global = -2.19616475645e-19, cumulative = -7.46683600213e-12
MULES: Solving for alpha1
Phase-1 volume fraction = 0.565726291588 Min(alpha1) = -1.20738285791e-16 Max(alpha1) = 1.00000000004
MULES: Solving for alpha1
Phase-1 volume fraction = 0.565726291587 Min(alpha1) = -1.20301624124e-16 Max(alpha1) = 1.00000000004
MULES: Solving for alpha1
Phase-1 volume fraction = 0.565726291587 Min(alpha1) = -1.19866388138e-16 Max(alpha1) = 1.00000000004
GAMG: Solving for p_rgh, Initial residual = 3.43466231407e-06, Final residual = 1.31067977375e-08, No Iterations 3
GAMG: Solving for p_rgh, Initial residual = 1.02334734839e-07, Final residual = 4.13903094003e-09, No Iterations 2
time step continuity errors : sum local = 5.99217004483e-13, global = -6.97485426546e-15, cumulative = -7.4738108564e-12
GAMG: Solving for p_rgh, Initial residual = 1.53852919472e-08, Final residual = 2.13417817015e-09, No Iterations 1
GAMGPCG: Solving for p_rgh, Initial residual = 3.47303297889e-09, Final residual = 3.47303297889e-09, No Iterations 0
time step continuity errors : sum local = 5.0280953268e-13, global = -5.66583750124e-15, cumulative = -7.4794766939e-12
smoothSolver: Solving for epsilon, Initial residual = 0.0226768400394, Final residual = 3.62783564708e-11, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.0477686014334, Final residual = 8.33785770767e-11, 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 1e-6;
relTol 0;
smoother FDICGaussSeidel;
nPreSweeps 1;
nPostSweeps 3;
nBottomSweeps 3;
cacheAgglomeration false;
nCellsInCoarsestLevel 10;
agglomerator faceAreaPair;
mergeLevels 1;
}
tolerance 1e-06;
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 1e-7;
relaxationFactor 0.7;
fixedLineCoeffs
{
refPoint (0 0 11.392);
direction (0 1 0);
}
}

fixedAxis1
{
sixDoFRigidBodyMotionConstraint fixedAxis;
tolerance 1e-07;
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);
}
}

[ngj]

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 cell-motion 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 face-decomposition (cannot recall for the cell-decomposition), 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 cell-motion 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 cell-motion method? This might cause a difference in the systems you are comparing in terms of complexity, accelerations of the object, etc.

Kind regards

Niels

[Ya_Squall2010]

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.6-ext, but not in cell-motion 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 W-cycle;
policy AAMG;
nPreSweeps 0;
nPostSweeps 2;
groupSize 4;
minCoarseEqns 30;
nMaxLevels 100;
scale on;
smoother GaussSeidel;

minIter 0;
maxIter 500;
tolerance 1e-7;
relTol 0.0;
};

I thus guess there's already been done by someone out there.

Best,
Xin

[Ya_Squall2010]

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!

[Ya_Squall2010]

Originally Posted by ngj

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 OF16-ext.

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)

Good luck,

Niels

Hi Niels,

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

[ngj]

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.6-ext are compiler _either_ with face or cell decomposition.

Kind regards,

Niels

[Ya_Squall2010]

Originally Posted by ngj

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.6-ext are compiler _either_ with face or cell decomposition.

Kind regards,

Niels

You are most welcome. As pointed out by Dr. Jasak however, the AMG solver is only available for series mode. And it's performance is questionable:

I got the best performance by using this configure:

motionU
{
solver amgSolver;
cycle V-cycle;
policy AAMG;
nPreSweeps 2;
nPostSweeps 2;
groupSize 4;
minCoarseEqns 20;
nMaxLevels 100;
scale on;
smoother DICGaussSeidel;

minIter 0;
maxIter 100;
tolerance 1e-7;
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.33224460301e-08, No Iterations 13
amgSolver: Solving for motionUy, Initial residual = 0.0123148790208, Final residual = 8.59546457804e-08, No Iterations 12
amgSolver: Solving for motionUz, Initial residual = 0.0389445726386, Final residual = 7.11929446358e-08, 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.53697128507e-08, No Iterations 57
DICPCG: Solving for motionUy, Initial residual = 0.0123160324063, Final residual = 9.36460962979e-08, No Iterations 56
DICPCG: Solving for motionUz, Initial residual = 0.038947341729, Final residual = 9.36730620876e-08, 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?

[ngj]

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

[ngj]

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 2-D 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 hard-coded fix to resolve one specific problem.

Have a nice weekend,

Niels

[Ya_Squall2010]

I just gave up on this face-decomposition method because of performance issue. I did try the cell-decomp 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.