Initial updates to create_global_domain_structure.F90 for nested overshields

Chris Smartt
1 parent 0cd31073
Showing 15 changed files with 221 additions and 31 deletions Show diff stats
DOCUMENTATION/USER_GUIDE/Imgs/edge_coupled_stripline-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_bundle_form_01-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_bundle_form_01a-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_bundle_form_01b-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_bundle_form_02-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_bundle_form_02a-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_bundle_form_02b-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_cable_form_01-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_cable_form_02-eps-converted-to.pdf
DOCUMENTATION/USER_GUIDE/Imgs/flex_cable_form_03-eps-converted-to.pdf
SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90
TEST_CASES/OVERSHIELD_COAX_FREE_SPACE_OVER_GROUND_PLANE/coax_over_ground_plane.bundle_spec
TEST_CASES/OVERSHIELD_COAX_FREE_SPACE_OVER_GROUND_PLANE/coax_over_ground_plane_ac.spice_model_spec
TEST_CASES/generate_spice_cable_bundle_model
clean_project
@@ -73,8 +73,10 @@
 !                    so that we can work out the is_shield flag properly in all circumstances.
 !     18/10/2017 CJS: include 8b. Copy the cable based conductor labels to the bundle structure
 !     16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions
-!     16/3/2018 CJS add y offset for ML_flex_cable
+!     16/3/2018 CJS  add y offset for ML_flex_cable
 !     22/10/2019 CJS Fix bug with numbering here: overshields were numbered using terminal conductor number not global domain number
+!     24/10/2019 CJS Allow mutliple layers of overshields (previously overshields could not be 'nested' 
+!                    (though they could cover shielded cables).
 !
  
 SUBROUTINE create_global_domain_structure(bundle)
@@ -97,7 +99,7 @@ USE maths
 ! local variables
  
 ! global cable, conductor and domain numbers and loop variables
-  integer :: cable
+  integer :: cable,cable_t
   integer :: tot_n_cables  
   integer :: tot_n_cables_without_ground_plane
  
@@ -115,7 +117,9 @@ USE maths
   integer :: first_external_conductor
  
 ! variables for looking at overshields and determining which domain cables are in
-  integer :: tot_n_overshields,n_overshield_domains,overshield
+  integer :: tot_n_overshields,n_overshield_domains,overshield,reference_overshield
+  real(dp)::radius_reference_overshield
+  
   integer,allocatable  :: overshield_shape(:)
   real(dp),allocatable :: overshield_x(:)
   real(dp),allocatable :: overshield_y(:)
@@ -125,11 +129,12 @@ USE maths
   real(dp),allocatable :: overshield_h(:)
   integer,allocatable  :: overshield_domain(:)
   integer,allocatable  :: overshield_reference_terminal_conductor(:)
+  integer,allocatable  :: overshield_to_cable_number(:)
  
   integer,allocatable  :: overshield_n_conductors(:)
   integer,allocatable  :: overshield_terminal_conductor(:,:) 
  
-  real(dp) :: cable_x,cable_y,dist_cable_to_overshield_centre
+  real(dp) :: cable_x,cable_y,dist_cable_to_overshield_centre,cable_r
   logical :: is_overshield_domain
  
 ! referencing arrays
@@ -160,7 +165,7 @@ USE maths
   integer :: dim
   integer :: row,col
   integer :: row_l,col_l
-  integer :: row_g,col_g
+  integer :: row_g,col_g,row_g_ref
  
   character(LEN=3) :: conductor_string  ! used to generate the conductor labels in the bundle
  
@@ -291,12 +296,14 @@ USE maths
     ALLOCATE( overshield_w2(1:tot_n_overshields)  )
     ALLOCATE( overshield_h(1:tot_n_overshields)  )
     ALLOCATE( overshield_reference_terminal_conductor(1:tot_n_overshields) )
+    ALLOCATE( overshield_to_cable_number(1:tot_n_overshields) )
  
     overshield=0
     conductor_count=0
     do cable=1,tot_n_cables
       if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_overshield) then 
         overshield=overshield+1
+        overshield_to_cable_number(overshield)=cable
         overshield_shape(overshield)=circle
         overshield_x(overshield)=bundle%cable_x_offset(cable)
         overshield_y(overshield)=bundle%cable_y_offset(cable)
@@ -305,6 +312,20 @@ USE maths
         overshield_w2(overshield)=0d0
         overshield_h(overshield)=0d0
         overshield_reference_terminal_conductor(overshield)=conductor_count+1
