domain/decomposition generators

matlab/generate_grid.m

@@ -0,0 +1,112 @@
+function [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
+          xu1,yu1,zu1,xv1,yv1,zv1,xw1,yw1,zw1,xp1,yp1,zp1,dx,dy,dz]=...
+            generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,x_periodic,y_periodic,z_periodic)
+  
+  if x_periodic==1
+    nxu=nxp;
+    nxv=nxp;
+    nxw=nxp;
+  else
+    nxu=nxp+1;
+    nxv=nxp;
+    nxw=nxp;
+  end
+  if y_periodic==1
+    nyu=nyp;
+    nyv=nyp;
+    nyw=nyp;
+  else
+    nyu=nyp;
+    nyv=nyp+1;
+    nyw=nyp;
+  end
+  if z_periodic==1
+    nzu=nzp;
+    nzv=nzp;
+    nzw=nzp;
+  else
+    nzu=nzp;
+    nzv=nzp;
+    nzw=nzp+1;
+  end
+  %grid step:
+  if x_periodic==1
+    dx=(b-a)/(nxp);
+  else
+    dx=(b-a)/(nxp-1);
+  end
+  if y_periodic==1
+    dy=(d-c)/(nyp);
+  else
+    dy=(d-c)/(nyp-1);
+  end
+  if z_periodic==1
+    dz=(f-e)/(nzp);
+  else
+    dz=(f-e)/(nzp-1);
+  end
+  %
+  if x_periodic==1
+    xu1=a+((1:nxu)-1)*dx;
+  else
+    xu1=a+((1:nxu)-3/2)*dx;
+  end
+  if y_periodic==1
+    yu1=c+((1:nyu)-1/2)*dy;
+  else
+    yu1=c+((1:nyu)-1)*dy;
+  end
+  if z_periodic==1
+    zu1=e+((1:nzu)-1/2)*dz;
+  else
+    zu1=e+((1:nzu)-1)*dz;
+  end
+  %
+  if x_periodic==1
+    xv1=a+((1:nxv)-1/2)*dx;
+  else
+    xv1=a+((1:nxv)-1)*dx;
+  end
+  if y_periodic==1
+    yv1=c+((1:nyv)-1)*dy;
+  else
+    yv1=c+((1:nyv)-3/2)*dy;
+  end
+  if z_periodic==1
+    zv1=e+((1:nzv)-1/2)*dz;
+  else
+    zv1=e+((1:nzv)-1)*dz;
+  end
+  %
+  if x_periodic==1
+    xw1=a+((1:nxw)-1/2)*dx;
+  else
+    xw1=a+((1:nxw)-1)*dx;
+  end
+  if y_periodic==1
+    yw1=c+((1:nyw)-1/2)*dy;
+  else
+    yw1=c+((1:nyw)-1)*dy;
+  end
+  if z_periodic==1
+    zw1=e+((1:nzw)-1)*dz;
+  else
+    zw1=e+((1:nzw)-3/2)*dz;
+  end
+  %
+  if x_periodic==1
+    xp1=a+((1:nxp)-1/2)*dx;
+  else
+    xp1=a+((1:nxp)-1)*dx;
+  end
+  if y_periodic==1
+    yp1=c+((1:nyp)-1/2)*dy;
+  else
+    yp1=c+((1:nyp)-1)*dy;
+  end
+  if z_periodic==1
+    zp1=e+((1:nzp)-1/2)*dz;
+  else
+    zp1=e+((1:nzp)-1)*dz;
+  end
+end

matlab/generate_procgrid.m

@@ -0,0 +1,51 @@
+function [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
+          ibegv,iendv,jbegv,jendv,kbegv,kendv,...
+          ibegw,iendw,jbegw,jendw,kbegw,kendw,...
+          ibegp,iendp,jbegp,jendp,kbegp,kendp]=generate_procgrid(...
+              nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,...
+              nxp,nyp,nzp,nxprocs,nyprocs,nzprocs)
+  %
+  % determine the 3d decomposition (pointers) of fields for u,v,w,p
+  %
+  chi='_';
+  % U:
+  [ibegu,iendu]=mpe_decomp1d(nxu,nxprocs,chi);
+  [jbegu,jendu]=mpe_decomp1d(nyu,nyprocs,chi);
+  [kbegu,kendu]=mpe_decomp1d(nzu,nzprocs,chi);
+  % V:
+  [ibegv,iendv]=mpe_decomp1d(nxv,nxprocs,chi);
+  [jbegv,jendv]=mpe_decomp1d(nyv,nyprocs,chi);
+  [kbegv,kendv]=mpe_decomp1d(nzv,nzprocs,chi);
+  % W:
+  [ibegw,iendw]=mpe_decomp1d(nxw,nxprocs,chi);
+  [jbegw,jendw]=mpe_decomp1d(nyw,nyprocs,chi);
+  [kbegw,kendw]=mpe_decomp1d(nzw,nzprocs,chi);
+  % P:
+  [ibegp,iendp]=mpe_decomp1d(nxp,nxprocs,chi);
+  [jbegp,jendp]=mpe_decomp1d(nyp,nyprocs,chi);
+  [kbegp,kendp]=mpe_decomp1d(nzp,nzprocs,chi);
+end
+
+function [s,e]=mpe_decomp1d(n,numprocs,chi)
+  %
+  % determines 1d decomposition (as in fortran dns code)
+  %
+  for myid=0:numprocs-1
+    nlocal  = floor(n/ numprocs);
+    s(myid+1) = myid * nlocal + 1;
+    deficit = mod(n,numprocs);
+    s(myid+1) = s(myid+1) + min(myid,deficit);
+    if myid<deficit
+      nlocal = nlocal + 1;
+    end
+    e(myid+1) = s(myid+1) + nlocal - 1;
+    if (e(myid+1)> n | myid==numprocs-1)
+      e(myid+1) = n;
+    end
+    if chi=='p'
+      %    /* finally: augment all counts by one */
+      s(myid+1) = s(myid+1) + 1;
+      e(myid+1) = e(myid+1) + 1;
+    end
+  end
+end