User Tools

Parametrized postprocessing with GMSH

Extracting ion flux onto rotatable targets

This section corresponds to 3D PIC-MC simulations of cylindrical sputter targets. In order to predict the sputter erosion profile, the normal component of the ion flux with respect to the target surface is required. The ion flux is defined as the product of ion density and velocity.

For a planar target, where the sputter direction is e.g. parallel to the Z axis, it is very trivial to obtain the averaged flux by following steps:

  • Merge the density_Arplus_xxxx.pos file into GMSH (–> View 0)
  • Merge the velocity_Arplus_xxxx.pos fils into GMSH (–> View 1)
  • Use the MathEval plugin of GMSH in order to create the product density and vz
  • Use the CutPlane plugin of GMSH in order to extract the plane just above the target surface

Alternatively, the absorption plot absorption_Arplus_xxxx.pos can be used in this case.

For cylindrical rotating targets, the absorption plot might be not sufficient, since - due to the use of a cartesian grid in PIC-MC - the cells adjacent to the target surface form a step-like structure rather than a smooth cylindrical surface.

Thus, the scripts introduced here have the following purposes:

  • Make a smooth projection of the ion flux onto the cylindrical surface
  • Perform a time average over the last N snapshots
  • Create POS files of the projected absorption on the curved target surface, which can be triangulated
  • Finally performing circular integration of the ion flux in order to predict the effective erosion profile.

First step: Create ASCII files of the ion flux in the relevant area

In order to work properly, the extraction script needs the ion fluxes as ASCII text file with following restrictions:

  • The data should be restricted to the region adjacent to the target surface of interest
  • The data points should have equal spacing in X, Y, and Z direction

In order to prepare the ASCII files, go into the folder of your PAR files and create the following directories:

  • ./extract
  • ./extract/flux

Within the ./extract directory, create a file make-flux.geo with a content like follows:

time = 0;
For time In {50 : 60 : 1.0} 
    Merge Sprintf("../frame-dual-3D/velocity_Arplus_%.2fus.pos", time);
    Merge Sprintf("../frame-dual-3D/density_Arplus_%.2fus.pos",  time);
    Plugin(MathEval).Expression0= "v0*w0";
    Plugin(MathEval).Expression1= "v0*w1";
    Plugin(MathEval).Expression2= "v0*w2";
    Plugin(MathEval).Expression3= "";
    Plugin(MathEval).Expression4= "";
    Plugin(MathEval).Expression5= "";
    Plugin(MathEval).Expression6= "";
    Plugin(MathEval).Expression7= "";
    Plugin(MathEval).Expression8= "";
    // Creating POS files is optional and can be omitted
    Save View[6] Sprintf("flux/flux_Arplus_%.2fus.pos", time);
    Save View[6] Sprintf("flux/flux_Arplus_%.2fus.txt", time);
    Delete View[6];
    Delete View[5];
    Delete View[4];
    Delete View[3];
    Delete View[2];
    Delete View[1];
    Delete View[0];

This script loads the density and the velocity plot of a given time step and performs the CutPlane plugin four times in order to restrict the resulting vector plot to the region around the target surface. It is recommended to perform these steps once in GMSH manually and to transfer the CutPlane parameters to the script above. Afterwards the script can be invoked with the Linux binary of GMSH in batch mode:

gmsh make-flux.geo -

After running the GMSH script, the flux files created under ./extract/flux should look like follows1):

Fig. 1: Cylindrical target geometry with extracted flux in region of interest

Second step: Make a projection of the flux onto the curved target surface

The following RIG-VM script assumes, that the target axis is parallel to the Z axis. If that is not the case in the actual geometry, there is the possibility to perform an internal coordinate rotation within the script.

Before working with the RIG-VM script, you have to find out the number of cells as well as the lowest and highest coordinate in each direction (X, Y, Z), which can be done by viewing one of the ASCII flux files within a text editor.

All of the required adjustments are to be made in the upper part of the script, as summarized in the following:

  • BASENAME: Template for input filenames. Change this e. g. if your ion species is different from Arplus
  • ColX, ColY, ColZ: Column offsets for the X, Y, and Z coordinate, should be either (0, 1, 2), (1, 2, 0), or (2, 0, 1). If the cylinder is already oriented along z-axis chose the first one (i. e. ColZ=2). If the cylinder is oriented along X axis in the GEO file, chose the combination, where ColZ=0. If it is oriented along Y axis chose ColZ=1.
  • IZ0, IZ1: The lowest and highest cell index along target axis, where the flux evaluation should be performed. Make sure, that the full racetrack is comprised within that range.
  • NX, NY, NZ: Number of cells in each coordinate direction. Important: If the coordinate system is rotated, that also applies for these values. Remember that the Z-axis is always parallel to the target axis…
  • X0, X1, Y0, Y1, Z0, Z1: The lowest and highest coordinate in each direction.
  • Alpha0, Alpha1, dAlpha: The lowest and highest angular value as well as the angular step in order to define the relevant circle arc on the target surface. Make sure that this region is fully comprised within the flux field.
  • XM, YM: The center of the cylinder in the XY plane.
  • R0: The radius of the cylindrical target
  • DR: A safety margin in the order of one half cell spacing in order to prevent, that empty cells are considered within the evaluation.