+        
+      else
+
+! This is not an overshield
+
+        if ( (bundle%cable(cable)%cable_type.NE.cable_geometry_type_ground_plane) .AND. (overshield.NE.0) ) then
+
+! check that oveshields are the last conductors defined in the bundle_spec file
+          run_status='ERROR, Overshields should be the last cables defined in the bundle specification'
+          CALL write_program_status()
+          STOP 1
+        
+        end if      
+      
       end if
       conductor_count=conductor_count+bundle%cable(cable)%tot_n_conductors
     end do
@@ -330,47 +351,175 @@ USE maths
   cable_reference_conductor(1:bundle%n_cables)=0
   n_overshield_domains=0     ! counter for overshield domains
  
+  write(*,*)'**********************************************************************'
+  
   do cable=1,tot_n_cables
  
+    write(*,*)
+    write(*,*)'Checking cable number',cable
+  
     if (bundle%cable(cable)%cable_type.NE.cable_geometry_type_overshield) then 
-! if this cable is not an overshield then see whether it is WITHIN an overshield
+    
+! IF THIS CABLE IS NOT AN OVERSHIELD THEN SEE WHETHER IT IS WITHIN AN OVERSHIELD
+  
+      write(*,*)'THIS CABLE IS NOT AN OVERSHIELD'
  
       cable_x=bundle%cable_x_offset(cable)
       cable_y=bundle%cable_y_offset(cable)
+
+! find the smallest overshield enclosing this cable      
+
+      reference_overshield=0
+      radius_reference_overshield=1d30
  
       do overshield=1,tot_n_overshields
-      
+            
         dist_cable_to_overshield_centre=sqrt( (cable_x-overshield_x(overshield))**2+  &
-                                              (cable_y-overshield_y(overshield))**2  )
+                                                (cable_y-overshield_y(overshield))**2  )
+                                              
         if ( dist_cable_to_overshield_centre.LT.overshield_r(overshield) ) then
-! we conclude that this cable is with the overshield
+! we conclude that this cable is within the overshield
+
+          write(*,*)'This cable is inside overshield number',overshield
  
-          if (overshield_domain(overshield).EQ.0) then 
+          if (radius_reference_overshield.GT.overshield_r(overshield)) then
+! this overshield is smaller than the current reference overshield so this becomes the new reference overshield
+            
+            reference_overshield=overshield
+            radius_reference_overshield=overshield_r(overshield)
+            
+            write(*,*)'This cable reference overshield has been updated to',reference_overshield
+
+          end if  ! is this the smallest overshield?
+          
+        end if  ! is cable within this overshield?
+                 
+      end do
+      
+      if (reference_overshield.NE.0) then
+! this cable has been identified as being within an overshield
+
+        if (overshield_domain(reference_overshield).EQ.0) then 
  
 ! this overshield domain is unallocated so create a new domain
  
-            tot_n_domains=tot_n_domains+1
-            n_overshield_domains=n_overshield_domains+1
-            if(verbose) write(*,*)'Creating overshield domain=',tot_n_domains
-            overshield_domain(overshield)=tot_n_domains
+          tot_n_domains=tot_n_domains+1
+          n_overshield_domains=n_overshield_domains+1
+          if(verbose) write(*,*)'Creating overshield domain=',tot_n_domains
+          overshield_domain(reference_overshield)=tot_n_domains
+            
+          global_domain_reference_conductor(tot_n_domains)=overshield_reference_terminal_conductor(reference_overshield)
  
-            global_domain_reference_conductor(tot_n_domains)=overshield_reference_terminal_conductor(overshield)
+          write(*,*)'This overshield domain is unallocated so create a new domain, domain=',tot_n_domains
  
-          end if ! this overshield domain is unallocated
+        end if ! this overshield domain is unallocated
+          
+        cable_reference_domain(cable)=overshield_domain(reference_overshield)
  
-          cable_reference_domain(cable)=overshield_domain(overshield)
+        write(*,*)'Setting the cable reference domain to',overshield_domain(reference_overshield)
  
-        end if   ! this cable is with the overshield
+      end if   ! this cable is within the overshield
+             
+      terminal_conductor=terminal_conductor+bundle%cable(cable)%tot_n_conductors ! update the terminal conductor count
  
-      end do   ! next overshield to check
+    else  ! THIS CABLE IS AN OVERSHIELD SO WE NEED TO CHECK WHETHER IT IS INSIDE ANY OTHER OVERSHIELDS   
+   
+      write(*,*)'THIS CABLE IS AN OVERSHIELD'
+       
+      cable_x=bundle%cable_x_offset(cable)
+      cable_y=bundle%cable_y_offset(cable)
+      cable_r=bundle%cable(cable)%parameters(1)  ! overshield radius
  
