CFD Online Discussion Forums
Page 1 of 2 1 2
Show 40 post(s) from this thread on one page

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Improve simpleFoam convergence (https://www.cfd-online.com/Forums/openfoam-solving/78354-improve-simplefoam-convergence.html)

Daniele111 July 19, 2010 14:08
Improve simpleFoam convergence
 
Hi
How can I improve simpleFoam convergenge? Now time simulation is 30min, can i improve it?

Thanks

matejfor July 20, 2010 02:40
Sure you can. Have a faster computer.
No, well, yes, you can make it run faster but you need to give us more information. Your question could be paraphrased as: "I have a car and it takes 30 minutes to get to work, could it be faster?"

Post here fvSchemes, fvSolution, and screenshot of the residual plot, tell us what is the case (turbulent? what is Re?,backward facing step or swirler or flow around car, plane, building, channel flow), what are the BCs and if the variables (like velocity) are oscillating (e.g. post a plot of a velocity probe).

matej

Daniele111 July 20, 2010 05:25
2 Attachment(s)
:) you're right!! And I like your metaphor! It is a turbulent case (k-epilon model). It is similar to a backward facing step, i have a atmospheric bl with a obstacle on bottom.

top bottom outlet inlet
k symmetry wallFunction zerogradient fixed value
U symmetry 0 zerogradient log profile
p symmetry zerogradient 0 zerogradient
eps symmetry wallfunction zerogradient imposed profile

Thanks

MartinB July 20, 2010 05:59
Hi Daniele,

I would reduce the relaxation factors to:
Code:
relaxationFactors
{
    p              0.3;
    U              0.7;
    k              0.5;
    epsilon        0.5;
    R              0.5;
    nuTilda        0.5;
}
And I would change the solvers to:
Code:
    p
    {
        solver          GAMG;
        tolerance        1e-7;
        relTol          0.001;
        minIter          5;
        maxIter          100;
        smoother        GaussSeidel; // DIC; //DICGaussSeidel; //FDIC;
        nPreSweeps      1;
        nPostSweeps      3;
        nFinestSweeps    3;
        scaleCorrection true;
        directSolveCoarsest false;
        cacheAgglomeration on;
        nCellsInCoarsestLevel 50;    // 500
        agglomerator    faceAreaPair;
        mergeLevels      1;    // 3
    }
   
    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance        1e-6;
        relTol          0.01;
        nSweeps          1;
        maxIter        100;
    };
The solvers for k, epsilon, R and nuTilda the same as for U.

Martin

Daniele111 July 20, 2010 06:44
1 Attachment(s)
Thanks
The execution time is reduced from 2239 sec to 662, but after the obstacle the solution it isn't accurate

Daniele

MartinB July 20, 2010 06:53
Hi Daniele,
your convergence criteria in fvSolution with 1e-3 is rather low...
Restart your simulation from the latest timestep with convergence criteria of 1e-5.
Martin

Daniele111 July 20, 2010 07:03
Thanks
I try it. Uhm In my project the obstacle on the bottom change its position, and then i calculate new solution at new position, this is do in a evolutive algorithm, so the iteration after the first must be very fast less than a minute. can i initialize the new mesh with the last time of previous simulation?

Daniele

MartinB July 20, 2010 07:25
I think, you can use the "mapFields" utility.
With "mapFields -help" you can see the parameters necessary.
In your "system" directory you need a file "mapFieldsDict":

