From b5ed1429afc941b802029235f71f61ceba74bcc4 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 23 Jun 2010 18:29:18 +0200 Subject: [PATCH] Topological partition almost works --- eplf.config | 1 + scripts/eplf_edit.py | 6 +- scripts/input_wrapper.py | 15 ++++ src/ezfio_interface.irp.f | 1 + src/grid.irp.f | 144 ++++++++++++++++++++++++++++++++++++-- src/main.irp.f | 4 ++ src/to_cube.irp.f | 1 + 7 files changed, 163 insertions(+), 9 deletions(-) diff --git a/eplf.config b/eplf.config index 920c11c..109adf6 100644 --- a/eplf.config +++ b/eplf.config @@ -51,6 +51,7 @@ grid_data density real (grid_num_x,grid_num_y,grid_num_z) 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) compute nproc integer diff --git a/scripts/eplf_edit.py b/scripts/eplf_edit.py index da32b93..3ed86c1 100755 --- a/scripts/eplf_edit.py +++ b/scripts/eplf_edit.py @@ -214,7 +214,7 @@ def write_main_file(file,inp): x = ' ' exec "if inp.%s: x = 'X'"%(p) print >>file, "(%s) %s"%(x,p[8:]) - print >>file, "" + ## Execute temporary input file ############################### @@ -262,21 +262,17 @@ def edit_temp_file(input_file,write_file,read_file,saved_file=None): def build_script(inp): - grid = "${EPLF_PATH}/bin/to_cube %s"%(inp.name) run_script = "run_"+inp.name file = open(run_script,'w') buffer = """ #!/bin/bash source ${HOME}/.eplfrc $RUNCOMMAND -$RUNGRID """.lstrip() command = "${EPLF_PATH}/bin/eplf %s"%(inp.name) if inp.has_mpi: command = "${MPIRUN} -np %d "%(inp.nproc)+command - grid = "${MPIRUN} -np 1 "+command buffer = buffer.replace("$RUNCOMMAND",command) - buffer = buffer.replace("$RUNGRID",grid) print >>file, buffer file.close() os.chmod(run_script,0744) diff --git a/scripts/input_wrapper.py b/scripts/input_wrapper.py index bab210b..3f488f5 100644 --- a/scripts/input_wrapper.py +++ b/scripts/input_wrapper.py @@ -13,6 +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 """.split() rw_data_full = rw_data + geom_data @@ -26,6 +27,9 @@ DEFAULT_COMPUTE_DENSITY_GRAD = False DEFAULT_COMPUTE_EPLF_LAPL = False DEFAULT_COMPUTE_ELF_LAPL = False DEFAULT_COMPUTE_DENSITY_LAPL = False +DEFAULT_COMPUTE_EPLF_PART = False +DEFAULT_COMPUTE_ELF_PART = False +DEFAULT_COMPUTE_DENSITY_PART = False DEFAULT_POINT_NUM = [80,80,80] ###################################################################### @@ -322,6 +326,17 @@ class InputFile(object): raise InputFileExn("Wrong type") ezfio.set_compute_density_lapl(value) + def get_compute_elf_partition(self): + if ezfio.has_compute_elf_partition(): + return ezfio.get_compute_elf_partition() + else: + return DEFAULT_COMPUTE_ELF_PART + + def set_compute_elf_partition(self,value): + if not isinstance(value,bool): + raise InputFileExn("Wrong type") + ezfio.set_compute_elf_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/ezfio_interface.irp.f b/src/ezfio_interface.irp.f index 02adf26..a47fa63 100644 --- a/src/ezfio_interface.irp.f +++ b/src/ezfio_interface.irp.f @@ -40,6 +40,7 @@ data = [ \ ("grid_data_density" , "real" , "(grid_x_num,grid_y_num,grid_z_num)" ), ("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)" ), ] data_no_set = [\ diff --git a/src/grid.irp.f b/src/grid.irp.f index af5f31a..bb56299 100644 --- a/src/grid.irp.f +++ b/src/grid.irp.f @@ -249,13 +249,149 @@ SUBST [ X ] END_TEMPLATE -IRP_IF caca -BEGIN_PROVIDER [ integer, grid_partition_$X, (grid_x_num,grid_y_num,grid_z_num) ] +BEGIN_TEMPLATE + +BEGIN_PROVIDER [ real, grid_$X_partition, (grid_x_num,grid_y_num,grid_z_num) ] implicit none BEGIN_DOC ! Create the topological partition of $X END_DOC -END_PROVIDER -IRP_ENDIF + integer :: ix, iy, iz + real, parameter :: UNDEFINED=123456789.123456789 + + if (mpi_master) then + + grid_$X_partition(1,1,1) = UNDEFINED + call get_grid_data_$X_partition(grid_$X_partition) + if (grid_$X_partition(1,1,1) /= UNDEFINED) then + return + endif + + do iz=1,grid_z_num + do iy=1,grid_y_num + do ix=1,grid_x_num + grid_$X_partition(ix,iy,iz) = 0. + enddo + enddo + enddo + + do iz=1,grid_z_num + do iy=1,grid_y_num + do ix=1,grid_x_num + call find_attribution(ix,iy,iz,grid_$X,grid_$X_partition) + enddo + enddo + enddo + + call set_grid_data_$X_partition(grid_$X_partition) + endif + +END_PROVIDER + +SUBST [ X ] + elf;; + +END_TEMPLATE + +BEGIN_PROVIDER [ real, current_partition_index ] + implicit none + BEGIN_DOC +! Current number for the space partitioning + END_DOC + current_partition_index = 1. +END_PROVIDER + +recursive subroutine find_attribution(ix,iy,iz,g,a) + implicit none + + integer, intent(in) :: ix, iy, iz + real, intent(in) :: g(grid_x_num,grid_y_num,grid_z_num) + real, intent(inout) :: a(grid_x_num,grid_y_num,grid_z_num) + + real :: buffer(-1:1,-1:1,-1:1) + + integer :: i,j,k + + if (a(ix,iy,iz) /= 0) then + return + else if (ix == 1) then + call find_attribution(2,iy,iz,g,a) + a(ix,iy,iz) = a(2,iy,iz) + else if (ix == grid_x_num) then + call find_attribution(ix-1,iy,iz,g,a) + a(ix,iy,iz) = a(ix-1,iy,iz) + else if (iy == 1) then + call find_attribution(ix,2,iz,g,a) + a(ix,iy,iz) = a(ix,2,iz) + else if (iy == grid_y_num) then + call find_attribution(ix,iy-1,iz,g,a) + a(ix,iy,iz) = a(ix,iy-1,iz) + else if (iz == 1) then + call find_attribution(ix,iy,2,g,a) + a(ix,iy,iz) = a(ix,iy,2) + else if (iz == grid_z_num) then + call find_attribution(ix,iy,iz-1,g,a) + a(ix,iy,iz) = a(ix,iy,iz-1) + + else + + ! Copy the environment of the current point + do i=-1,1 + do j=-1,1 + buffer(-1,j,i) = g(ix-1, iy+j, iz+i) + buffer( 0,j,i) = g(ix , iy+j, iz+i) + buffer( 1,j,i) = g(ix+1, iy+j, iz+i) + enddo + enddo + + ! Avoid strict equality between points + buffer(0,0,0) *= 1.000001 + + ! Find the index of the maximum value + real :: m,mm + integer :: im(3) + mm = buffer(-1,-1,-1) + m = buffer(-1,-1,-1) + do i=-1,1 + do j=-1,1 + do k=-1,1 + if (buffer(k,j,i) > m) then + m = buffer(k,j,i) + im(1) = k + im(2) = j + im(3) = i + endif + mm = min(buffer(k,j,i),mm) + enddo + enddo + enddo + + + 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 + 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 + +end diff --git a/src/main.irp.f b/src/main.irp.f index c31ad0c..6d48700 100644 --- a/src/main.irp.f +++ b/src/main.irp.f @@ -17,6 +17,10 @@ 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 diff --git a/src/to_cube.irp.f b/src/to_cube.irp.f index 7dedc1a..dd33184 100644 --- a/src/to_cube.irp.f +++ b/src/to_cube.irp.f @@ -52,6 +52,7 @@ program to_cube density_grad;4;,4;; elf_grad;4;,4;; eplf_grad;4;,4;; + elf_partition;3; ;; END_TEMPLATE 10 format (2X,I3,3(2X,F10.6)) 11 format (2X,I3,4(2X,F10.6))