-      terminal_conductor=terminal_conductor+bundle%cable(cable)%tot_n_conductors ! update the terminal conductor count
+! find the smallest overshield enclosing this overshield if there is one   
+
+      reference_overshield=0
+      radius_reference_overshield=1d30
+      
+      do overshield=1,tot_n_overshields
+      
+        cable_t=overshield_to_cable_number(overshield)
+        if (cable_t.NE.cable) then   ! don't check an overshield against itself
+      
+          dist_cable_to_overshield_centre=sqrt( (cable_x-overshield_x(overshield))**2+  &
+                                                (cable_y-overshield_y(overshield))**2  )
+                                              
+          if ( dist_cable_to_overshield_centre.LT.overshield_r(overshield) ) then
+! we conclude that this cable is within the overshield
+
+            write(*,*)'This overshield centre is inside overshield number',overshield
+
+!! Check that can be removed when we are happy with nexted overshields code            
+!            run_status='ERROR, Nested overshields cannot be used in SACAMOS at the moment'
+!            CALL write_program_status()
+!            STOP 1
+
+            if (radius_reference_overshield.GT.overshield_r(overshield)) then
+! this overshield is smaller than the current reference overshield so this becomes the new reference overshield
+            
+              if (overshield_r(overshield).GT.cable_r) then
+! the cable is smaller than the overshield and hence is inside it
+
+! Check that the lowest number overshield is the smaller. Overshields should be numbered inside to out
+
+                if (cable_t.LT.cable) then
+              
+                  run_status='ERROR, Nested overshields should be numbered in the order inside to out'
+                  CALL write_program_status()
+                  STOP 1
+              
+                end if
+
+                write(*,*)'This overshield is smaller than the reference overshield and is therefore inside it'
+           
+                reference_overshield=overshield
+                radius_reference_overshield=overshield_r(overshield)
+            
+                write(*,*)'This cable reference overshield has been updated to',reference_overshield
+              
+              end if  ! This overshield is smaller than the reference overshield and is therefore inside it
+              
+            end if  ! is this the smallest overshield?
+          
+          end if  ! is cable centre within this overshield?
+          
+        end if  ! don't check an overshield against itself
+       
+      end do
  
+      if (reference_overshield.NE.0) then
+! this cable has been identified as being within an overshield
+
+        if (overshield_domain(reference_overshield).EQ.0) then 
+          
+! this overshield domain is unallocated so create a new domain
+  
+          tot_n_domains=tot_n_domains+1
+          n_overshield_domains=n_overshield_domains+1
+          if(verbose) write(*,*)'Creating overshield domain=',tot_n_domains
+          overshield_domain(reference_overshield)=tot_n_domains
+            
+          global_domain_reference_conductor(tot_n_domains)=overshield_reference_terminal_conductor(reference_overshield)
+            
+          write(*,*)'This overshield domain is unallocated so create a new domain, domain=',tot_n_domains
+                        
+        end if ! this overshield domain is unallocated
+          
+        cable_reference_domain(cable)=overshield_domain(reference_overshield)
+          
+        write(*,*)'Setting the cable reference domain to',overshield_domain(reference_overshield)
+          
+      end if   ! this cable is within the overshield
+             
+      terminal_conductor=terminal_conductor+bundle%cable(cable)%tot_n_conductors ! update the terminal conductor count
+     
+       
     end if ! not an overshield
  
   end do ! next cable
  
   if (verbose) write(*,*) 'n_overshield_domains=',n_overshield_domains
+  
+  write(*,*)'**********************************************************************'
  
 ! 3c part 2. loop over the overshield domains and count the number of external conductors in each one  
  
@@ -736,6 +885,7 @@ USE maths
     end do ! next local domain
  
   end do ! next cable
+  
  
 ! 6.  Now that we have counted the number of domains we can allocate the 
 !     domain based number of conductors and L and C matrices plus
@@ -854,6 +1004,12 @@ USE maths
   terminal_conductor=0       ! this is the external conductor number
   domain_count=0
  
+  if (verbose) then
+    write(*,*)''
+    write(*,*)'SET MI, MV matrices'
+    write(*,*)''
+  end if
+  
   do cable=1,tot_n_cables_without_ground_plane
  
     local_n_domains=bundle%cable(cable)%tot_n_domains
@@ -864,6 +1020,7 @@ USE maths
 ! copy the cable based MI and MV matrices to the global structure. The cable based
 ! MI and MV matrices are on the basis of terminal_conductor numbering
     if (verbose) then
