Improved topological partition. Added eplf and density partition

eplf.config

@@ -52,6 +52,8 @@ grid_data
  density_grad        real (grid_num_x,grid_num_y,grid_num_z,4)
  density_lapl        real (grid_num_x,grid_num_y,grid_num_z)
  elf_partition       real (grid_num_x,grid_num_y,grid_num_z)
+ eplf_partition      real (grid_num_x,grid_num_y,grid_num_z)
+ density_partition   real (grid_num_x,grid_num_y,grid_num_z)
 
 compute
  nproc               integer
@@ -65,4 +67,6 @@ compute
  density_grad        logical
  density_lapl        logical
  elf_partition       logical
+ eplf_partition      logical
+ density_partition   logical
 

scripts/eplf_edit.py

@@ -165,18 +165,38 @@ class GridEditor(Editor):
     buffer = lines.pop(0).lower().split()
     if buffer[0] != "default":
       self.inp.point_num = map(int,buffer)
+    else:
+      try:
+       os.remove("%s/grid/point_num.gz"%self.inp.name)
+      except:
+       pass
 
     buffer = lines.pop(0).lower().split()
     if buffer[0] != "default":
       self.inp.origin = map(float,buffer)
+    else:
+      try:
+       os.remove("%s/grid/origin.gz"%self.inp.name)
+      except:
+       pass
 
     buffer = lines.pop(0).lower().split()
     if buffer[0] != "default":
       self.inp.opposite = map(float,buffer)
+    else:
+      try:
+       os.remove("%s/grid/opposite.gz"%self.inp.name)
+      except:
+       pass
 
     buffer = lines.pop(0).lower().split()
     if buffer[0] != "default":
       self.inp.step_size = self.step_size
+    else:
+      try:
+       os.remove("%s/grid/step_size.gz"%self.inp.name)
+      except:
+       pass
 
       
 

scripts/input_wrapper.py

@@ -13,7 +13,7 @@ point_num step_size origin opposite nproc
 compute_elf compute_eplf compute_density
 compute_elf_grad compute_eplf_grad compute_density_grad
 compute_elf_lapl compute_eplf_lapl compute_density_lapl
