Tutorial On Richardson Zaki Cylinder Batch Sedimentation Experiment

1. Generate the Mesh

Since the grid for a cylinder is intrinsically unstructured, it requires a few steps to generate the kind of easy-to-partition mesh necessary for FCPU to make use of the loadNP2Sieve feature. If one wants to simply generate the whole big grid with GMSH, many of the following steps are not necessary. But this usually means very large memory request and very slow partitioning procedure. What is descried below is very efficient and thus highly recommended.

1.1 GMSH

Use the following script to generate the original mesh

[sjin1@moondance H8NP2]$ cat R10dx0.4Layer1.geo
//Only Change the following two lines
radius = 10;
element_dim = 0.4; // control the elements' size

//Do NOT Change below
length = element_dim;
extrusion_number = 1; // number of elements in the extrusion direction, = (number of planes -1),=length/element_dim
//Do NOT Change below
Mesh.MshBinary = 0;
Mesh.Format = 1;
Mesh.SaveAll = 0; //Ignore Physical definitions and save all elements
Mesh.ElementOrder = 1;
Mesh.SecondOrderLinear = 1;
Mesh.Optimize = 1; //Optimize the mesh using Netgen to improve the quality of tetrahedral elements
Mesh.Algorithm = 2; //1 by Antoine (1=meshadapt, 2=delaunay)
Mesh.Algorithm3D = 1;//4 by Antoine //(1=delaunay, 4=netgen, 5=tetgen)

////////////   No Need to Change Beyond ////////////////////////////
// points used to create quarter of circle
Point(1) = {radius, 0, 0, element_dim} ;
Point(2) = {0, radius, 0, element_dim} ;
Point(3) = {0, 0, 0, element_dim};

Circle(1) = {1,3,2};
//Delete { Point {3}; }

// copy and rotate the quarter circle to create a full one
Rotate {{0,0,1}, {0,0,0}, Pi/2} {
  Duplicata { Line{1}; }
}

Rotate {{0,0,1}, {0,0,0}, Pi} {
  Duplicata { Line{1}; }
}

Rotate {{0,0,1}, {0,0,0}, 3*Pi/2} {
  Duplicata { Line{1}; }
}

Line Loop(1) = {1,2,3,4};
Plane Surface(1) = {1};

// extrude in third direction
e1[]=Extrude { 0 , 0 , length }
{
        Surface{ 1 };
        Layers
        {
                { extrusion_number },
                { 1. }
        };
};

Physical Volume("domain") = {e1[1]};

Note that only the first a few lines need to be changed for a specific cylinder. To explain what we are doing exactly, we need to read this code. Most importantly,

length = element_dim;
extrusion_number = 1;

these two lines says that we are only generating mesh for a thin cylinder with only one element size high.
The reason we are doing this is that GMSH numbers the elements and vertexes quite randomly (it seems to me random but maybe there is a specific algorithm to compute them. At least they are not in a simple continuous manner), which makes it impossible to take advantage of the repetitive nature of the cylinder mesh and use the loadNP2Sieve options.
Here we will first generate the mesh for a single layer and use our own tools to generate the repetitive grid.

Once the .geo file is created, use gmsh to generate the .msh file. For example,

gmsh -3 R10dx0.4Layer1.geo

THe -3 option (means 3D) will not invoke the GUI of Gmsh and will simply generate the cylinderR10H8Dx0.4.msh file on the command line.
The generated mesh can be view in GUI by

gmsh cylinderR10H8Dx0.4.msh

1.2 Convert GMSH to FCPU

Then the gmsh .msh file should be converted using our conversion tool gmsh2fcpu to FCPU recognizable input mesh files

./gmsh2fcpu cylinderR10H8Dx0.4.msh

where

[sjin1@moondance H8NP2]$ ls -l gmsh2fcpu
lrwxrwxrwx 1 sjin1 users 55 Jan 7 00:56 gmsh2fcpu -> /home/sjin1/FCPU-XP-DEV/trunk/tools/Gmsh2fcpu/gmsh2fcpu

The gmsh2fcpu code can be found in the FCPU SVN code under trunk/tools/Gmsh2fcpu/.
You may see some warnings like

[sjin1@moondance dx0.4]$ ./gmsh2fcpu R10dx0.4Layer1.msh
R10dx0.4Layer1.msh is a gmsh file, version 2, ASCII format.
WARNING ! MALLOC : passed size is 0 !4876 nodes
14154 elements
0 boundary conditions sets
[sjin1@moondance dx0.4]$

This is because at end of the .geo file, we have only

Physical Volume("domain") = {e1[1]};

and don't have any boundary domains specified. This way gmsh will only output 3D elements (tetrahedra), not the 1D or 2D elements.

This way, the 3d.lcon 3d.mybc 3d.nodes files are generated, where the 3d.mybc file is not useful since we didn't specify boundary conditions. This is OK since the FCPU code can automatically generate BC for cylinders on the fly using the autoBC option.

1.3 Get Ready for Parallel Partition

Here we need to use our own tools in FCPU-XP-DEV/trunk/tools/sortMesh to prepare for parallel partition. There are two codes here (once "make" is done, they should appear): sortMesh and repeatLayer.
First, we use sortMesh to renumber of gmsh-generated mesh to have the elements and vertexes show up in order.
Simply run

