From 3ef75ef8b8b3f05a942a7f52d4778c9f136da7a5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 24 Jun 2010 10:40:51 +0200 Subject: [PATCH] Improved topological partition. Added eplf and density partition --- eplf.config | 4 +++ scripts/eplf_edit.py | 20 ++++++++++++++ scripts/input_wrapper.py | 24 +++++++++++++++- src/compute.irp.f | 2 ++ src/elf_function.irp.f | 2 +- src/ezfio_interface.irp.f | 4 +++ src/grid.irp.f | 58 +++++++++++++++++++++++---------------- src/main.irp.f | 17 +++++++++--- src/to_cube.irp.f | 2 ++ 9 files changed, 103 insertions(+), 30 deletions(-) diff --git a/eplf.config b/eplf.config index 109adf6..33c656f 100644 --- a/eplf.config +++ b/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 diff --git a/scripts/eplf_edit.py b/scripts/eplf_edit.py index 3ed86c1..d6df965 100755 --- a/scripts/eplf_edit.py +++ b/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 diff --git a/scripts/input_wrapper.py b/scripts/input_wrapper.py index 3f488f5..90ce3ff 100644 --- a/scripts/input_wrapper.py +++ b/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) diff --git a/src/compute.irp.f b/src/compute.irp.f index 0d23769..7e32ec4 100644 --- a/src/compute.irp.f +++ b/src/compute.irp.f @@ -20,4 +20,6 @@ SUBST [ X ] density_grad;; density_lapl;; elf_partition;; + eplf_partition;; + density_partition;; END_TEMPLATE diff --git a/src/elf_function.irp.f b/src/elf_function.irp.f index 3e9b74f..601ada0 100644 --- a/src/elf_function.irp.f +++ b/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 diff --git a/src/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index a47fa63..c6305f8 100644 --- a/src/ezfio_interface.irp.f +++ b/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 = [\ diff --git a/src/grid.irp.f b/src/grid.irp.f index bb56299..0c79bb8 100644 --- a/src/grid.irp.f +++ b/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,33 +375,30 @@ 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 - current_partition_index += 1. - TOUCH current_partition_index + if (a(ix,iy,iz) == 0.) then + current_partition_index += 1. + 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 + a(ix,iy,iz) = a(ix+im(1), iy+im(2), iz+im(3)) endif endif diff --git a/src/main.irp.f b/src/main.irp.f index 6d48700..4b7fed8 100644 --- a/src/main.irp.f +++ b/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 diff --git a/src/to_cube.irp.f b/src/to_cube.irp.f index dd33184..effb9ae 100644 --- a/src/to_cube.irp.f +++ b/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))