Merge pull request #258 from OpenSEMBA/MPI-problem-with-MTLN

[WIP] Investigating possible changes to MTL MPI . Problems with cables along the frontier between MPI layers

Author: Alberto-o

CMakeLists.txt

@@ -28,6 +28,11 @@ option(SEMBA_FDTD_MAIN_LIB "Compiles main library" ON)
 option(SEMBA_FDTD_COMPONENTS_LIB "Compiles components library" ON)
 option(SEMBA_FDTD_OUTPUTS_LIB "Compiles outputs library" ON)
 # Compilation defines.
+if(CMAKE_BUILD_TYPE MATCHES "Release" OR CMAKE_BUILD_TYPE MATCHES "release" )
+	add_definitions(-DCompileWithRelease)
+else()
+	add_definitions(-DCompileWithDebug)
+endif()
 if(SEMBA_FDTD_ENABLE_SMBJSON)
 	add_definitions(-DCompileWithSMBJSON)
 endif()

src_main_pub/observation.F90

@@ -4038,6 +4038,9 @@ subroutine FlushMTLNObservationFiles(nEntradaRoot, mtlnProblem)
       unit = 2000
       do i = 1, size(mtln_solver%bundles)
         do j = 1, size(mtln_solver%bundles(i)%probes)
+#ifdef CompileWithMPI
+          if (.not. mtln_solver%bundles(i)%probes(j)%in_layer) cycle
+#endif          
           path = trim(trim(nEntradaRoot)//"_"//trim(mtln_solver%bundles(i)%probes(j)%name)//".dat")
           open (unit=unit, file=trim(path))
           write (*, *) 'name: ', trim(mtln_solver%bundles(i)%probes(j)%name)

src_main_pub/timestepping.F90

@@ -1139,14 +1139,14 @@ end subroutine initializeLumped
 
       subroutine initializeWires()
          real (kind=rkind_tiempo) :: dtcritico, newdtcritico
-         ! real (kind=rkind) :: dtcritico, newdtcritico
          character(len=BUFSIZE) :: dubuf, buff
          logical :: l_auxinput, l_auxoutput
 #ifdef CompileWithMPI
          integer(kind=4) :: ierr
 #endif
 
          dtcritico=this%sgg%dt
+#ifndef CompileWithMTLN         
          if ((trim(adjustl(this%control%wiresflavor))=='holland') .or. &
             (trim(adjustl(this%control%wiresflavor))=='transition')) then
 #ifdef CompileWithMPI
@@ -1234,10 +1234,21 @@ subroutine initializeWires()
             endif
          endif
 #endif
+
+
+#else 
+! else of #ifndef CompileWithMTLN          
+#ifdef CompileWithMPI
+         call MPI_Barrier(SUBCOMM_MPI,ierr)
+#endif
+         write(dubuf,*) 'Init MTLN Wires...';  call print11(this%control%layoutnumber,dubuf)
+         call InitWires_mtln(this%sgg,Ex,Ey,Ez,this%eps0, this%mu0, this%mtln_parsed,this%thereAre%MTLNbundles, dtcritico)
+#endif
+
+
       !!!sincroniza el dtcritico
 #ifdef CompileWithMPI
          newdtcritico = 0.0
-         ! call MPI_AllReduce( dtcritico, newdtcritico, 1_4, REALSIZE, MPI_MIN, SUBCOMM_MPI, ierr)
          call MPI_AllReduce( dtcritico, newdtcritico, 1_4, REALSIZE_tiempo, MPI_MIN, SUBCOMM_MPI, ierr)
          dtcritico=newdtcritico
 #endif
@@ -1246,7 +1257,11 @@ subroutine initializeWires()
             if ((this%control%layoutnumber==0).and.this%control%verbose) call WarnErrReport(buff)
          else
             if (.not.(this%control%resume.and.this%control%permitscaling)) then !no abortasr solo advertir si permittivity scaling
+#ifdef CompileWithMTLN
+               write(buff,'(a,e10.2e3)')  'WIR_ERROR: Possibly UNSTABLE dt, make dt < ',dtcritico
+#else
                write(buff,'(a,e10.2e3)')  'WIR_ERROR: Possibly UNSTABLE dt, decrease wire radius, number of parallel WIREs, use -stableradholland or make dt < ',dtcritico
+#endif
                if (this%control%layoutnumber==0) call WarnErrReport(buff,.true.)
             else
                write(buff,'(a,e10.2e3)')  'WIR_WARNING: Resume and Pscaling with wires. Possibly UNSTABLE dt, decrease wire radius, number of parallel WIREs: dt is over ',dtcritico
@@ -1255,15 +1270,6 @@ subroutine initializeWires()
          endif
       !!!
 !!
-#ifdef CompileWithMTLN
-#ifdef CompileWithMPI
-         call MPI_Barrier(SUBCOMM_MPI,ierr)
-#endif
-         write(dubuf,*) 'Init MTLN Wires...';  call print11(this%control%layoutnumber,dubuf)
-         call InitWires_mtln(this%sgg,Ex,Ey,Ez,Idxh,Idyh,Idzh,this%eps0, this%mu0, this%mtln_parsed,this%thereAre%MTLNbundles)
-#else
-         write(buff,'(a)') 'WIR_ERROR: Executable was not compiled with MTLN modules.'
-#endif
 
       end subroutine initializeWires
 

src_mtln/mtl.F90

@@ -14,10 +14,14 @@ module mtl_mod
     integer, parameter :: COMM_SEND = 1
     integer, parameter :: COMM_RECV = -1
     integer, parameter :: COMM_NONE = 0
+    integer, parameter :: COMM_FIELD = 1
+    integer, parameter :: COMM_V = 2
+    integer, parameter :: COMM_BOTH = 3
 
     type, public :: communicator_t
         integer :: field_index = -1, v_index = -1
         integer :: comm_task = COMM_NONE
+        integer :: comm_type = COMM_NONE
         integer :: delta_rank = 0
     end type
     type, public :: comm_t
@@ -180,7 +184,7 @@ function mtl_unshielded(lpul, cpul, rpul, gpul, &
             call res%initStepSizeAndFieldSegments(step_size, segments, layer_indices)
             call res%initCommunicators(alloc_z)
             res%layer_indices = layer_indices
-            res%bundle_in_layer = bundle_in_layer
+            res%bundle_in_layer = bundle_in_layer 
         else
             res%step_size =  step_size
             allocate(res%layer_indices(0,0))
@@ -350,6 +354,7 @@ subroutine initStepSizeAndFieldSegments(this, step_size, segments, layer_indices
             if (j /= size(layer_indices,1)) then
                 this%step_size(n + layer_indices(j,2) - layer_indices(j,1) + 1) = this%step_size(n + layer_indices(j,2) - layer_indices(j,1))
                 this%segments(n + layer_indices(j,2) - layer_indices(j,1) + 1) = this%segments(n + layer_indices(j,2) - layer_indices(j,1))
+                this%segments(n + layer_indices(j,2) - layer_indices(j,1) + 1)%orientation = -1
                 n = n + 1
             end if
             n = n + layer_indices(j,2) - layer_indices(j,1) + 1
@@ -361,7 +366,7 @@ subroutine initStepSizeAndFieldSegments(this, step_size, segments, layer_indices
     subroutine initCommunicators(this, alloc_z)
         class(mtl_t) :: this
         integer (kind =4), dimension(2) :: alloc_z
-        integer :: j, n
+        integer :: j, n, z
         integer :: rank, ierr
         integer (kind =4) :: z_init, z_end
         type(communicator_t), dimension(:), allocatable :: aux_comm
@@ -374,7 +379,33 @@ subroutine initCommunicators(this, alloc_z)
         z_end = alloc_z(2)
 
         do j = 1, size(this%segments)
+            if (this%segments(j)%orientation == -1) cycle
 
+            z = this%segments(j)%z
+            if (.not. isSegmentZOriented(j) .and. ((z == z_end) .or. (z == z_init + 1))) then 
+
+                n = size(this%mpi_comm%comms)
+                deallocate(aux_comm)
+                allocate(aux_comm(n+1))
+                aux_comm(1:n) = this%mpi_comm%comms
+                
+                aux_comm(n+1)%field_index = j
+                aux_comm(n+1)%comm_type = COMM_FIELD
+                aux_comm(n+1)%v_index = -1
+                if (z == z_end) then 
+                    aux_comm(n+1)%delta_rank = 1
+                    aux_comm(n+1)%comm_task = COMM_RECV
+                else if (z == z_init + 1) then 
+                    aux_comm(n+1)%delta_rank = -1
+                    aux_comm(n+1)%comm_task = COMM_SEND
+                end if
+
+                deallocate(this%mpi_comm%comms)
+                allocate(this%mpi_comm%comms(n+1))
+                this%mpi_comm%comms = aux_comm
+
+
+            end if
             if (isSegmentZOriented(j) .and. &
                (isSegmentNextToLayerEnd(j,z_end) .or. isSegmentNextToLayerInit(j,z_init))) then
 
@@ -383,6 +414,7 @@ subroutine initCommunicators(this, alloc_z)
                 allocate(aux_comm(n+1))
                 aux_comm(1:n) = this%mpi_comm%comms
                 aux_comm(n+1)%field_index = j
+                aux_comm(n+1)%comm_type = COMM_BOTH
 
                 if (isSegmentNextToLayerEnd(j,z_end)) then 
                     aux_comm(n+1)%delta_rank = 1
@@ -438,7 +470,12 @@ logical function isSegmentZOriented(j)
 
     logical function isSegmentZPositive(j)
         integer, intent(in) :: j
-        isSegmentZPositive = (this%segments(j)%orientation > 0)
+        isSegmentZPositive = (this%segments(j)%orientation == 3)
+    end function
+
+    logical function isSegmentZNegative(j)
+        integer, intent(in) :: j
+        isSegmentZNegative = (this%segments(j)%orientation == -3)
     end function
 
     logical function isSegmentBeforeLayerEnd(j, z_end)
@@ -473,14 +510,14 @@ logical function isSegmentNextToLayerEnd(j, z_end)
         integer, intent(in) :: j, z_end
         integer :: z
         z = this%segments(j)%z
-        isSegmentNextToLayerEnd = (abs(z-z_end)<= 1)
+        isSegmentNextToLayerEnd = (z==z_end) .or. (z==z_end-1)
     end function
 
     logical function isSegmentNextToLayerInit(j, z_init)
         integer, intent(in) :: j, z_init
         integer :: z
         z = this%segments(j)%z
-        isSegmentNextToLayerInit = (abs(z-z_init-1) <= 1)
+        isSegmentNextToLayerInit = (z==z_init) .or. (z==z_init+1)
     end function
 
 

src_mtln/mtl_bundle.F90

@@ -403,6 +403,7 @@ subroutine Comm_MPI_V(this)
         call MPI_COMM_RANK(SUBCOMM_MPI, rank, ierr)
         number_of_conductors = size(this%v,1)
         do i = 1, size(this%mpi_comm%comms)
+            if (this%mpi_comm%comms(i)%comm_type == COMM_V .or. this%mpi_comm%comms(i)%comm_type == COMM_BOTH) then 
                 if (this%mpi_comm%comms(i)%comm_task == COMM_SEND) then 
                     do c = 1, number_of_conductors
                         call MPI_send(this%v(c, this%mpi_comm%comms(i)%v_index),1, REALSIZE, & 
@@ -419,6 +420,7 @@ subroutine Comm_MPI_V(this)
                                       SUBCOMM_MPI, status, ierr)
                     end do
                 end if
+            end if
         end do
 
 
@@ -429,17 +431,19 @@ subroutine Comm_MPI_Fields(this)
         integer :: i, ierr, rank, status(MPI_STATUS_SIZE)
         call MPI_COMM_RANK(SUBCOMM_MPI, rank, ierr)
         do i = 1, size(this%mpi_comm%comms)
-            if (this%mpi_comm%comms(i)%comm_task == COMM_SEND) then 
-                call MPI_send(this%external_field_segments(this%mpi_comm%comms(i)%field_index)%field, 1, REALSIZE, & 
-                              rank+this%mpi_comm%comms(i)%delta_rank, & 
-                              100*(rank+this%mpi_comm%comms(i)%delta_rank+1), & 
-                              SUBCOMM_MPI, ierr)
-            end if
-            if (this%mpi_comm%comms(i)%comm_task == COMM_RECV) then 
-                call MPI_recv(this%external_field_segments(this%mpi_comm%comms(i)%field_index)%field,1, REALSIZE, & 
-                              rank+this%mpi_comm%comms(i)%delta_rank, &
-                              100*(rank+1), &
-                              SUBCOMM_MPI, status, ierr)
+            if (this%mpi_comm%comms(i)%comm_type == COMM_FIELD .or. this%mpi_comm%comms(i)%comm_type == COMM_BOTH) then 
+                if (this%mpi_comm%comms(i)%comm_task == COMM_SEND) then 
+                    call MPI_send(this%external_field_segments(this%mpi_comm%comms(i)%field_index)%field, 1, REALSIZE, & 
+                                    rank+this%mpi_comm%comms(i)%delta_rank, & 
+                                    100*(rank+this%mpi_comm%comms(i)%delta_rank+1), & 
+                                    SUBCOMM_MPI, ierr)
+                end if
+                if (this%mpi_comm%comms(i)%comm_task == COMM_RECV) then 
+                    call MPI_recv(this%external_field_segments(this%mpi_comm%comms(i)%field_index)%field,1, REALSIZE, & 
+                                    rank+this%mpi_comm%comms(i)%delta_rank, &
+                                    100*(rank+1), &
+                                    SUBCOMM_MPI, status, ierr)
+                end if
             end if
         end do
     end subroutine

src_mtln/mtln_solver.F90

@@ -16,7 +16,6 @@ module mtln_solver_mod
         type(network_manager_t) :: network_manager
         type(probe_t), allocatable, dimension(:) :: probes
         integer :: number_of_bundles
-        ! logical :: has_multiwires
         integer :: number_of_steps
     contains
 
@@ -70,8 +69,6 @@ function mtlnCtor(parsed, alloc) result(res)
             return
         end if
 
-        ! res%has_multiwires = parsed%has_multiwires
-        
         res%dt = pre%dt
         res%time  = 0.0
         res%final_time = pre%final_time
@@ -126,6 +123,10 @@ subroutine step_alone(this)
     subroutine setExternalLongitudinalField(this)
         class(mtln_t) :: this
         integer :: i
+#ifdef CompileWithMPI
+        integer :: ierr
+      call MPI_Barrier(SUBCOMM_MPI,ierr)
+#endif
         do i = 1, this%number_of_bundles
             if (this%bundles(i)%bundle_in_layer) call this%bundles(i)%setExternalLongitudinalField()
         end do
@@ -150,10 +151,10 @@ subroutine advanceNWVoltage(this)
         integer :: i,j
         integer ::b, c, v_idx, i_idx
         integer :: n
-#ifdef CompileWithMPI
-        integer (kind=4) :: ierr
-        call mpi_barrier(subcomm_mpi, ierr)
-#endif
+! #ifdef CompileWithMPI
+!         integer (kind=4) :: ierr
+!         call mpi_barrier(subcomm_mpi, ierr)
+! #endif
         if (this%number_of_bundles /= 0) then 
             do i = 1, size(this%network_manager%networks)
                 do j = 1, size(this%network_manager%networks(i)%nodes)
@@ -191,6 +192,10 @@ subroutine advanceNWVoltage(this)
     subroutine advanceBundlesCurrent(this)
         class(mtln_t) :: this
         integer :: i
+#ifdef CompileWithMPI
+        integer :: ierr
+        call mpi_barrier(subcomm_mpi, ierr)
+#endif
         do i = 1, this%number_of_bundles
             if (this%bundles(i)%bundle_in_layer) call this%bundles(i)%advanceCurrent()
 

src_mtln/network_manager.F90

@@ -83,13 +83,16 @@ function network_managerCtor(networks, description, final_time, dt) result(res)
         character(*), dimension(:), intent(in) :: description
         real, intent(in) :: final_time, dt
         type(network_manager_t) :: res
-
+        logical :: printInput = .true.
         res%dt = dt
         res%time = 0.0
         res%networks = networks
         call res%circuit%init(copy_node_names(networks), copy_sources(networks))
         res%circuit%dt = dt
-        call res%circuit%readInput(description, .true.)
+#ifdef CompileWithRelease
+        printInput = .false.
+#endif        
+        call res%circuit%readInput(description, printInput)
         call res%circuit%setModStopTimes(dt)
 
     end function