Code:
/*--------------------------------*- C++ -*----------------------------------*\\
| =========                |                                                |
| \\      /  F ield        | OpenFOAM: The Open Source CFD Toolbox          |
|  \\    /  O peration    | Version:  1.6                                  |
|  \\  /    A nd          | Web:      www.OpenFOAM.org                      |
|    \\/    M anipulation  |                                                |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version    2.0;
    format      ascii;
    class      dictionary;
    location    "system";
    object      mapFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// This dictionary is necessary to map fields from another simulation to this
// one. The results from the source will be interpolated to the mesh in this
// case directory.
// Usually patchMap and cuttingPatches will be empty.
// In the case that you want to use the output of one simulation as the
// input for the next one, define patchMap as (inlet outlet).

patchMap        (  );  // usually empty
cuttingPatches  (  );

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Navigate to your new (!) case directory.

Use the mapFields command like this:
mapFields /home/daniele/workinghere/firstRun -sourceTime 1234

Check the mapped results with paraFoam.

Start the simulation in your new case directory.

Martin

Daniele111 July 20, 2010 09:54
Hi
If I start new simulation with last time of the previous simulation time execution is 93s... I would about 30s 40s; if I use mapfield time execution is 400s. There is a thing that I dont understand: I start a simulation, then I copy last time directory in 0 directory, and restart the same case, and time execution is 30s, but if I start from solution because time is so high?

MartinB July 20, 2010 10:52
It's not necessary to copy last time directory to "0" directory, you can edit the "system/controlDict" to start from your latest time directory:
Code:
startFrom      startTime;              //  firstTime, startTime, latestTime
startTime      0;                      //  set > 0 to continue
Replace "startFrom startTime;" with "startFrom latestTime;" and the simulation continues.
Or replace "startTime 0;" with "startTime 1000;" to start from iteration 1000.

You can use the mapField command with the parameter "-consistent" if your mesh is exactly the same. Otherwise if the meshes are different, the values are interpolated and the solver needs some time steps to adjust the solution again.

May be I didn't understand your questions right... you don't speak german, do you?

Martin

Daniele111 July 20, 2010 11:28
I understand. I try it. Two mesh are very similar, top inlet and outlet are identical, the obstacle position on bottom si traslate of some cm.
I did a simulation, and then I start a new simulation with the same mesh but startTime is latestTime of previous sim, and time execution is 40s, is normal? I don't change nothing, and i start from convergence solution, time execution should not be faster?

p.s i don't speak german, and my english is very bad...you are very kind, and your advice are very useful.

MartinB July 20, 2010 12:23
It's difficult to say if 40s is normal: usually one can observe that the time necessary for a single iteration step (i.e. solving for Ux, Uy, Uz, p ...) is higher in the beginning of the simulation (at least for a steady state simpleFoam simulation). With reaching convergence the time for each single iteration step drops significantly.

Next to the time used for a single iteration step the number of iterations used to reach convergence is important. For my simulations here (simpleFoam without turbulence) I get convergence in a "fresh" simulation after 800 iterations. When using mapFields I can reduce the number of iterations to 400 iterations, although the first few single iteration steps do need the same execution time.

Another important issue in your case is the question: when is your simulation converged good enough? I think a convergence tolerance of 1e-3 is not enough, but may be, 1e-4 is fine and 1e-5 is an overkill? You can save quite a lot of iterations as soon as you know which accuracy you really need...

Some more ideas to speed up your simulations:
Are you running your simulation in parallel mode? If not, and if you have a multicore or multiprocessor machine, use the parallelism!
An important parameter which reduce the number of iterations is the relaxation factor: you can try to increase these factors!
If you are familiar with OpenFOAM's c++ programming language you might try to use an dynamically increasing relaxation factor during runtime. There is at least one thread in this forum where this idea is discussed.

If you can't speed up the simulations time directly you can try to use your hardware more efficiently: create a shell script or a python script that runs dozens of simulations over night. I use this technique quite often... automatic mesh creation with blockMesh is not hard, and with the "sample" utility you can do automatic postprocessing, too.

Hope these ideas are useful

Martin

MartinB July 20, 2010 14:06
Here are some more ideas:
I suppose you are using hexahedra elements. If not, try switching to hexahedra.

Reduce the mesh size, i.e. the number of elements. Run your simulation until convergence. Then use a finer mesh and map the previous results...

Optimize your mesh so it really fits the needs of your case: a fine mesh in the important areas (behind the obstacle for example), a coarser mesh elsewhere.

With the refineMesh utility you can refine the whole mesh immediately. With refineWallLayer you can split the cells next to a patch. If you have hexahedra and you want to refine locally only, I can upload an utility...

Although this is usually done to avoid divergence in complex cases: run simpleFoam without turbulence first, at least for let's say 200 iterations. Then switch on turbulence and continue the simulation.

For my previous suggestion of using mapFields you can try this: just use U and p for the mapping and leave the turbulence fields.

Of course you can combine some of these strategies, too...

Martin

matejfor July 21, 2010 02:50
Daniele,
martin gave you extensive advice to which there's not much to add.
Just a few bits:

a) the precision of the result - using upwind scheme for convection terms will gain you speed, but not precision. The good news is, you may change the scheme on the fly using some utility or script.
Other issue with precision might be the use of k-eps model and consequently wall functions in OpenFOAM.

b) I understand that you have solved solution, then you move the obstacle, map thesolution and run. well your solution is converged, but the flowfield is different as the obstacle is moved so the solution will take a little longer at the beginning. Alternativelly to Martin's suggestion with switching off the turbulence, I would make it bit easier with changing the viscosity which means changing the turbulence level.

You are after two rabbits here. The precision one is going one direction while the speedy one is running slightly other direction. It's hard to get both at the same time.


good luck
matej

Daniele111 July 21, 2010 03:00
Dear Martin

Thanks for everything. You were very patient and kind to me. I try different and valuable advice you gave me.

Daniele

Daniele111 July 22, 2010 07:08
Hi
I refine my mesh only where is necessary and i reduce my domain size. I calculated a first solution with a convergence tolerance of 1e-4, and solution is very good. I try to increase tolerance (1e-3,1e-2) for the following simulation that I inizialized with mapFields, and the time simulation is very good 20s, and 10s, and the solution quality is good. Is it wrong calculate the new solution from a solution more precise,but raising tolerance?

Daniele

MartinB July 22, 2010 07:24
Hi Daniele,
it is not wrong, but you should be aware of the results quality!
Try this: write your results every 50th iteration. You can do this with the following entry in "system/controlDict":
writeInterval 50;

Then check the results in paraview, i.e. have a careful look at how much the results change from let's say timestep 200 to timestep 250, or timestep 200 to 800... it is useful to reduce the number of colors to 10 (in Color Scale Editor, enable "Discrete Colors" with Solution of 10). You can see the amount of change in the fields... then decide for your particular problem, if it is ok or not.

I am doing flow channel shape optimization here, and I can say quite early if the new geometry is better or worse than the latter one. But when reaching the optimum, I must decrease the convergence tolerances more and more.

Martin

matejfor July 22, 2010 07:28
It is completely wrong. You've chosen the wrong rabbit to chase. The tolerance we are talking about is the tolerance with which the solver solves the linearized set of N-S equations. When you release the tolerance, the solver will be faster, but the price....
You can play with a tolerance but when you want reliable solution you should lower it back towards the end of the solution of each of your case.

good luck
matej

Daniele111 July 22, 2010 07:38
Dear Matey
But if i inizialized my simulation with a field very precise, I must not use necessarily a very low tolerance or not?

Daniele111 July 22, 2010 08:02
1 Attachment(s)
Post my solution. The source for mapfields utility is the same for all. I change tolerance convergence.

MartinB July 22, 2010 09:52
Quote:
Dear Matey
But if i inizialized my simulation with a field very precise, I must not use necessarily a very low tolerance or not?
Yes, you should do! The precise field for initialization gets "unprecise" immediately on the new mesh.

Quote:
Post my solution. The source for mapfields utility is the same for all. I change tolerance convergence.
Although your results for tau/tau_0 look stable after reaching 1e-2, it really is a very high tolerance criteria! I use at least 1e-4 for simpleFoam. The other results (U for example) will differ much more, I suppose. Well, I hope you are not running nuclear power plants or passenger air planes with your simulations ;)

Daniele111 July 22, 2010 10:58
No I'm not simulating a nuclear power plant! :D But also u profile is good, i think that I'll use 1e-3 tolerance convergence. I know the importance of tolerance i made structural simulation for a wing with nastran and I know the its importance

Daniele111 July 22, 2010 12:15
Hi
Do you know a way to interface my openfoam simulation with a matlab script that use OpenFoam results?

MartinB July 22, 2010 14:46
The best interface is interaction via files, I think...
I found this thread where you had problems starting OpenFOAM from Matlab:
http://www.cfd-online.com/Forums/ope...de-matlab.html

Do you still have problems in doing so? If so, try
Code:
    unix('foamExec simpleFoam')
since this should set your environment variables properly.

Otherwise
Code:
    unix('source path/to/OpenFOAMs/bashrc; simpleFoam')
could work.

For the data you need in Matlab:
Do you only need values for shear stress on the wall?
Is it enough to have coordinates for each shear stress value, or do you need the full mesh structure?
Do you use a more sophisticated transport model than newtonian flow?
Are you using OpenFOAM in parallel?

I have some code snippets here that would write a simple ";"-separated .csv file which you can parse in Matlab.

Daniele111 July 22, 2010 15:08
If I use your command I have this:

unix('foamExec blockMesh')
/bin/bash: foamExec: command not found


I only need values for shear stress on the wall, and it is enough to have coordinates for each shear stress value;

"Do you use a more sophisticated transport model than newtonian flow?" ->> No
Are you using OpenFOAM in parallel? ->> No (maybe later)

Thanks
Daniele

MartinB July 22, 2010 15:24
What about the other command: unix('source path/to/OpenFOAMs/bashrc; simpleFoam') ?

"path/to/OpenFOAMs" must be replaced with the full path to your OpenFOAM's installation directory, than the "bashrc" file follows, which includes the environment for OpenFOAM. Than there is a ";", followed by another command (simpleFoam) to be executed in the just set environment...

What's the name of the patch, you need the shear stress values for? Is it "bottom"?

Daniele111 July 22, 2010 15:33
I do:
unix('/home/acconcia/OpenFOAM/OpenFOAM-1.6.x/etc/bashrc; blockMesh')
and

/bin/bash: /home/acconcia/OpenFOAM/OpenFOAM-1.6.x/etc/bashrc: permission denied
/bin/bash: blockMesh: command not found

Patch name is "bottom"

How long do you work with OpenFoam? Do you use it for university or work? I'm star to use is one month ago, and I like it, it's very flexible.

Daniele111 July 22, 2010 16:37
unix('/home/acconcia/OpenFOAM/OpenFOAM-1.6.x/etc/bashrc; blockMesh')

Now It works, I changed matlab permission

Thanks

MartinB July 22, 2010 16:50
Hi Daniele,

here we go:

Go to the simpleFoam source directory:
"~/OpenFOAM/OpenFOAM-1.6.x/applications/solvers/incompressible/simpleFoam/"
The version number might be different for your system.

Create a new file there with name "write.H". Paste this code in there:
Code:
// write matlab file
{
    Info << "Writing wall shear stress for matlab!" << endl;

    word wallPatchName = "bottom";
    Info << "Searching for patch " << wallPatchName << endl;

    // Check, if the patch name exists:
    label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
    if (patchID < 0)
    {
        Info << "No patch found with name \"" << wallPatchName << "\"!"
            << endl;
    }
    else
    {
        /* for laminar case:
        Info<< "Reading viscosity nu\n" << endl;
        dimensionedScalar nu_
        (
            laminarTransport.lookup("nu")
        );
        */
        // calculating wall shear rate
        volVectorField wallGradU
        (
                IOobject
                (
                        "wallGradU",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                mesh,
                dimensionedVector
                (
                        "wallGradU",
                        U.dimensions()/dimLength,
                        vector::zero
                )
        );

        forAll(wallGradU.boundaryField(), patchi)
        {
            wallGradU.boundaryField()[patchi] =
                -U.boundaryField()[patchi].snGrad();
        }

        // calculating wall shear stress
        volVectorField tau
        (
                IOobject
                (
                        "tau",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                //laminar case: nu_ * wallGradU
                turbulence->nuEff() * wallGradU
                // turbulence->nu() * wallGradU
        );

        // Pointer to the patch we want to have shear stress for:
        const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
        label nFaces = cPatch.size();

        pointField faceCenterCoords(nFaces, vector(0,0,0));
        pointField faceAreas(nFaces, vector(0,0,0));
        pointField wallShearStress(nFaces, vector(0,0,0));

        label runFaces = 0;

        // collecting the interesting data
        for (int i=0; i<nFaces; i++)
        {
            faceCenterCoords[runFaces] = mesh.Cf().boundaryField()[patchID][i];
            faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
            wallShearStress[runFaces] = wallGradU.boundaryField()[patchID][i];
            runFaces++;
        }

        //  make a cute filename
        fileName casePath = cwd();//casePath.name();
        fileName caseName = casePath.name();
        fileName myFile = caseName+"_wallShearStress.csv";
        Info << "Data is written to file " << myFile << endl;

        // create a new file
        OFstream *myStream = NULL;
        myStream = new OFstream(myFile);

        for (int i=0; i<nFaces; i++)
        {
            *myStream << i << "; "
                    // coordinates on face centers:
                    << mesh.Cf().boundaryField()[patchID][i].component(0) << "; "
                    << mesh.Cf().boundaryField()[patchID][i].component(1) << "; "
                    << mesh.Cf().boundaryField()[patchID][i].component(2) << "; "

                    // wall shear stress as a vector:
                    << tau.boundaryField()[patchID][i].component(0) << "; "
                    << tau.boundaryField()[patchID][i].component(1) << "; "
                    << tau.boundaryField()[patchID][i].component(2) << "; "

                    // magnitude of wall shear stress:
                    << mag(tau.boundaryField()[patchID][i]) << endl;
        }

    }
}
Then open the simpleFoam.C file. Add this line to the includes at the top:
Code:
#include "OFstream.H"
Behind the runTime loop of simpleFoam paste the following two lines, so that it looks like:
Code:
        ...
        #include "convergenceCheck.H"
    }
    // writing additional data for matlab
