Octree grid structure and grid density function in Palabos v2

Hello community,

I am trying to understand the octree grid structure from the new release of Palabos, found in examples/showCases/gridRefinement3d. This is the call to the constructor and to the function generateGridStructure() in the example:

[code=“cpp”]
void createOctreeGridStructure(SimulationParameters &param)
{
bool useSamples = false;
plint numSamples = -1;
plint maxIter = 100;
bool removeBlocks = true;
bool fineToCoarse = true;
int numLevelsToGroupBlocks = 0;
int numLevelsToGroupOverlaps = -1;
bool strongGrouping = false;
bool verbose = true;
bool stlOutput = true;
std::string stlBaseName = “octree”;
T gridDensityScaleFactor = (T)1;

OctreeGridGenerator<T> octreeGridGenerator(
        param.fullDomain,
        param.gridDensityFunctionFile,
        param.minLeafLevel,
        param.maxLeafLevel,
        param.nBlock,
        global::mpi().getSize(),
        gridDensityScaleFactor,
        useSamples,
        numSamples,
        maxIter,
        removeBlocks,
        fineToCoarse,
        numLevelsToGroupBlocks,
        numLevelsToGroupOverlaps,
        strongGrouping,
        param.outDir,
        verbose,
    stlOutput,
        stlBaseName);

param.ogs = octreeGridGenerator.generateOctreeGridStructure();
param.fullDomain = octreeGridGenerator.getFullDomain();
param.dxFinest = octreeGridGenerator.getDxFinestLevel();

}



In the grid density function file, gridDensityFunction.dat, the first line is the domain in physical units(?), the second line is dx in physical units(?), third line is the Nx, Ny and Nz of the outer cuboid. The fourth line contains Nx*Ny*Nz points between 0 and 1 (or negative values in regions where the user does not want blocks to be allocated) that gives the octree grid generator the criteria to refine the domain or not. 

I still do not understand how to generate the values for the file gridDensityFunction.dat, and how to manipulate it to generate my own octree grid structure. Can I just give values of 1 where I want the refinement and 0 where I don't want it?

For the output, I also do not understand some other things. What does maxOutputLevel mean? How can we output 100 levels while when we generated the octree we specified 3 levels with the parameter numLevels?

Can octrees be implemented for moving geometries with immersed boundary method? For the moment I would like to implement the case where the finest grid envelops the whole region of rotation. 

Thank you very much for your help!

Patricia

Hello Patricia,

Here is the basic code I use to create the grid density function file. This version allows to recreate the mesh needed for the given example.
In can be adapted to suit other simulations.

Hope it can save some people time.

(Probably some useless includes^^)

[code=“cpp”]
#include
#include
#include
#include
#include
#include
#include
#include <dirent.h>
#include <errno.h>
#include

using namespace std;
typedef double T;