If all of these values are properly adjusted, the script, which is shown in the following, can be run with

rvm extract_flux.r
# extract_flux.r
# 3D flux files should be present under ./flux/flux_Arplus_*.txt
integer IZ0 = 13,   IZ1 = 86;   # Mininum and maximum Z coordinate index
float   XM  = 167.5, YM = 247;  # Center coordinate of cylinder in XY plane
string BASENAME : "flux_Arplus_$(time%,2f)us";  # Scheme of input filenames  
integer ColX=1, ColY=2, ColZ=0; # Column offsets for rotating the coordinate system
                                # If the cylindricat target in the geo file is already
                                # oriented in Z direction, use ColX=0, ColY=1, ColZ=2
# Number of entries in X, Y, Z
# Lowest and highest coordinate in X, Y, Z, respectively
integer NX = 70;  float X0 = 47.03214285714284, X1 = 258.4678571428572, DX = (X1-X0)/(NX-1);
integer NY = 37;  float Y0 = 267,               Y1 = 339,               DY = (Y1-Y0)/(NY-1);
integer NZ = 100; float Z0 = 3.75,              Z1 = 746.25,            DZ = (Z1-Z0)/(NZ-1);
float Alpha0 = -30, Alpha1=50, dAlpha=1;  # Extracted angular range
float R0     = 69;                        # Target diameter plus safety margin
float DR     = 2;                         # Safety margin in the order of at least one half cell spacing
# =======================================================================================================
# No more adjustments needed from this point on
# =======================================================================================================
string str12(x=float) = return "$("$(x)"%12s)";
integer dim_fn = (IZ1-IZ0+1)*((Alpha1-Alpha0)/dAlpha+1);
float FN_Sum[dim_fn];
float realX(x=float, y=float, z=float) = return ColX==0 ? x ! ColY==0 ? y ! z;
float realY(x=float, y=float, z=float) = return ColX==1 ? x ! ColY==1 ? y ! z;
float realZ(x=float, y=float, z=float) = return ColX==2 ? x ! ColY==2 ? y ! z;
block extract_flux(BASENAME=string, sum->var) : {
  float bilinear(f00=float, f10=float, f01=float, f11=float, x=float, y=float, x0=float, y0=float, x1=float, y1=float) = {
    float F = 1/((x1-x0)*(y1-y0));
    return F*(f00*(x1-x)*(y1-y) + f10*(x-x0)*(y1-y) + 
	      f01*(x1-x)*(y-y0) + f11*(x-x0)*(y-y0));
  float index(x=float, y=float, z=float) = 
    return floor((x-X0)/DX) + NX*floor((y-Y0)/DY) + NX*NY*floor((z-Z0)/DZ);
  float simeq(a=float, b=float) = return abs(a-b)<1e-6 ? 1 ! 0;
  console "** Reading $(BASENAME).txt **\n";
    startline=1; separator=" ";
  integer iz;
  float   zplane -> Z0 + DZ*iz;
  console "** Generating flux fields **\n";
  float FX[NX*NY*NZ], FY[NX*NY*NZ], FZ[NX*NY*NZ];
  integer i, ind;
  for(i=0; i<M.n_rows; i=i+1) {
    float x = floor(([4+ColX][i]-X0)/DX+0.5)*DX+X0+DX/10,
          y = floor(([4+ColY][i]-Y0)/DY+0.5)*DY+Y0+DY/10,
          z = floor(([4+ColZ][i]-Z0)/DZ+0.5)*DZ+Z0+DZ/10;
    float fx=[7+ColX][i],
    ind = index(x, y, z);  
    FX[ind] = fx; FY[ind] = fy; FZ[ind] = fz;
  for(iz=IZ0; iz<=IZ1; iz=iz+1) {
    console "** Processing plane ix=$(iz) **\n";
    if(iz==IZ0) then file "$(BASENAME)_EXT.pos" << out "View\"$(BASENAME)\" {\n";
    append "$(BASENAME)_EXT.pos" << {
      float alpha, pi = 4*atan(1), x, y, x0, y0, z0, nx, ny, tx, ty;
      integer ind;
      for(alpha=Alpha0; alpha<=Alpha1+1e-6; alpha=alpha+dAlpha) {
	x =  XM+(R0+DR)*sin(alpha*pi/180);
	y =  YM+(R0+DR)*cos(alpha*pi/180);
	x0 = floor((x-X0)/DX)*DX + X0;
	y0 = floor((y-Y0)/DY)*DY + Y0;
	ind = index(x, y, zplane+1e-3);
        fx = bilinear(FX[ind], FX[ind+1], FX[ind+NX], FX[ind+NX+1], x, y, x0, y0, x0+DX, y0+DY);
        fy = bilinear(FY[ind], FY[ind+1], FY[ind+NX], FY[ind+NX+1], x, y, x0, y0, x0+DX, y0+DY);
        fz = bilinear(FZ[ind], FZ[ind+1], FZ[ind+NX], FZ[ind+NX+1], x, y, x0, y0, x0+DX, y0+DY);
	# Normal unit vector
        nx = (x-XM)/sqrt((x-XM)^2+(y-YM)^2);
        ny = (y-YM)/sqrt((x-XM)^2+(y-YM)^2);
	# Tangential unit vector
	tx = -ny;
	ty = nx;
	fn = fx*nx + fy*ny;    # Normal flux component
	ft = fx*tx + fy*ty;    # Tangential flux component
	out "SP($(realX(x, y, zplane)), $(realY(x, y, zplane)), $(realZ(x, y, zplane))){$(-fn)};\n";
        # fn auf externes Array addieren
        integer isum = (alpha-Alpha0)/dAlpha + ((Alpha1-Alpha0)/dAlpha+1)*(iz-IZ0);
        sum[isum] = sum[isum] + fn;
    if(iz==IZ1) then append "$(BASENAME)_EXT.pos" << out "};\n";
  return 0;
integer count;
for(time=50.0; time<=60.0+1e-6; time=time+1) {
  extract_flux(BASENAME, FN_Sum);
for(i=0; i<FN_Sum.len; i=i+1) FN_Sum[i] = -(FN_Sum[i])/count;  # Flux unit is 1/m^2s
# Output of averaged flux as POS and TXT
file "AVG_flux_Arplus.pos" << {
  file "AVG_flux_Arplus.txt" << out "$("alpha"%12s)  $("z"%12s)  $("fn"%12s)\n";
  out "View\"AVG_Flux_Arplus\" {\n";
  integer iz;
  float   zplane -> Z0 + DZ*iz;
  for(iz=IZ0; iz<=IZ1; iz=iz+1) {
    float alpha, pi = 4*atan(1), x, y, x0, y0, z0, nx, ny, tx, ty;
    for(alpha=Alpha0; alpha<=Alpha1+1e-6; alpha=alpha+dAlpha) {
      x =  XM+(R0+DR)*sin(alpha*pi/180);
      y =  YM+(R0+DR)*cos(alpha*pi/180);
      integer isum = (alpha-Alpha0)/dAlpha + ((Alpha1-Alpha0)/dAlpha+1)*(iz-IZ0);
      out "SP($(realX(x, y, zplane)), $(realY(x, y, zplane)), $(realZ(x, y, zplane))){$FN_Sum[isum]};\n";
      append "AVG_flux_Arplus.txt" << out "$(str12(alpha))  $(str12(zplane))  $(str12(FN_Sum[isum]))\n";
  out "};\n";

After running the script, the normal fluxes extracted for the target surface for each time step can be found as POS file flux_Arplus_XX.XXus_EXT.pos within the ./extract directory. Furthermore, two files, AVG_flux_Arplus.pos and AVG_flux_Arplus.txt, contain the projected flux averaged over the given time interval. The resulting representation in GMSH may look as in the following Fig. 2. Since the points are in regular arrangement, they can be triangulated in GMSH even if they are not located within a plane surface.

Fig. 2: Averaged projected flux on cylindrical target surface

Creating animations

As first step, multiple successive views have to be loaded into gmsh and combined into a sequence with the command Combine Timesteps→From All Views (see Fig. 3). The command From Visible Views would only consider the views with an active checkbox in front. Another possibility to combine plots would be to use one of the write_pos_time_steps() functions from the RVM postprocessing module.

Figure 3: GMSH menu entry for creating image sequences

The next step is to create a GMSH script which subsequently saves each image of the sequence as PNG file. A possible example of such a script is given in the following:

Print.Background       =  1; // Use this, if the background color should be included
General.GraphicsWidth  =  1280;
General.GraphicsHeight =  720;
Steps = View[0].NbTimeStep; // Automatic determination of the number of time steps
For Num In {0:Steps-1}
  View[0].TimeStep = Num;
  Print Sprintf("ani-%02g.png", Num);

With Windows Movie Maker these PNG files can be combined into a WMV film sequence. Another possibility is the use of Blender which is described step-by-step in this tutorial. For creating animated GIFs, the linux tool convert from the ImageMagick package (maybe not installed by default on your cluster) can be used as follows:

convert -delay 20 -loop 0 ani*.png output.gif

With the -delay switch, the time delay between frames can be specified in 1/100 seconds. The switch -loop 0 creates an endless loop.

The POS files are not required in the following operations and can be deleted