#  include "write.H" // <--- this is the important one
    ...
In the simpleFoam directory call
Code:
wclean
wmake
Now if you are running simpleFoam an extra file with name "your_case_name_wallShearStress.csv" is written. It's format is:
number; coordX; coordY; coordZ; wallShearStressX; wallShearStressY; wallShearStressZ; magnitude(wallShearStress)

I haven't tested this code, so use it with care! Check and double check, if the wall shear stress is calculated really correctly!!!
Have a look a line 62, where the wall shear stress is calculated:
Code:
turbulence->nuEff() * wallGradU
is just a guess, since I don't have experience with turbulence. It might be nu, or nut or something totally different there...

Have fun :-)

Martin

p.s.: I'm using OpenFOAM for seven month now... but I'm still learning, of course ;-) And I'm using it for work.

Daniele111 July 23, 2010 04:04
Dear Martin
Thanks for your code, is a gold mining! There are many things I want to learn to. I'll change the computation of wallShearStress, I think that i'll use the openfoam code for this and I add it to write.h.

Daniele

aldo.iannetti July 26, 2010 11:05
MRF simple foam linear upwind: convergence troubles
 
1 Attachment(s)
Hi,
I'm trying to set a linearUpwind MRF simpleFoam simulation of a cooling fan. The mesh has been imported from GAMBIT (8.5 e+6 cells, not bad quality) and the simulation diverges.Here attached: log, fvSchemes,fvSolution. Could you please have a look in order to give me advices about settings?
thanks
Aldo