-compute_elf_partition
+compute_elf_partition compute_eplf_partition compute_density_partition
 """.split()
 rw_data_full = rw_data + geom_data
 
@@ -337,6 +337,28 @@ class InputFile(object):
       raise InputFileExn("Wrong type")
     ezfio.set_compute_elf_partition(value)
 
+  def get_compute_eplf_partition(self):
+    if ezfio.has_compute_eplf_partition():
+      return ezfio.get_compute_eplf_partition()
+    else:
+      return DEFAULT_COMPUTE_ELF_PART
+
+  def set_compute_eplf_partition(self,value):
+    if not isinstance(value,bool):
+      raise InputFileExn("Wrong type")
+    ezfio.set_compute_eplf_partition(value)
+
+  def get_compute_density_partition(self):
+    if ezfio.has_compute_density_partition():
+      return ezfio.get_compute_density_partition()
+    else:
+      return DEFAULT_COMPUTE_DENSITY_PART
+
+  def set_compute_density_partition(self,value):
+    if not isinstance(value,bool):
+      raise InputFileExn("Wrong type")
+    ezfio.set_compute_density_partition(value)
+
   # Build the corresponding properties
   for i in ro_data:
     exec "%s = property(fget=get_%s,fset=None)"%(i,i)

src/compute.irp.f

@@ -20,4 +20,6 @@ SUBST [ X ]
  density_grad;;
  density_lapl;;
  elf_partition;;
+ eplf_partition;;
+ density_partition;;
 END_TEMPLATE

src/elf_function.irp.f

@@ -62,7 +62,7 @@ BEGIN_PROVIDER [ double precision, elf_value_p ]
  Ds = kinetic_energy_p - 0.125d0 *  &
   ( (modulus_a2/density_alpha_value_p) + &
     (modulus_b2/density_beta_value_p) )
- D0 = 2.d0*Cf*density_value_p**(5./3.)
+ D0 = (2.d0*Cf*density_value_p**(5./3.))
  elf_value_p = 1.d0/ (1.d0 + (Ds/D0)**2)
 
 END_PROVIDER

src/ezfio_interface.irp.f

@@ -31,6 +31,8 @@ data = [ \
 ("compute_density_grad"            , "logical"       , ""                      ),
 ("compute_density_lapl"            , "logical"       , ""                      ),
 ("compute_elf_partition"           , "logical"       , ""                      ),
+("compute_density_partition"       , "logical"       , ""                      ),
+("compute_eplf_partition"          , "logical"       , ""                      ),
 ("grid_data_eplf"                  , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
 ("grid_data_eplf_grad"             , "real"          , "(grid_x_num,grid_y_num,grid_z_num,4)" ),
 ("grid_data_eplf_lapl"             , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
@@ -41,6 +43,8 @@ data = [ \
 ("grid_data_density_grad"          , "real"          , "(grid_x_num,grid_y_num,grid_z_num,4)" ),
 ("grid_data_density_lapl"          , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
 ("grid_data_elf_partition"         , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
+("grid_data_eplf_partition"        , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
+("grid_data_density_partition"     , "real"          , "(grid_x_num,grid_y_num,grid_z_num)" ),
 ]
 
 data_no_set = [\

src/grid.irp.f

@@ -97,6 +97,10 @@ BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
   enddo
  enddo
 
+ if (mpi_master) then
+   print *, 'Computing $X' 
+ endif
+
  integer :: icount
  icount = mpi_size
  do iz=1,grid_z_num
@@ -181,6 +185,10 @@ END_PROVIDER
   enddo
  enddo
 
+ if (mpi_master) then
+   print *, 'Computing $X' 
+ endif
+
  integer :: icount
  icount = mpi_size
  do iz=1,grid_z_num
@@ -286,11 +294,14 @@ BEGIN_PROVIDER [ real, grid_$X_partition, (grid_x_num,grid_y_num,grid_z_num) ]
 
    call set_grid_data_$X_partition(grid_$X_partition)
  endif
+ print *,  int(current_partition_index), ' basins found'
 
 END_PROVIDER
 
 SUBST [ X ]
  elf;;
+ eplf;;
+ density;;
 
 END_TEMPLATE
 
@@ -299,7 +310,7 @@ BEGIN_PROVIDER [ real, current_partition_index ]
  BEGIN_DOC  
 ! Current number for the space partitioning
  END_DOC
- current_partition_index = 1.
+ current_partition_index = 0.
 END_PROVIDER
 
 recursive subroutine find_attribution(ix,iy,iz,g,a)
@@ -313,7 +324,7 @@ recursive subroutine find_attribution(ix,iy,iz,g,a)
 
  integer :: i,j,k
 
- if (a(ix,iy,iz) /= 0) then
+ if (a(ix,iy,iz) /= 0.) then
   return
  else if (ix == 1) then
   call find_attribution(2,iy,iz,g,a)
@@ -345,14 +356,16 @@ recursive subroutine find_attribution(ix,iy,iz,g,a)
     enddo
    enddo 
 
-   ! Avoid strict equality between points
-   buffer(0,0,0) *= 1.000001
-
    ! Find the index of the maximum value
-   real :: m,mm
+   real :: m
    integer :: im(3)
-   mm = buffer(-1,-1,-1)
-   m  = buffer(-1,-1,-1)
+   m  = buffer(0,0,0)
+   real :: average, variance
+   im(1) = 0
+   im(2) = 0
+   im(3) = 0
+   average = 0.
+   variance = 0.
    do i=-1,1
     do j=-1,1
      do k=-1,1
@@ -362,34 +375,31 @@ recursive subroutine find_attribution(ix,iy,iz,g,a)
        im(2) = j
        im(3) = i
       endif
-      mm = min(buffer(k,j,i),mm)
+      average += buffer(k,j,i)
+      variance += buffer(k,j,i)**2
      enddo
     enddo
    enddo 
+   variance = variance / 27.
+   average = average / 27.
+   variance = (variance - average**2)
+
+   ! Avoid strict equality between points
+   buffer(0,0,0) += sqrt(variance/26.)
 
 
    if ( (im(1) == 0).and.(im(2) == 0).and.(im(3) == 0) ) then
      ! If the maximum is at the center, a new maximum is found
-     if ( (m-mm)/(m+mm+1.e-12) > 0.01) then
-       do i=-1,1
-        do j=-1,1
-         do k=-1,1
-          if (g(ix+k,iy+j,iz+i) == g(ix,iy,iz)) then
-            a(ix+k,iy+j,iz+i) = current_partition_index
-          endif
-         enddo
-        enddo
-       enddo
+     if (a(ix,iy,iz) == 0.) then
          current_partition_index += 1.
-       TOUCH current_partition_index 
+         a(ix,iy,iz) = current_partition_index
+         print '(3I3,2(2X,F10.7))', ix,iy,iz,average,sqrt(variance)
      endif
    else
      ! Otherwise, get the partition index of the max value
      call find_attribution(ix+im(1), iy+im(2), iz+im(3),g,a)
-     if (a(ix,iy,iz) == 0) then
      a(ix,iy,iz) = a(ix+im(1), iy+im(2), iz+im(3))
    endif
-   endif
 
  endif
 

src/main.irp.f

@@ -17,10 +17,6 @@ program eplf
    PROVIDE grid_elf
    FREE grid_elf
  endif
- if (comp_elf_partition) then
-   PROVIDE grid_elf_partition
-   FREE grid_elf_partition
- endif
  if (comp_elf_grad.or.comp_elf_lapl) then
    PROVIDE grid_elf_grad
    PROVIDE grid_elf_lapl
@@ -38,6 +34,19 @@ program eplf
    FREE grid_density_grad
    FREE grid_density_lapl
  endif
+
+ if (comp_elf_partition) then
+   PROVIDE grid_elf_partition
+   FREE grid_elf_partition
+ endif
+ if (comp_eplf_partition) then
+   PROVIDE grid_eplf_partition
+   FREE grid_eplf_partition
+ endif
+ if (comp_density_partition) then
+   PROVIDE grid_density_partition
+   FREE grid_density_partition
+ endif
  call finish()
 end
 

src/to_cube.irp.f

@@ -53,6 +53,8 @@ program to_cube
    elf_grad;4;,4;;
    eplf_grad;4;,4;;
    elf_partition;3; ;;
+   eplf_partition;3; ;;
+   density_partition;3; ;;
   END_TEMPLATE
   10   format (2X,I3,3(2X,F10.6))
   11   format (2X,I3,4(2X,F10.6))