+      write(*,*)''
       write(*,*)'Copy MI, MV, cable=',cable,' of', tot_n_cables_without_ground_plane
       write(*,*)'Matrix dimension',n_cable_conductors+1
       write(*,*)'intial terminal_conductor count=',terminal_conductor
@@ -880,7 +1037,7 @@ USE maths
         else
           row_g=terminal_conductor_to_reference_global_domain_conductor(terminal_conductor+row_l-1) ! last conductor for this cable and the reference    
         end if
-      
+     
         if (col_l.NE.n_cable_conductors+1) then
           col_g=terminal_conductor+col_l
         else
@@ -890,12 +1047,44 @@ USE maths
         bundle%global_MI%mat(row_g,col_g)=bundle%cable(cable)%MI%mat(row_l,col_l)
         bundle%global_MV%mat(row_g,col_g)=bundle%cable(cable)%MV%mat(row_l,col_l)
  
-        if (verbose) write(*,*)'Copy MV element',row_g,col_g,bundle%global_MV%mat(row_g,col_g)
+!      if (verbose) write(*,*)'Copy MV element',row_g,col_g,bundle%global_MV%mat(row_g,col_g)
  
+          if (verbose) then
+          write(*,*)'Copy MI element',' row_g=',row_g,' col_g=',col_g,      &
+                                    ' row_l=',row_l,' col_l=',col_l,bundle%cable(cable)%MI%mat(row_l,col_l)
+        end if
+        
       end do ! next col
  
     end do ! next row
  
+    if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_overshield) then 
+! This is an overshield so invoke the special process for overshields
+
+! Get the global domain number of the shield conductor
+      row_g=terminal_conductor_to_global_domain_conductor(terminal_conductor+1)
+      
+! Get the global domain number of the shield reference conductor
+      row_g_ref=terminal_conductor_to_reference_global_domain_conductor(terminal_conductor+1)
+      
+! Add all the non-zero elements of the MI matrix shield row to the shield reference row
+
+      write(*,*)'Overshield cable',cable,' Copying row',row_g,' to row',row_g_ref
+                      
+      do col_g=1,bundle%global_MI%dim
+
+! In the global MI, MV structure rows correspond to the global domain conductor number, cols correspond to the terminal conductor number 
+                
+        if (bundle%global_MI%mat(row_g,col_g).NE.0d0) then
+        
+          bundle%global_MI%mat(row_g_ref,col_g)=bundle%global_MI%mat(row_g,col_g)
+        
+        end if
+                
+      end do ! next col
+      
+    end if ! overshield
+    
     do local_domain=1,local_n_domains  ! exclude the external domain for now
  
       if (local_domain.NE.local_n_domains) then
@@ -1363,6 +1552,7 @@ write(*,*)&#39;****************************************************************&#39;
   if (allocated( overshield_h )) DEALLOCATE( overshield_h)
   if (allocated( overshield_domain )) DEALLOCATE( overshield_domain)
   if (allocated( overshield_reference_terminal_conductor)) DEALLOCATE( overshield_reference_terminal_conductor )
+  if (allocated( overshield_to_cable_number)) DEALLOCATE( overshield_to_cable_number )
   if (allocated( global_domain_reference_conductor )) DEALLOCATE( global_domain_reference_conductor )
   DEALLOCATE( cable_reference_conductor  )
   DEALLOCATE( cable_reference_domain  )
@@ -3,9 +3,9 @@
 #MOD_bundle_lib_dir
 ./
 2    # Number of cables in bundle, cable list follows: name then x,y of cable
-overshield
-0.0     0.01   ! x and y coordinates of the cable centre
 inner_wire
 0.0     0.01   ! x and y coordinates of the cable centre
+overshield
+0.0     0.01   ! x and y coordinates of the cable centre
 ground_plane
  
@@ -15,23 +15,23 @@ coax_over_ground_plane
 90.0 0.0	ktheta kphi (degrees)
 -1.0  0.0	Etheta Ephi
 # End 1 termination model
-1.0      End 1 voltage source list
+0.0      End 1 voltage source list
+1.0
 0.0
-0.0
-10.0  	End 1 impedance list
-100.0
+100.0  	End 1 impedance list
+10.0
 0.0
 # End 2 termination model
 0.0        End 2 voltage source list
 0.0
 0.0
-20.0  	   End 2 impedance list
-25.0
+25.0  	   End 2 impedance list
+20.0
 0.0
 # Type of analysis
 AC
 log          # frequency scale (log or lin)
 1e3 1e8 1000 # fmin fmax number_of_frequencies
 # Output conductor number and end number
-1    1
+2    1
 lin   # output type (lin or dB)