Daniele111 August 21, 2010 06:24
But if I want extract also k with code? How can i modify the code post previously?

MartinB August 21, 2010 08:37
Hi Daniele,

what about adding
Code:
<< k.boundaryField()[patchID][i] << "; "
to the "myStream" loop?
Martin

Daniele111 August 21, 2010 08:47
Hi Martin
It is that I do...I did a syntax error probably.
Thanks!

Daniele

Daniele111 November 4, 2010 19:12
Hi
Now I try to use a cluster to increase my speed, and so i'm traying tu use OpenFoam parallel utility. But when i run decomposePar i have this error:

FOAM FATAL IO ERROR:

Cannot find 'value' entry on patch bottom of field epsilon in file "/u/acconcia/OpenFOAM/acconcia-1.6/run/blayer4/0/epsilon"
which is required to set the values of the generic patch field.
(Actual type epsilonWallFunction)

Please add the 'value' entry to the write function of the user-defined boundary-condition
or link the boundary-condition into libfoamUtil.so

file: /u/acconcia/OpenFOAM/acconcia-1.6/run/blayer4/0/epsilon::boundaryField::bottom from line 49 to line 52.

From function genericFvPatchField<Type>::genericFvPatchField(con st fvPatch&, const Field<Type>&, const dictionary&)
in file genericFvPatchField/genericFvPatchField.C at line 71.