int main()
{

/////////////////////////////////////////////////////////////////////////////////////////////
// Initialization
/////////////////////////////////////////////////////////////////////////////////////////////
// Physical coordinates of the zone we want to impose refinement
const T x0 = -4.0;
const T x1 = 11.0;
const T y0 = -4.0;
const T y1 = 4.0;
const T z0 = -4.0;
const T z1 = 4.0;

// Physical length spacing.
const T dx = 0.1;

// Number of nodes of the grid density function, in each direction.
const int nx = (x1-x0)/dx + 2;
const int ny = (y1-y0)/dx + 2;
const int nz = (z1-z0)/dx + 2;

// Temporary variable definition.
T x01; int nx01;
T x11; int nx11;
T y01; int ny01;
T y11; int ny11;
T z01; int nz01;
T z11; int nz11;

// Definition of the grid density array.
// It is set to zero (no refinement) by default.
T ***myGridFunction = new T**[nx];
for (int i=0; i<nx; i++){
    myGridFunction[i] = new T*[ny];
    for (int j=0; j<ny; j++){
        myGridFunction[i][j] = new T[nz];
        for (int k=0; k<nz; k++){
            myGridFunction[i][j][k] = 0.0;
        }
    }
}

/////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////////////////////
// Zone of single refinement
/////////////////////////////////////////////////////////////////////////////////////////////
// Physical coordinates of the refined domain.
x01 = -4.0; y01 = -4.0; z01 = -4.0;
x11 = 11.0; y11 = 4.0; z11 = 4.0;

// Definition of first and last node of refinement, in each direction.
nx01 = (x01-x0)*nx / (x1-x0); ny01 = (y01-y0)*ny / (y1-y0); nz01 = (z01-z0)*nz / (z1-z0); 
nx11 = (x11-x0)*nx / (x1-x0); ny11 = (y11-y0)*ny / (y1-y0); nz11 = (z11-z0)*nz / (z1-z0); 
// Loop on this refinement domain.
for (int iX=nx01; iX<nx11; iX++){
    for (int iY=ny01; iY<ny11; iY++){
        for (int iZ=nz01; iZ<nz11; iZ++){
            // Assign grid density function to 0.5, ensuring medium refinement.
            myGridFunction[iX][iY][iZ] = 0.5;
        }
    }
}

/////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////////////////////
// Zone of double refinement
/////////////////////////////////////////////////////////////////////////////////////////////
// Physical coordinates of the refined domain.
x01 = -4.0; y01 = -4.0; z01 = -4.0;
x11 = 4.0; y11 = 4.0; z11 = 4.0;

// Definition of first and last node of refinement, in each direction.
nx01 = (x01-x0)*nx / (x1-x0); ny01 = (y01-y0)*ny / (y1-y0); nz01 = (z01-z0)*nz / (z1-z0); 
nx11 = (x11-x0)*nx / (x1-x0); ny11 = (y11-y0)*ny / (y1-y0); nz11 = (z11-z0)*nz / (z1-z0); 
// Loop on this refinement domain.
for (int iX=nx01; iX<nx11; iX++){
    for (int iY=ny01; iY<ny11; iY++){
        for (int iZ=nz01; iZ<nz11; iZ++){
            // Assign grid density function to 1.0, ensuring max refinement.
            myGridFunction[iX][iY][iZ] = 1.0;
        }
    }
}

/////////////////////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////////////////////
// Output to file
/////////////////////////////////////////////////////////////////////////////////////////////
string nameOutFile = “gridDensityFunction.dat”;
ofstream outFile;
outFile.open(nameOutFile.c_str());
outFile << x0 << “\t” << x1 << “\t”<< y0 << “\t” << y1 << “\t” << z0 << “\t” << z1 << “\n”;
outFile << dx << “\n”;
outFile << nx << “\t”<< ny << “\t”<< nz << “\n”;
for (int iX=0; iX<nx; iX++){
for (int iY=0; iY<ny; iY++){
for (int iZ=0; iZ<nz; iZ++){
outFile << myGridFunction[iX][iY][iZ] << “\t”;
}
}
outFile << “\n”;
}
outFile.close();
/////////////////////////////////////////////////////////////////////////////////////////////

return 0;

}

1 Like

Hi Philippe,

I ended up writting a Python code that does pretty much the same, writing “-1” inside the sphere for voxelizing the lattice in that region, and writing “1” around, where we want the refinement. I have added “0.5” for the wake were we want medium refinement:

[code=“python”]
import numpy as np
import matplotlib.pyplot as plt

Input parameters

l = np.array([-4, 11, -4, 4, -4, 4])

Nx = 152 ; Ny = 82 ; Nz = 82

dx = 0.1

fname = “gridDensityFunction.dat”

N = np.array([Nx,Ny,Nz])
rho = np.zeros((Nx,Ny,Nz))
density = []
dx_ = (l[1]-l[0])/Nx
dy_ = (l[3]-l[2])/Ny
dz_ = (l[5]-l[4])/Nz

Initialize

