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