My boundary condition are the same of the case that I use not in parallel and they haven't give me problem.

matejfor November 4, 2010 23:50
The best think you can do is to follow what the error message clearly says so (especially when it does say it so clearly). On patch "bottom" you have not prescribed value of epsilon. What is the BC for epsilon on bottom? It is probably some which needs a value.

Now why you do not need it in single but need it in parallel?
My guess is that you do not use k-eps but k-omega where BC for epsilon is not needed.
The decompose utility does not know what you will use for the run so it does decompose everything and checks for you BCs which is very kind of it.
But it's only my guess.


good luck

matej

Daniele111 November 5, 2010 03:45
1 Attachment(s)
These are my BC. And I'm sure that I'm using k-eps model, this BC have never been a problem.

Thanks

matejfor November 5, 2010 04:04
OK. then you have 2 possibilities:

(1) add the value uniform number; line to your bottom BC for the decomposer to be happy
(end remove it afterwards as the value is not needed by the wall functions (you can check it in doxygen)

(2) remove the need of it from the decomposer, recompile it and if it works report it


matej

Daniele111 November 5, 2010 10:40
Hi
I solved my problem, paralle is happy with a value!
To calculate wall shear stress with single processor I used Martin's code of post 29, but in parallel I can't use it: is possible a similar code to have a output at simulation stop? A tau.dat like output then i'll manipulate with matlab.

MartinB November 5, 2010 14:16
Hi Daniele,

not fully tested, but it should be a point to start from:

Code:
// write matlab file
{
    Info << "Writing wall shear stress for matlab!" << endl;

    word wallPatchName = "bottom";
    Info << "Searching for patch " << wallPatchName << endl;

    // Check, if the patch name exists:
    label patchID = mesh.boundaryMesh().findPatchID(wallPatchName);
    if (patchID < 0)
    {
        Info << "No patch found with name \"" << wallPatchName << "\"!"
            << endl;
    }
    else
    {
        /* for laminar case:
        Info<< "Reading viscosity nu\n" << endl;
        dimensionedScalar nu_
        (
            laminarTransport.lookup("nu")
        );
        */
        // calculating wall shear rate
        volVectorField wallGradU
        (
                IOobject
                (
                        "wallGradU",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                mesh,
                dimensionedVector
                (
                        "wallGradU",
                        U.dimensions()/dimLength,
                        vector::zero
                )
        );

        forAll(wallGradU.boundaryField(), patchi)
        {
            wallGradU.boundaryField()[patchi] =
                -U.boundaryField()[patchi].snGrad();
        }

        // calculating wall shear stress
        volVectorField tau
        (
                IOobject
                (
                        "tau",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                ),
                //laminar case: nu_ * wallGradU
                turbulence->nuEff() * wallGradU
                // turbulence->nu() * wallGradU
        );
       
        // Pointer to the patch we want to have shear stress for:
        const polyPatch& cPatch = mesh.boundaryMesh()[patchID];
        label nFaces = cPatch.size();

        pointField faceCenterCoords(nFaces, vector(0,0,0));
        pointField faceAreas(nFaces, vector(0,0,0));
        pointField wallShearStress(nFaces, vector(0,0,0));

        label runFaces = 0;

        // collecting the interesting data
        for (int i=0; i<nFaces; i++)
        {
            faceCenterCoords[runFaces] = mesh.Cf().boundaryField()[patchID][i];
            faceAreas[runFaces] = mesh.Sf().boundaryField()[patchID][i];
            wallShearStress[runFaces] = tau.boundaryField()[patchID][i];
            runFaces++;
        }

        // for parallel
        reduce(nFaces, sumOp<label>());
        Info << "Number of global faces: " << nFaces << endl;
       

        //  make a cute filename
        fileName casePath = cwd();//casePath.name();
        fileName caseName = casePath.name();
        fileName myFile = caseName+"_wallShearStress.csv";
        Info << "Data is written to file " << myFile << endl;

        // create a new file
        OFstream *myStream = NULL;
        if (Pstream::master())
        {
            myStream = new OFstream(myFile);
        }

        // slaves send data to master, master writes data to file
        if (Pstream::parRun())
        {
            if (Pstream::myProcNo() != Pstream::masterNo())
            {
                OPstream toMaster(Pstream::blocking, Pstream::masterNo());
                toMaster << faceCenterCoords;
                toMaster << faceAreas;
                toMaster << wallShearStress;
            }
            else // master part
            {
                // first write own data
                for (int i=0; i<faceCenterCoords.size(); i++)
                {
                    *myStream << i << ";"
                        << faceCenterCoords[i].component(0) << "; "
                        << faceCenterCoords[i].component(1) << "; "
                        << faceCenterCoords[i].component(2) << "; "

                        << wallShearStress[i].component(0) << "; "
                        << wallShearStress[i].component(1) << "; "
                        << wallShearStress[i].component(2) << "; "

                        << mag(wallShearStress[i]) << endl;
                }

                label run = faceCenterCoords.size();
               
                // then receive slave data and write
                pointField bufferFaceCenterCoords;
                pointField bufferFaceAreas;
                pointField bufferWallShearStress;

               
                for (int slave=Pstream::firstSlave(); slave <= Pstream::lastSlave(); slave++)
                {
                    IPstream fromSlave(Pstream::blocking, slave);
                    fromSlave >> bufferFaceCenterCoords;
                    fromSlave >> bufferFaceAreas;
                    fromSlave >> bufferWallShearStress;
                   
                    for (int i=0; i<bufferFaceCenterCoords.size(); i++)
                    {
                        *myStream << run << ";"
                            << bufferFaceCenterCoords[i].component(0) << "; "
                            << bufferFaceCenterCoords[i].component(1) << "; "
                            << bufferFaceCenterCoords[i].component(2) << "; "

                            << bufferWallShearStress[i].component(0) << "; "
                            << bufferWallShearStress[i].component(1) << "; "
                            << bufferWallShearStress[i].component(2) << "; "

                            << mag(bufferWallShearStress[i]) << endl;
                        run++;
                    }
                }
            }
        }
        else // single processor only
        {
            for (int i=0; i<nFaces; i++)
            {
                *myStream << i << "; "
                        // coordinates on face centers:
                        << mesh.Cf().boundaryField()[patchID][i].component(0) << "; "
                        << mesh.Cf().boundaryField()[patchID][i].component(1) << "; "
                        << mesh.Cf().boundaryField()[patchID][i].component(2) << "; "

                        // wall shear stress as a vector:
                        << tau.boundaryField()[patchID][i].component(0) << "; "
                        << tau.boundaryField()[patchID][i].component(1) << "; "
                        << tau.boundaryField()[patchID][i].component(2) << "; "

                        // magnitude of wall shear stress:
                        << mag(tau.boundaryField()[patchID][i]) << endl;
            }
        }

    }
}
Martin


All times are GMT -4. The time now is 08:07.
Page 1 of 2 1 2
Show 40 post(s) from this thread on one page