/home/sjin1/FCPU-XP-DEV/trunk/tools/sortMesh/sortMesh

and you will find two new files generated:

sorted.lcon sorted.nodes

These two are nothing but a re-order of 3d.lcon and 3d.nodes.

Secondly, we call the repeatLayer program to generate larger domain meshes for NP=2 parallel partition.

[sjin1@moondance dx0.4]$ /home/sjin1/FCPU-XP-DEV/trunk/tools/sortMesh/repeatLayer  20
repeat 20 times
51198 new nodes
283080 new elements

Note that here you need to supply an integer as the argument, which dictates how many times the single layer generated above is to be repeated along the axis.
You need to choose this number based on the height for NP=2 parition, not the final production domain size. For example, we want to have H=4 on each process in the final parallel job, so H=8 for the NP=2 partitioning run (see next). Therefore, H/dx=8/0.4=20.

Upon completion, a new directory is created and three files can be found inside

[sjin1@moondance dx0.4]$ ls 20layersH8NP2/
3d.aux 3d.lcon 3d.nodes

1.4 Run NP2 outputSieveOnly

Now get into the 20layersH8NP2 directory generated above, and first check the contents of 3d.aux to make sure the radius is correct since my code does not check that (R=1 always).

[sjin1@moondance dx0.4]$ cd 20layersH8NP2/
[sjin1@moondance 20layersH8NP2]$ cat 3d.aux
dx=0.4
$GEO-TYPE cylinder
xc=0,yc=0,zc=0
R=10,H=8

Now run the following script

[sjin1@moondance H8NP2]$ cat run2
#!/bin/bash
CODE=/home/sjin1/FCPU-XP-DEV/trunk/fcpu
nohup mpiexec  -n 2 -l  $CODE -Re 1.0 -Fr 0.0098 -dt 0.05 -NSteps 0 \
  -outputSieveOnly 1 -subSteps 1  -curvedBC 1\
        -identicalVelStructure 1 -dumpStep 10   &> out2 &
[sjin1@moondance H8NP2]$

and the 3d_* files are created.

2. Production Run

Suppose we decide to run a cylinder sedimentation of H=32 with NP=8. Usually this is done in a different directory.

2.1 Particle Input File

First make sure a new 3d.aux is represent with

dx=0.4
$GEO-TYPE cylinder
xc=0,yc=0,zc=0
R=10,H=32

and then run

/home/sjin1/FCPU-XP-DEV/trunk/tools/randPar3D 1.1 10000

You may need to go to /home/sjin1/FCPU-XP-DEV/trunk/tools/ and run make to build the randPar3D code first.
The first option 1.1 specifies the minimal center-to-center distance between spheres and it should be larger than 1.
The second option specifies the number of spheres to generate.
For the first time, we don't have a good idea of how many can be randomly placed inside the given cylinder, so a large number such as 10000 is used.
From its out, we can obtain an idea on the maximum amount (you may need to kill the run with Ctrl+C when it is giving more output) and then redo it with proper amount such as

/home/sjin1/FCPU-XP-DEV/trunk/tools/randPar3D 1.1 4800

Once finished, a file named D1.1NP4800.xyz is generated with x, y and z coordinates for each particle listed per line. 4800 particle with diameter 1 in this cylinder means a solid volume fraction of 25%.
Then invoke the following parInput3D.m Matlab script

function parInput3D(filename)
R=0.5;rhoa=1.2;rhob=1.2;
p3d=load(filename);
nSp=length(p3d(:,1));
fp=fopen('3d.particles','w');
fprintf(fp,'#nSp=%d\n',nSp);
fprintf(fp,'#rho_f=1, t=0\n');
fprintf(fp,'#1-id 2-radius 3-rho_p 4-fixed 5-x 6-y 7-z 8-u 9-v 10-z 11-w1 12-w2 13-w3 14-theta 15-phi\n');
rhoRand=round(rand(nSp,1));
for i=1:nSp
        if(rhoRand(i)<0.5)
                rhop=rhoa;
        else
                rhop=rhob;
        end
 
        fprintf(fp, '%d %f %f %d %f %f %f %f %f %f %f %f %f %f %f\n' ...
                ,i-1,R,rhop,0,p3d(i,1),p3d(i,2),p3d(i,3),0,0,0,0,0,0,0,0);
end
fclose(fp);

In Matlab, do

»parInput3D('D1.1NP4800.xyz')

and 3d.particles is created. Note that the values of rhoa and rhob can be changed to have two solid species with different solid densities while in this test, we only deal with a single specie.
Now putting all required 3d* files together, you should have

[sjin1@moondance H32NP8]$ ls 3d*
3d.aux                      3d_IsieveNP2Rank0.txt  3d_renumberingNP2Rank0.txt
3d_coordinatesNP2Rank0.txt  3d_IsieveNP2Rank1.txt  3d_renumberingNP2Rank1.txt
3d_coordinatesNP2Rank1.txt  3d.particles
[sjin1@moondance H32NP8]$

Usually the production runs are done on a different system so the above files need to be copied over.

2.2 Run Production Jobs

For example, the script used on the terminus.ucalgary.ca cluster looks like the following

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License