x0 = l[0]
y0 = l[2]
z0 = l[4]

x = np.zeros(Nx)
y = np.zeros(Ny)
z = np.zeros(Nz)

for i in range(0,Nx):
x[i] = x0 + i * dx_
for j in range(0,Ny):
y[j] = y0 + j * dy_
for k in range(0,Nz):
z[k] = z0 + k * dz_
if((-0.5 <= y[j] <= 0.5 and -0.5 <= z[k] <= 0.5) and -0.5 <= x[i] <= 0.5 ):
rho[i,j,k] = -1
elif((-0.65 <= y[j] <= 0.65 and -0.65 <= z[k] <= 0.65) and -0.65 <= x[i] <= 0.65 ):
rho[i,j,k] = 1
elif((-0.8 <= y[j] <= 0.8 and -0.8 <= z[k] <= 0.8) and 0.65 <= x[i] <= 0.8 ):
rho[i,j,k] = 0.5
else:
rho[i,j,k] = 0
density.append(rho[i,j,k])

Write output for palabos in fname

with open(fname,‘w’) as f:
f.write(str(l)[1:-1]) # 1:-1 to get rid of []
f.write(’\n’)
f.write(str(dx))
f.write(’\n’)
f.write(str(N)[1:-1])
f.write(’\n’)
s = str(density)[1:-1]
f.write(s.replace(’,’,’’))




Thanks for sharing!

Patricia
1 Like

Hello to both of you,

thank you for sharing your codes for the generation of the grid density function. You are right that 1 stands for maximal refinement and 0 for places we want the lowest resolution (and all the values in between are used for intermediate refinement depending on the number of levels one wants). Finally -1 is for places where one would like blocks not to be allocated.

Note that the resolution of your grid density function is not necessarily related with the resolution of your simulation. Furthermore if the external cuboid of your density function is smaller than the domain of your simulation a 0 value is assumed there.

Finally you may notice some differences between you grid density function and the actual octree that is generated. This is usually due to the fact that for consistency of the refinement algorithm, one must ensure that each refined grid is only refined by a factor two and that each level must be composed of at least one block. Therefore the interfaces between levels (especially when one has many) may be larger that what you would expect.

For your other two questions:

The maxOutputLevel flag is for writing the vtk files output up to a certain level (if it is larger than the maximum refinement level it is automatically taken as the maximum level).

Moving boundaries can be used with the grid refinement. You must ensure that all your boundary surface is on the same refinement level (boundaries crossing the refinement interface are forbidden). “Simply” put your boundaries on the finest level and it should work without too much problems.

Do not hesitate to file any bug reports.

Best regards
Orestis

Hello Orestis,

Thank you for these precisions.

I am working on a domain much longer than wide and high, but I would like to keep cubic cuboids.
For that I am deallocating external cuboids. It works well but implies to define my grid density function on the whole domain.
To properly feet my geometries, the grid density function resolution needs to be of the order of the smaller cuboid length.
These 2 conditions result in very heavy grid density function, too heavy when reaching 8th refinement level.

Is there any way of deallocating external cuboids without using the grid density function?

Best regards,

Philippe

Hello Philippe,

For the moment there is not “automatic” way to deallocate any part of the domain, but it could be an interesting feature to add at some point.

Nevertheless, the octree generation is performed on some domain which does not have to be the same as the domain on which a grid density function is defined. For example the octree grid generation can be performed on a [-200,200]^3 domain while the grid density function is only defined on a [-10,10]x[-5,5]x[-2,2] domain. Would it do the trick?

Cheers,
Orestis

Hallo everybody,
thanks a lot for sharing and the useful informations!
Does anyone know if it is possible to implement a free surface field in the case of Octree multigrid? Have you got experience with the topic?
Best regards and thanks for the help!
Raffaele

Hello Patricia,

Can you use the octrees (3D refinement )for moving geometries with immersed boundary method in Palabos?
best wishes
Lin