From 44c11f06bad21a8b0e1bf76d968102ea52580a6e Mon Sep 17 00:00:00 2001 From: chris.smartt Date: Tue, 19 Jun 2018 15:17:27 +0100 Subject: [PATCH] Include software updates for infinite ground plane model, new flex cable and iterative solver in Laplace' --- SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 | 7 +++++-- SRC/BUNDLE_DOMAIN_CREATION/plot_bundle_cross_section.F90 | 9 ++++++++- SRC/BUNDLE_DOMAIN_CREATION/set_conductor_positions_for_Einc.F90 | 7 ++++--- SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 | 6 ++++-- SRC/CABLE_MODULES/ML_flex_cable.F90 | 471 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/CABLE_MODULES/Makefile | 1 + SRC/CABLE_MODULES/cable_checks.F90 | 136 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/CABLE_MODULES/cable_module.F90 | 38 ++++++++++++++++++++++++++++++++++---- SRC/CABLE_MODULES/flex_cable.F90 | 2 +- SRC/GENERAL_MODULES/constants.F90 | 2 ++ SRC/GENERAL_MODULES/general_module.F90 | 9 +++++++++ SRC/GENERAL_MODULES/generate_shapes.F90 | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/GENERAL_MODULES/plot_geometry.F90 | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/MATHS_MODULES/cmatrix.F90 | 3 +++ SRC/MATHS_MODULES/dmatrix.F90 | 2 ++ SRC/PUL_PARAMETER_CALCULATION/Aprod.F90 | 207 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/PUL_PARAMETER_CALCULATION/CG_solve.F90 | 382 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SRC/PUL_PARAMETER_CALCULATION/Laplace.F90 | 179 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------- SRC/PUL_PARAMETER_CALCULATION/Makefile | 2 +- SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 | 442 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------- SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 | 95 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------- SRC/cable_bundle_model_builder.F90 | 15 +++++++++++++++ SRC/cable_model_builder.F90 | 35 ++++++++++++++++++++++++++++++----- TEST_CASES/FLEX_CABLE_3_AND_SINGLE_WIRE_EINC/flex_cable_3.cable_spec | 18 +++++++++++------- TEST_CASES/FLEX_CABLE_3_CONDUCTOR/flex_cable_3.cable_spec | 35 +++++++++++++++++++++++++++-------- TEST_CASES/FLEX_CABLE_3_EINC/flex_cable_3.cable_spec | 27 +++++++++++++++++++++++++++ TEST_CASES/FLEX_CABLE_3_GROUND_PLANE/flex_cable_3.cable_spec | 26 ++++++++++++++++++++++++++ TEST_CASES/FLEX_CABLE_3_OVERSHIELD/flex_cable_3.cable_spec | 26 ++++++++++++++++++++++++++ TEST_CASES/FLEX_CABLE_DIELECTRIC_3_CONDUCTOR/flex_cable_3.cable_spec | 52 ++++++++++++++++++++++++++++++++++++++++++++++++++++ TEST_CASES/FLEX_CABLE_LOSSY_2_CONDUCTOR/flex_cable_2.cable_spec | 26 ++++++++++++++++++++++++++ TEST_CASES/generate_spice_cable_bundle_model | 1 - clean_project | 16 ++++++++++++++++ 32 files changed, 2229 insertions(+), 187 deletions(-) create mode 100644 SRC/CABLE_MODULES/ML_flex_cable.F90 create mode 100644 SRC/PUL_PARAMETER_CALCULATION/Aprod.F90 create mode 100644 SRC/PUL_PARAMETER_CALCULATION/CG_solve.F90 create mode 100755 clean_project diff --git a/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 b/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 index c5438bd..be8b032 100644 --- a/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 +++ b/SRC/BUNDLE_DOMAIN_CREATION/create_global_domain_structure.F90 @@ -72,6 +72,7 @@ ! 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 ! SUBROUTINE create_global_domain_structure(bundle) @@ -994,7 +995,8 @@ USE maths PUL%y(conductor)=bundle%cable_y_offset(cable) PUL%rtheta(conductor)=bundle%cable_angle(cable) - PUL%o(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_ox + PUL%ox(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_ox + PUL%oy(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_oy PUL%r(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_radius PUL%rw(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_width PUL%rw2(conductor)=bundle%cable(cable)%external_model(local_conductor)%conductor_width2 @@ -1002,7 +1004,8 @@ USE maths PUL%rd(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_radius PUL%rdw(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_width PUL%rdh(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_height - PUL%rdo(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_ox + PUL%rdox(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_ox + PUL%rdoy(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_oy PUL%epsr(conductor)=bundle%cable(cable)%external_model(local_conductor)%dielectric_epsr end if ! not a ground plane diff --git a/SRC/BUNDLE_DOMAIN_CREATION/plot_bundle_cross_section.F90 b/SRC/BUNDLE_DOMAIN_CREATION/plot_bundle_cross_section.F90 index 720f8b4..c581422 100644 --- a/SRC/BUNDLE_DOMAIN_CREATION/plot_bundle_cross_section.F90 +++ b/SRC/BUNDLE_DOMAIN_CREATION/plot_bundle_cross_section.F90 @@ -45,7 +45,7 @@ ! 14/6/2016 CJS Create the arrays of x and y coordinates of each conductor which are used in the incident field excitation ! 5/9/2016 CJS include additional shielded cable types with surface impedance loss models ! 20/4/2017 CJS Separate the set_conductor_positions_for_Einc process from plot_bundle_cross_section and make use of generate_shapes.F90 -! +! 16/3/2018 CJS plot ML_flex_cable ! SUBROUTINE plot_bundle_cross_section(bundle) @@ -151,6 +151,13 @@ USE maths bundle%cable_x_offset(cable), & bundle%cable_y_offset(cable), & bundle%cable_angle(cable),xmin,xmax,ymin,ymax ) + + else if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_ML_flex_cable) then + + CALL ML_flex_cable_plot( bundle%cable(cable), & + bundle%cable_x_offset(cable), & + bundle%cable_y_offset(cable), & + bundle%cable_angle(cable),xmin,xmax,ymin,ymax ) else if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_dconnector) then diff --git a/SRC/BUNDLE_DOMAIN_CREATION/set_conductor_positions_for_Einc.F90 b/SRC/BUNDLE_DOMAIN_CREATION/set_conductor_positions_for_Einc.F90 index bec215a..c9beee1 100644 --- a/SRC/BUNDLE_DOMAIN_CREATION/set_conductor_positions_for_Einc.F90 +++ b/SRC/BUNDLE_DOMAIN_CREATION/set_conductor_positions_for_Einc.F90 @@ -46,7 +46,7 @@ ! 14/6/2016 CJS Create the arrays of x and y coordinates of each conductor which are used in the incident field excitation ! 5/9/2016 CJS include additional shielded cable types with surface impedance loss models ! 20/4/2017 CJS Separate the set_conductor_positions_for_Einc process from plot_bundle_cross_section -! +! 16/3/2018 CJS add y offset for ML_flex_cable ! SUBROUTINE set_conductor_positions_for_Einc(bundle) @@ -98,7 +98,8 @@ USE maths do cable=1,bundle%n_cables ! loop over the conductors in the cable - if (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_flex_cable) then + if ((bundle%cable(cable)%cable_type.EQ.cable_geometry_type_flex_cable).OR. & + (bundle%cable(cable)%cable_type.EQ.cable_geometry_type_ML_flex_cable) ) then ! We include a fix here for flex cables which have to have their conductor offsets calculated here ! work out the centre of this conductor before rotation @@ -107,7 +108,7 @@ USE maths do i=1,bundle%cable(cable)%tot_n_conductors FC_x0=bundle%cable(cable)%external_model(i)%conductor_ox - FC_y0=0.0 + FC_y0=bundle%cable(cable)%external_model(i)%conductor_oy ! work out the centre of this conductor when the flex cable is rotated and offset conductor=conductor+1 bundle%conductor_x_offset(conductor)=FC_x0*cos(theta)-FC_y0*sin(theta)+bundle%cable_x_offset(cable) diff --git a/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 b/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 index f36591e..059e0d1 100644 --- a/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 +++ b/SRC/CABLE_BUNDLE_MODULES/cable_bundle_module.F90 @@ -892,7 +892,8 @@ CONTAINS theta=bundle%cable_angle(cable1) type1=bundle%cable(cable1)%cable_type - if (bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable) then + if ((bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable).AND. & + (bundle%cable(cable1)%cable_type.NE.cable_geometry_type_ML_flex_cable) ) then nec1=bundle%cable(cable1)%n_external_conductors else nec1=1 @@ -953,7 +954,8 @@ CONTAINS theta=bundle%cable_angle(cable2) type2=bundle%cable(cable2)%cable_type - if (bundle%cable(cable2)%cable_type.NE.cable_geometry_type_flex_cable) then + if ((bundle%cable(cable1)%cable_type.NE.cable_geometry_type_flex_cable).AND. & + (bundle%cable(cable1)%cable_type.NE.cable_geometry_type_ML_flex_cable) ) then nec2=bundle%cable(cable2)%n_external_conductors else nec2=1 diff --git a/SRC/CABLE_MODULES/ML_flex_cable.F90 b/SRC/CABLE_MODULES/ML_flex_cable.F90 new file mode 100644 index 0000000..519ae9f --- /dev/null +++ b/SRC/CABLE_MODULES/ML_flex_cable.F90 @@ -0,0 +1,471 @@ +! +! This file is part of SACAMOS, State of the Art CAble MOdels in Spice. +! It was developed by the University of Nottingham and the Netherlands Aerospace +! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK. +! +! Copyright (C) 2016-2017 University of Nottingham +! +! SACAMOS is free software: you can redistribute it and/or modify it under the +! terms of the GNU General Public License as published by the Free Software +! Foundation, either version 3 of the License, or (at your option) any later +! version. +! +! SACAMOS is distributed in the hope that it will be useful, but +! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! A copy of the GNU General Public License version 3 can be found in the +! file GNU_GPL_v3 in the root or at . +! +! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to +! the GNU Lesser General Public License. A copy of the GNU Lesser General Public +! License version can be found in the file GNU_LGPL in the root of EISPACK +! (/SRC/EISPACK ) or at . +! +! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk +! +! File Contents: +! SUBROUTINE ML_flex_cable_set_parameters +! SUBROUTINE ML_flex_cable_set_internal_domain_information +! SUBROUTINE ML_flex_cable_plot +! +! NAME +! ML_flex_cable_set_parameters +! +! AUTHORS +! Chris Smartt +! +! DESCRIPTION +! Set the overall parameters for a ML_flex_cable cable +! +! COMMENTS +! +! +! HISTORY +! +! started 26/9/2016 CJS +! 16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions +! +! +SUBROUTINE ML_flex_cable_set_parameters(cable) + +USE type_specifications + +IMPLICIT NONE + +! variables passed to subroutine + + type(cable_specification_type),intent(INOUT) :: cable + +! local variables + +! START + + cable%cable_type=cable_geometry_type_ML_flex_cable + cable%tot_n_conductors=0 ! this is set in the cable specification file + cable%tot_n_domains=1 + cable%n_external_conductors=0 ! this is set in the cable specification file + cable%n_internal_conductors=0 + cable%n_internal_domains=0 + cable%n_parameters=0 + cable%n_dielectric_filters=1 + cable%n_transfer_impedance_models=0 + +END SUBROUTINE ML_flex_cable_set_parameters +! +! NAME +! ML_flex_cable_set_internal_domain_information +! +! AUTHORS +! Chris Smartt +! +! DESCRIPTION +! Set the overall parameters for a ML_flex_cable cable +! +! COMMENTS +! Set the dimension of the domain transformation matrices to include an external reference conductor for the cable +! +! +! HISTORY +! +! started 13/4/2016 CJS +! 8/5/2017 CJS: Include references to Theory_Manual +! +! +SUBROUTINE ML_flex_cable_set_internal_domain_information(cable) + +USE type_specifications +USE general_module +USE constants + +IMPLICIT NONE + +! variables passed to subroutine + + type(cable_specification_type),intent(INOUT) :: cable + +! local variables + + integer :: nc + integer :: dim + integer :: i,prow,prow2,row2,row,col,conductivity_row + + real(dp) :: full_width,ML_flex_cable_xmin + +! variables for cable parameter checks + logical :: cable_spec_error + real(dp) :: wd ! dielectric width + real(dp) :: hd ! dielectric height + integer :: nrows ! number of rows of conductors + + real(dp) :: w ! conductor width + real(dp) :: h ! conductor height + real(dp) :: s ! conductor separation + real(dp) :: ox ! conductor offset x + real(dp) :: oy ! conductor offset y + integer :: ncrow ! number of conductors in this row + integer :: conductor ! counter for conductors + + real(dp) :: w1,h1,w2,h2,ox1,ox2,oy1,oy2 + + real(dp) :: sigma + + type(Sfilter) :: epsr + + integer :: nclocal + + character(LEN=error_message_length) :: message + + character(LEN=2) :: conductor_string + + integer :: check_type + +! START + + nc=cable%tot_n_conductors + +! set the information related ot the number of conductors + + cable%n_external_conductors=nc + +! Set the domain decomposition matrices ! Theory_Manual_Eqn 6.17, 6.18 + +! The dimension of the domain transformation matrices is the number of conductors+1 + dim=nc+1 + + cable%MI%dim=dim + ALLOCATE(cable%MI%mat(dim,dim)) + cable%MI%mat(:,:)=0d0 + + cable%MV%dim=dim + ALLOCATE(cable%MV%mat(dim,dim)) + cable%MV%mat(:,:)=0d0 + + do i=1,nc + row=i + col=row + cable%MI%mat(row,col)=1d0 + end do + do i=1,dim + row=dim + col=i + cable%MI%mat(row,col)=1d0 + end do + + do i=1,nc + row=i + col=row + cable%MV%mat(row,col)=1d0 + cable%MV%mat(row,dim)=-1d0 + end do + cable%MV%mat(dim,dim)=1d0 + +! Set the local reference conductor numbering + ALLOCATE( cable%local_reference_conductor(nc) ) + cable%local_reference_conductor(1:nc)=0 ! external domain conductor, reference not known + +! Set the local domain information: include a reference conductor in the count + ALLOCATE( cable%local_domain_n_conductors(1:cable%tot_n_domains) ) + cable%local_domain_n_conductors(1)=nc+1 ! external domain + +! Read the parameter list and set the conductor information + + nrows=NINT(cable%parameters(3)) + +! check that the number of conductors is consistent + prow=3 + nclocal=0 + + do row=1,nrows ! loop over rows of conductors + + ncrow=NINT(cable%parameters(prow+6)) + nclocal=nclocal+ncrow + prow=prow+6 + + end do ! next row of conductors + +! check that the number of conductors is consistent + if (nclocal.NE.nc) then + write(message,*)' Nc in .cable file=',nc,' conductor count=',nclocal + run_status='ERROR in cable_model_builder, inconsistent conductor count for flex_cable:' & + //trim(cable%cable_name)//'. '//trim(message) + CALL write_program_status() + STOP 1 + end if + +! conductivity for the conductor loss models + conductivity_row=prow+1 + sigma=cable%parameters(conductivity_row) ! conductivity + + ALLOCATE( cable%external_model(cable%n_external_conductors) ) + + prow=3 + conductor=0 ! reset the conductor count + + do row=1,nrows ! loop over rows of conductors + +! get the parameters for this row of conductors + ox=cable%parameters(prow+1) ! conductor offset in x + oy=cable%parameters(prow+2) ! conductor offset in y + w=cable%parameters(prow+3) ! conductor width + h=cable%parameters(prow+4) ! conductor height + s=cable%parameters(prow+5) ! conductor separation + ncrow=NINT(cable%parameters(prow+6)) + + full_width=w*ncrow+s*(ncrow-1) + ML_flex_cable_xmin=-full_width/2d0 + + do i=1,ncrow ! loop over the conductors in this row + + conductor=conductor+1 + +! set the conductor impedance model for this conductor + cable%conductor_impedance(conductor)%impedance_model_type=impedance_model_type_rectangular_with_conductivity + cable%conductor_impedance(conductor)%width=w + cable%conductor_impedance(conductor)%height=h + cable%conductor_impedance(conductor)%conductivity=sigma + +! Set the external domain conductor information + + CALL reset_external_conductor_model(cable%external_model(conductor)) + cable%external_model(conductor)%conductor_type=rectangle + cable%external_model(conductor)%conductor_width=w + cable%external_model(conductor)%conductor_width2=w + cable%external_model(conductor)%conductor_height=h + +! work out the offset of the ith conductor + cable%external_model(conductor)%conductor_ox=ox+ML_flex_cable_xmin+(w+s)*(i-1)+w/2d0 + cable%external_model(conductor)%conductor_oy=oy + + end do ! next conductor in this row + + prow=prow+6 + + end do ! next row of conductors + +! dielectric data + + epsr=cable%dielectric_filter(1) + wd=cable%parameters(1) ! dielectric width + hd=cable%parameters(2) ! dielectric height + +! do some consistency checks + cable_spec_error=.FALSE. ! assume no errors initially + message='' + +! Do some intersection checks on the flex cable + + prow=3 + do row=1,nrows ! loop over rows of conductors + +! get the parameters for this row of conductors + ox1=cable%parameters(prow+1) ! conductor offset in x + oy1=cable%parameters(prow+2) ! conductor offset in y + w=cable%parameters(prow+3) ! conductor width + h1=cable%parameters(prow+4) ! conductor height + s=cable%parameters(prow+5) ! conductor separation + ncrow=NINT(cable%parameters(prow+6)) + + w1=w*ncrow+s*(ncrow-1) ! w1 is the full width of the row of conductors + +! check whether this row of conductors is well specified - conductor width, height and spearation >0 + CALL flex_cable_check(ncrow,w,h1,s,0d0,0d0,cable_spec_error,cable%cable_name,message) + + if (cable_spec_error) then + run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message) + CALL write_program_status() + STOP 1 + end if + +! check whether this row of conductors intersects the dielectric boundary of the flex cable + check_type=2 + CALL ML_flex_cable_check(w1,h1,ox1,oy1,wd,hd,0d0,0d0,check_type,cable_spec_error,cable%cable_name,message) + + if (cable_spec_error) then + run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message) + CALL write_program_status() + STOP 1 + end if + +! loop over all the other rows to check for intersections between the rows of conductors + + prow=prow+6 + + prow2=prow + + do row2=row+1,nrows + +! get the parameters for this row of conductors + ox2=cable%parameters(prow2+1) ! conductor offset in x + oy2=cable%parameters(prow2+2) ! conductor offset in y + w=cable%parameters(prow2+3) ! conductor width + h2=cable%parameters(prow2+4) ! conductor height + s=cable%parameters(prow2+5) ! conductor separation + ncrow=NINT(cable%parameters(prow2+6)) + + w2=w*ncrow+s*(ncrow-1) ! w2 is the full width of the second row of conductors + + +! check whether the two rows of conductors intersect + check_type=1 + CALL ML_flex_cable_check(w1,h1,ox1,oy1,w2,h2,ox2,oy2,check_type,cable_spec_error,cable%cable_name,message) + + if (cable_spec_error) then + run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message) + CALL write_program_status() + STOP 1 + end if + + prow2=prow2+6 + + end do ! next row2 for checking intersection between rows of conductors + + end do ! next row of conductors to check + + CALL dielectric_check(epsr,cable_spec_error,cable%cable_name,message) + + if (cable_spec_error) then + run_status='ERROR in cable_model_builder, error on parameters for cable:'//trim(cable%cable_name)//'. '//trim(message) + CALL write_program_status() + STOP 1 + end if + +! add a dielectric region to the first conductor which encloses the whole cable + +! write the dielectric which is offset from the conductors +! The dielectric is centred at the cable centre + + cable%external_model(1)%dielectric_width=wd + cable%external_model(1)%dielectric_height=hd + cable%external_model(1)%dielectric_ox=0d0 + cable%external_model(1)%dielectric_oy=0d0 + cable%external_model(1)%dielectric_epsr=epsr + + CALL deallocate_Sfilter(epsr) + + ALLOCATE( cable%conductor_label(1:cable%tot_n_conductors) ) + do i=1,cable%tot_n_conductors + write(conductor_string,'(I2)')i + cable%conductor_label(i)='Cable name: '//trim(cable%cable_name)// & + '. type: '//trim(cable%cable_type_string)//'. conductor '//conductor_string//' : ML_flex_cable conductor' + end do + +END SUBROUTINE ML_flex_cable_set_internal_domain_information +! +! NAME +! ML_flex_cable_plot +! +! AUTHORS +! Chris Smartt +! +! DESCRIPTION +! plot ML_flex_cable cable +! +! COMMENTS +! +! +! HISTORY +! +! started 23/9/2016 CJS +! +! +SUBROUTINE ML_flex_cable_plot(cable,x_offset,y_offset,theta,xmin,xmax,ymin,ymax) + +USE type_specifications +USE general_module + +IMPLICIT NONE + +! variables passed to subroutine + + type(cable_specification_type),intent(IN) :: cable + + real(dp),intent(IN) :: x_offset,y_offset,theta + real(dp),intent(INOUT) :: xmin,xmax,ymin,ymax + +! local variables + + integer nc + + real(dp) :: full_width,ML_flex_cable_xmin + real(dp) :: xoff,yoff,x,y,w,h,s,wd,hd,ox,oy + + integer nrows,ncrow,row,prow + integer i + +! START + +! plot ML_flex_cable conductor + + nc=cable%tot_n_conductors + + nrows=NINT(cable%parameters(3)) + + prow=3 + + do row=1,nrows ! loop over rows of conductors + +! get the parameters for this row of conductors + ox=cable%parameters(prow+1) ! conductor offset in x + oy=cable%parameters(prow+2) ! conductor offset in y + w=cable%parameters(prow+3) ! conductor width + h=cable%parameters(prow+4) ! conductor height + s=cable%parameters(prow+5) ! conductor separation + ncrow=NINT(cable%parameters(prow+6)) + + full_width=w*ncrow+s*(ncrow-1) + ML_flex_cable_xmin=-full_width/2d0 + + do i=1,ncrow ! loop over the conductors in this row + +! work out the centre of this conductor before rotation + xoff=ox+ML_flex_cable_xmin+(w+s)*(i-1)+w/2d0 + yoff=oy + +! work out the centre of this conductor when the flex cable is rotated and offset + x=x_offset+xoff*cos(theta)-yoff*sin(theta) + y=y_offset+xoff*sin(theta)+yoff*cos(theta) + +! write the conductor + CALL write_rectangle(x,y,w,h,theta,conductor_geometry_file_unit,xmin,xmax,ymin,ymax) + + end do ! next conductor in this row + + prow=prow+6 + + end do ! next row of conductors + +! write the dielectric which is offset from the conductors + + wd=cable%external_model(1)%dielectric_width + hd=cable%external_model(1)%dielectric_height + CALL write_rectangle(x_offset,y_offset,wd,hd,theta,dielectric_geometry_file_unit,xmin,xmax,ymin,ymax) + + RETURN + +END SUBROUTINE ML_flex_cable_plot + + + + diff --git a/SRC/CABLE_MODULES/Makefile b/SRC/CABLE_MODULES/Makefile index ca85276..f0b00a9 100644 --- a/SRC/CABLE_MODULES/Makefile +++ b/SRC/CABLE_MODULES/Makefile @@ -8,6 +8,7 @@ $(OBJ_MOD_DIR)/cable_module.o: cable_module.F90 $(TYPE_SPEC_MODULE) $(GENERAL_M twinax.F90 \ spacewire.F90 \ flex_cable.F90 \ + ML_flex_cable.F90 \ Dconnector.F90 \ overshield.F90 \ ground_plane.F90 \ diff --git a/SRC/CABLE_MODULES/cable_checks.F90 b/SRC/CABLE_MODULES/cable_checks.F90 index 2885e16..8efc41e 100644 --- a/SRC/CABLE_MODULES/cable_checks.F90 +++ b/SRC/CABLE_MODULES/cable_checks.F90 @@ -627,6 +627,142 @@ RETURN END SUBROUTINE flex_cable_check ! ! NAME +! ML_flex_cable_check +! +! AUTHORS +! Chris Smartt +! +! DESCRIPTION +! check that the parameters defining a multi-layer flex cable are consistent +! This check takes the width, height and offset of two rectangles and checks whether they intersect +! It is used for checking intersection of rows of conductors and conducotrs and dielectric in flex cable models +! If the check fails the cable_spec_error flag is set on return, otherwise it is left unchanged +! +! COMMENTS +! +! +! HISTORY +! +! +! started 12/6/2018 CJS, based on flex_cable_check +! +! +SUBROUTINE ML_flex_cable_check(w1,h1,ox1,oy1,w2,h2,ox2,oy2,check_type,cable_spec_error,cable_name,message) + + real(dp),intent(IN) :: w1,h1 + real(dp),intent(IN) :: ox1,oy1 + real(dp),intent(IN) :: w2,h2 + real(dp),intent(IN) :: ox2,oy2 + + integer,intent(IN) :: check_type + + logical,intent(INOUT) :: cable_spec_error + + character(LEN=line_length),intent(IN) :: cable_name + character(LEN=error_message_length),intent(INOUT) :: message + +! local variables + +logical :: intersect,internal + +real(dp) :: p1,p2,p3,p4 +logical :: internalp3,internalp4 + +logical :: intersectx,internalx,externalx +logical :: intersecty,internaly,externaly + +! START + +if (cable_spec_error) RETURN ! return if an error has already been flagged + +intersect=.FALSE. +internal=.FALSE. +intersectx=.FALSE. +internalx=.FALSE. +externalx=.FALSE. +intersecty=.FALSE. +internaly=.FALSE. +externaly=.FALSE. + +p1=ox1-w1/2.0 +p2=ox1+w1/2.0 +p3=ox2-w2/2.0 +p4=ox2+w2/2.0 + +internalp3=.FALSE. +internalp4=.FALSE. + +if ((p3.GE.p1).AND.(p3.LE.p2)) internalp3=.TRUE. +if ((p4.GE.p1).AND.(p4.LE.p2)) internalp4=.TRUE. + +if (internalp3.AND.internalp4) then + intersectx=.TRUE. + internalx=.TRUE. ! rectangle 2 is inside rectangel1 +else if (internalp3.OR.internalp4) then + intersectx=.TRUE. + internalx=.FALSE. +else +! p3 and p4 are outside the p1-p2 range, we now nead to decide if p1-p2 is inside the range p3-p4 + if ((p3.LT.p1).AND.(p4.GT.p2)) then + intersectx=.TRUE. + externalx=.TRUE. ! rectangle 2 is outside rectangel1 + end if +end if + +p1=oy1-h1/2.0 +p2=oy1+h1/2.0 +p3=oy2-h2/2.0 +p4=oy2+h2/2.0 + +internalp3=.FALSE. +internalp4=.FALSE. + +if ((p3.GE.p1).AND.(p3.LE.p2)) internalp3=.TRUE. +if ((p4.GE.p1).AND.(p4.LE.p2)) internalp4=.TRUE. + +if (internalp3.AND.internalp4) then + intersecty=.TRUE. + internaly=.TRUE. ! rectangle 2 is inside rectangel1 +else if (internalp3.OR.internalp4) then + intersecty=.TRUE. + internaly=.FALSE. +else +! p3 and p4 are outside the p1-p2 range, we now nead to decide if p1-p2 is inside the range p3-p4 + if ((p3.LT.p1).AND.(p4.GT.p2)) then + intersecty=.TRUE. + externaly=.TRUE. ! rectangle 2 is outside rectangel1 + end if +end if + +if (intersectx.AND.intersecty) intersect=.TRUE. +if (internalx.AND.internaly) internal=.TRUE. +if (externalx.AND.externaly) internal=.TRUE. ! rectangle 2 is outside rectangel1 + +if ( (check_type.EQ.1).AND.(intersect) ) then +! we are checking conductor intersections + write(*,*)'Error in cable:',trim(cable_name) + message='conductors intersect' + write(*,'(A)')trim(message) + cable_spec_error=.TRUE. + RETURN + +end if + +! this is the chek for dielectrc being outside the condustor row, the dielectric is shape 2 when +! the subroutine is called +if ( (check_type.EQ.2).AND.( .NOT.(externalx.AND.externaly) ) ) then + write(*,*)'Error in cable:',trim(cable_name) + message='Dielectric intersects conductors' + write(*,'(A)')trim(message) + cable_spec_error=.TRUE. + RETURN +end if + +RETURN + +END SUBROUTINE ML_flex_cable_check +! +! NAME ! dielectric_check ! ! AUTHORS diff --git a/SRC/CABLE_MODULES/cable_module.F90 b/SRC/CABLE_MODULES/cable_module.F90 index a8e33e2..1782e01 100644 --- a/SRC/CABLE_MODULES/cable_module.F90 +++ b/SRC/CABLE_MODULES/cable_module.F90 @@ -92,6 +92,10 @@ ! flex_cable.F90: SUBROUTINE flex_cable_set_internal_domain_information ! flex_cable.F90: SUBROUTINE flex_cable_plot ! +! ML_flex_cable.F90: SUBROUTINE ML_flex_cable_set_parameters +! ML_flex_cable.F90: SUBROUTINE ML_flex_cable_set_internal_domain_information +! ML_flex_cable.F90: SUBROUTINE ML_flex_cable_plot +! ! Dconnector.F90: SUBROUTINE Dconnector_set_parameters ! Dconnector.F90: SUBROUTINE Dconnector_set_internal_domain_information ! Dconnector.F90: SUBROUTINE Dconnector_plot @@ -120,6 +124,7 @@ ! Include general conductor impedance model 12/05/2016 CJS ! put the external condcutor model information into its own structure 6/10/2016 CJS ! 16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions +! 16/3/2018 CJS Add ML_flex_cable ! MODULE cable_module @@ -171,10 +176,12 @@ TYPE:: external_conductor_model real(dp) :: conductor_width2 real(dp) :: conductor_height real(dp) :: conductor_ox + real(dp) :: conductor_oy real(dp) :: dielectric_radius real(dp) :: dielectric_width real(dp) :: dielectric_height real(dp) :: dielectric_ox + real(dp) :: dielectric_oy type(Sfilter) :: dielectric_epsr END TYPE external_conductor_model @@ -261,6 +268,7 @@ integer,parameter :: cable_geometry_type_twinax =7 integer,parameter :: cable_geometry_type_flex_cable =8 integer,parameter :: cable_geometry_type_Dconnector =9 integer,parameter :: cable_geometry_type_ground_plane =10 +integer,parameter :: cable_geometry_type_ML_flex_cable =11 integer,parameter :: impedance_model_type_PEC =0 integer,parameter :: impedance_model_type_cylindrical_with_conductivity =1 @@ -280,6 +288,7 @@ include 'shielded_twisted_pair.F90' include 'spacewire.F90' include 'twinax.F90' include 'flex_cable.F90' +include 'ML_flex_cable.F90' include 'Dconnector.F90' include 'overshield.F90' include 'ground_plane.F90' @@ -329,10 +338,12 @@ include 'conductor_impedance_model.F90' ext_model%conductor_width2=0d0 ext_model%conductor_height=0d0 ext_model%conductor_ox=0d0 + ext_model%conductor_oy=0d0 ext_model%dielectric_radius=0d0 ext_model%dielectric_width=0d0 ext_model%dielectric_height=0d0 ext_model%dielectric_ox=0d0 + ext_model%dielectric_oy=0d0 ext_model%dielectric_epsr=1d0 END SUBROUTINE reset_external_conductor_model @@ -372,6 +383,7 @@ END SUBROUTINE reset_external_conductor_model integer :: matrix_dimension logical :: file_exists + character(len=line_length) :: line ! START @@ -502,11 +514,27 @@ END SUBROUTINE reset_external_conductor_model read(file_unit,*) cable%external_model(i)%conductor_width read(file_unit,*) cable%external_model(i)%conductor_width2 read(file_unit,*) cable%external_model(i)%conductor_height - read(file_unit,*) cable%external_model(i)%conductor_ox + +! format change. This way of reading allows backward compatibility + read(file_unit,'(A)')line + read(line,*,err=100) cable%external_model(i)%conductor_ox,cable%external_model(i)%conductor_oy + GOTO 110 +100 read(line,*) cable%external_model(i)%conductor_ox + cable%external_model(i)%conductor_oy=0d0 +110 CONTINUE + read(file_unit,*) cable%external_model(i)%dielectric_radius read(file_unit,*) cable%external_model(i)%dielectric_width read(file_unit,*) cable%external_model(i)%dielectric_height - read(file_unit,*) cable%external_model(i)%dielectric_ox + +! format change. This way of reading allows backward compatibility + read(file_unit,'(A)')line + read(line,*,err=200) cable%external_model(i)%dielectric_ox,cable%external_model(i)%dielectric_oy + GOTO 210 +200 read(line,*) cable%external_model(i)%dielectric_ox + cable%external_model(i)%dielectric_oy=0d0 +210 CONTINUE + CALL read_Sfilter(cable%external_model(i)%dielectric_epsr,cable_file_unit) end do @@ -677,11 +705,13 @@ END SUBROUTINE reset_external_conductor_model write(file_unit,*) cable%external_model(i)%conductor_width ,' conductor_width ' write(file_unit,*) cable%external_model(i)%conductor_width2 ,' conductor_width2 ' write(file_unit,*) cable%external_model(i)%conductor_height ,' conductor_height ' - write(file_unit,*) cable%external_model(i)%conductor_ox ,' conductor_ox ' + write(file_unit,*) cable%external_model(i)%conductor_ox,cable%external_model(i)%conductor_oy & + ,' conductor_ox, conductor_oy ' write(file_unit,*) cable%external_model(i)%dielectric_radius ,' dielectric_radius' write(file_unit,*) cable%external_model(i)%dielectric_width ,' dielectric_width ' write(file_unit,*) cable%external_model(i)%dielectric_height ,' dielectric_height' - write(file_unit,*) cable%external_model(i)%dielectric_ox ,' dielectric_ox ' + write(file_unit,*) cable%external_model(i)%dielectric_ox,cable%external_model(i)%dielectric_oy & + ,' dielectric_ox, dielectric_oy ' CALL write_Sfilter(cable%external_model(i)%dielectric_epsr,cable_file_unit) end do diff --git a/SRC/CABLE_MODULES/flex_cable.F90 b/SRC/CABLE_MODULES/flex_cable.F90 index fafac68..aef51fb 100644 --- a/SRC/CABLE_MODULES/flex_cable.F90 +++ b/SRC/CABLE_MODULES/flex_cable.F90 @@ -249,7 +249,7 @@ IMPLICIT NONE do i=1,cable%tot_n_conductors write(conductor_string,'(I2)')i cable%conductor_label(i)='Cable name: '//trim(cable%cable_name)// & - '. type: '//trim(cable%cable_type_string)//'. conductor '//conductor_string//' : flex_cable conductor' + '. type: '//trim(cable%cable_type_string)//'. conductor '//conductor_string//' : original_flex_cable conductor' end do END SUBROUTINE flex_cable_set_internal_domain_information diff --git a/SRC/GENERAL_MODULES/constants.F90 b/SRC/GENERAL_MODULES/constants.F90 index 22d236b..4fdb8a5 100644 --- a/SRC/GENERAL_MODULES/constants.F90 +++ b/SRC/GENERAL_MODULES/constants.F90 @@ -104,6 +104,8 @@ IMPLICIT NONE real(dp) :: Laplace_surface_mesh_constant=3d0 ! the mesh edge length on boundaries is calculated as radius/Laplace_surface_mesh_constant ! if Laplace_surface_mesh_constant=5 we have just over 30 elements on the circumference + + real(dp) :: max_mesh_edge_length=1d30 ! the maximum mesh edge length on internal boundaries real(dp) :: Twisted_pair_equivalent_radius=1.5d0 ! The twisted pair commmon mode interaction is calculated by treating the ! common mode as being carried on an 'equivalent cylindrical conductor' diff --git a/SRC/GENERAL_MODULES/general_module.F90 b/SRC/GENERAL_MODULES/general_module.F90 index 15d61a9..c0ed5ac 100644 --- a/SRC/GENERAL_MODULES/general_module.F90 +++ b/SRC/GENERAL_MODULES/general_module.F90 @@ -132,6 +132,8 @@ integer,parameter :: Pspice =3 integer,parameter :: circle =1 integer,parameter :: rectangle=2 integer,parameter :: Dshape=3 +integer,parameter :: semicircle=4 +integer,parameter :: semicircle_diameter=5 integer :: spice_version @@ -168,6 +170,13 @@ logical :: verbose=.FALSE. logical :: plot_propagation_correction_filter_fit_data=.FALSE. +logical :: direct_solver=.FALSE. + +logical :: inf_gnd=.TRUE. + +!logical :: use_ABC=.TRUE. +logical :: use_ABC=.FALSE. + ! There are some differences between windows and unix relating to ! file separators and making directories. The different formats for ! the two ooperating systems supported are included below diff --git a/SRC/GENERAL_MODULES/generate_shapes.F90 b/SRC/GENERAL_MODULES/generate_shapes.F90 index 4d7fb97..ebdfdd5 100644 --- a/SRC/GENERAL_MODULES/generate_shapes.F90 +++ b/SRC/GENERAL_MODULES/generate_shapes.F90 @@ -30,6 +30,7 @@ !SUBROUTINE generate_rectangle_points !SUBROUTINE generate_Dshape_points !SUBROUTINE generate_arc_points +!SUBROUTINE generate_semicircle_points ! !SUBROUTINE generate_circle_points ! @@ -545,3 +546,76 @@ IMPLICIT NONE RETURN END SUBROUTINE generate_line_points +! +!SUBROUTINE generate_semicircle_points +! +! NAME +! SUBROUTINE generate_semicircle_points +! +! DESCRIPTION +! write a semicircle with specified x,y centre and radius to file for plotting with gnuplot +! the semicircle is in the region y>=ycentre +! +! +! COMMENTS +! return the extent of the plotting area... +! +! HISTORY +! +! started 10/05/2013 CJS +! +! +SUBROUTINE generate_semicircle_points(npts,shape_x,shape_y,x,y,r) + +USE type_specifications +USE constants + +IMPLICIT NONE + + integer,intent(OUT) :: npts + real(dp),allocatable,intent(INOUT) :: shape_x(:) + real(dp),allocatable,intent(INOUT) :: shape_y(:) + real(dp),intent(IN) :: x,y,r ! centre x and y coordinates and radius + +! local variables + + real(dp) x1,y1 + real(dp) x2,y2 + real(dp) x3,y3 + + integer :: loop + +! START + +! write the semicircle as two arcs and a straight line to close the shape +! the straght line is formed by the first and last points of the arc + + x1=x+r + y1=y + + x2=x + y2=y+r + + x3=x-r + y3=y + + do loop=1,2 ! first pass to count the points, second pass to set the point coordinates + + npts=0 + + CALL generate_arc_points(npts,shape_x,shape_y,loop,x,y,x1,y1,x2,y2) + + CALL generate_arc_points(npts,shape_x,shape_y,loop,x,y,x2,y2,x3,y3) + + CALL generate_line_points(npts,shape_x,shape_y,loop,x3,y3,x1,y1) + + if (loop.EQ.1) then + ALLOCATE( shape_x(1:npts) ) + ALLOCATE( shape_y(1:npts) ) + end if + + end do ! next loop + + RETURN + +END SUBROUTINE generate_semicircle_points diff --git a/SRC/GENERAL_MODULES/plot_geometry.F90 b/SRC/GENERAL_MODULES/plot_geometry.F90 index 95058ff..42db2d0 100644 --- a/SRC/GENERAL_MODULES/plot_geometry.F90 +++ b/SRC/GENERAL_MODULES/plot_geometry.F90 @@ -29,6 +29,7 @@ !SUBROUTINE write_circle !SUBROUTINE write_rectangle !SUBROUTINE write_Dshape +!SUBROUTINE write_semicircle ! !SUBROUTINE write_circle ! @@ -220,3 +221,67 @@ IMPLICIT NONE RETURN END SUBROUTINE write_Dshape +! +!SUBROUTINE write_semicircle +! +! NAME +! SUBROUTINE write_semicircle +! +! DESCRIPTION +! write a semicircle with specified x,y centre and radius to file for plotting with gnuplot +! the semicircle is in the region y>=0 +! +! +! COMMENTS +! return the extent of the plotting area... +! +! HISTORY +! +! started 10/05/2013 CJS +! using generate_shapes.F90 20/4/2017 CJS +! +! +SUBROUTINE write_semicircle(x,y,r,unit,xmin,xmax,ymin,ymax) + +USE type_specifications +USE constants + +IMPLICIT NONE + + real(dp),intent(IN) :: x,y,r ! centre x and y coordinates and radius + real(dp),intent(INOUT) :: xmin,xmax,ymin,ymax ! extent of the plotting area + integer, intent(IN) :: unit ! unit to write to + +! local variables + + integer :: npts + real(dp),allocatable :: shape_x(:) + real(dp),allocatable :: shape_y(:) + + integer :: i + +! START + + CALL generate_semicircle_points(npts,shape_x,shape_y,x,y,r) + + do i=1,npts + + write(unit,8000)shape_x(i),shape_y(i) +8000 format (4E14.6) + + xmin=min(xmin,shape_x(i)) + xmax=max(xmax,shape_x(i)) + ymin=min(ymin,shape_y(i)) + ymax=max(ymax,shape_y(i)) + + end do + + write(unit,*) + write(unit,*) + + DEALLOCATE( shape_x ) + DEALLOCATE( shape_y ) + + RETURN + +END SUBROUTINE write_semicircle diff --git a/SRC/MATHS_MODULES/cmatrix.F90 b/SRC/MATHS_MODULES/cmatrix.F90 index 2f768bd..c94cffb 100644 --- a/SRC/MATHS_MODULES/cmatrix.F90 +++ b/SRC/MATHS_MODULES/cmatrix.F90 @@ -32,6 +32,7 @@ ! SUBROUTINE write_cmatrix_abs(Mat,dim,unit) ! SUBROUTINE cinvert_Gauss_Jordan(A,n,AI,dim) ! SUBROUTINE c_condition_number (A,n,condition_number,dim) +! ! 16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions ! ! NAME @@ -454,3 +455,5 @@ IMPLICIT NONE RETURN END SUBROUTINE c_condition_number + + diff --git a/SRC/MATHS_MODULES/dmatrix.F90 b/SRC/MATHS_MODULES/dmatrix.F90 index b2ea97e..d08202a 100644 --- a/SRC/MATHS_MODULES/dmatrix.F90 +++ b/SRC/MATHS_MODULES/dmatrix.F90 @@ -28,6 +28,7 @@ ! SUBROUTINE dwrite_matrix(a,ar,ac,dim,unit) ! SUBROUTINE dread_matrix(a,ar,ac,dim,unit) ! SUBROUTINE dinvert_Gauss_Jordan(A,ar,AI,dim) +! !_____________________________________________________________________ ! ! NAME @@ -283,3 +284,4 @@ IMPLICIT NONE RETURN END SUBROUTINE dinvert_Gauss_Jordan + diff --git a/SRC/PUL_PARAMETER_CALCULATION/Aprod.F90 b/SRC/PUL_PARAMETER_CALCULATION/Aprod.F90 new file mode 100644 index 0000000..cb940c1 --- /dev/null +++ b/SRC/PUL_PARAMETER_CALCULATION/Aprod.F90 @@ -0,0 +1,207 @@ +! This file is part of SACAMOS, State of the Art CAble MOdels in Spice. +! It was developed by the University of Nottingham and the Netherlands Aerospace +! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK. +! +! Copyright (C) 2016-2017 University of Nottingham +! +! SACAMOS is free software: you can redistribute it and/or modify it under the +! terms of the GNU General Public License as published by the Free Software +! Foundation, either version 3 of the License, or (at your option) any later +! version. +! +! SACAMOS is distributed in the hope that it will be useful, but +! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! A copy of the GNU General Public License version 3 can be found in the +! file GNU_GPL_v3 in the root or at . +! +! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to +! the GNU Lesser General Public License. A copy of the GNU Lesser General Public +! License version can be found in the file GNU_LGPL in the root of EISPACK +! (/SRC/EISPACK ) or at . +! +! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk +! +! File Contents: +! +! SUBROUTINE Aprod +! SUBROUTINE ATprod +! SUBROUTINE zAprod +! SUBROUTINE zATprod +! +! NAME +! SUBROUTINE Aprod +! +! DESCRIPTION +! Matrix vector multiplication for sparse, real matrices +! +! COMMENTS +! +! +! HISTORY +! 19/3/2018 CJS +! +! +SUBROUTINE Aprod ( n, x, y ) + +! form the vector y=K*x + +USE type_specifications + +IMPLICIT NONE + +integer n +real(dp) :: x(n), y(n) + +integer :: i,row,col + +! START + +y(:)=0d0 + +do i=1,n_entry + row=i_K(i) + col=j_K(i) + +! if (row.GT.n) then +! write(*,*)'Error in Aprod, row=',row,' n=',n +! end if +! if (col.GT.n) then +! write(*,*)'Error in Aprod, col=',col,' n=',n +! end if + + y(row)=y(row)+s_K_re(i)*x(col) + +end do + +RETURN + +END SUBROUTINE Aprod +! +! NAME +! SUBROUTINE ATprod +! +! DESCRIPTION +! Transpose Matrix vector multiplication for sparse, real matrices +! +! COMMENTS +! +! +! HISTORY +! 19/3/2018 CJS +! +! +SUBROUTINE ATprod ( n, x, y ) + +! form the vector y=AT*x + +USE type_specifications + +IMPLICIT NONE + +integer n +real(dp) :: x(n), y(n) + +integer :: i,row,col + +! START + +y(:)=0d0 + +do i=1,n_entry + row=j_K(i) + col=i_K(i) + + y(row)=y(row)+s_K_re(i)*x(col) + +end do + +RETURN + +END SUBROUTINE ATprod +! +! NAME +! SUBROUTINE zAprod +! +! DESCRIPTION +! Matrix vector multiplication for sparse, real matrices +! +! COMMENTS +! +! +! HISTORY +! 19/3/2018 CJS +! +! +SUBROUTINE zAprod ( n, x, y ) + +! form the vector y=A*x + +USE type_specifications + +IMPLICIT NONE + +integer n +complex(dp) :: x(n), y(n) + +integer :: i,row,col + +! START + +y(:)=0d0 + +do i=1,n_entry + row=i_K(i) + col=j_K(i) + + y(row)=y(row)+s_K(i)*x(col) + +end do + +RETURN + +END SUBROUTINE zAprod +! +! NAME +! SUBROUTINE zATprod +! +! DESCRIPTION +! Transpose Matrix vector multiplication for sparse, real matrices +! +! COMMENTS +! +! +! HISTORY +! 19/3/2018 CJS +! +! +SUBROUTINE zATprod ( n, x, y ) + +! form the vector y=AT*x + +USE type_specifications + +IMPLICIT NONE + +integer n +complex(dp) :: x(n), y(n) + +integer :: i,row,col + +! START + +y(:)=0d0 + +do i=1,n_entry + row=j_K(i) + col=i_K(i) + + y(row)=y(row)+(s_K(i))*x(col) + +end do + +RETURN + +END SUBROUTINE zATprod diff --git a/SRC/PUL_PARAMETER_CALCULATION/CG_solve.F90 b/SRC/PUL_PARAMETER_CALCULATION/CG_solve.F90 new file mode 100644 index 0000000..8da0821 --- /dev/null +++ b/SRC/PUL_PARAMETER_CALCULATION/CG_solve.F90 @@ -0,0 +1,382 @@ +! +! This file is part of SACAMOS, State of the Art CAble MOdels in Spice. +! It was developed by the University of Nottingham and the Netherlands Aerospace +! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK. +! +! Copyright (C) 2016-2017 University of Nottingham +! +! SACAMOS is free software: you can redistribute it and/or modify it under the +! terms of the GNU General Public License as published by the Free Software +! Foundation, either version 3 of the License, or (at your option) any later +! version. +! +! SACAMOS is distributed in the hope that it will be useful, but +! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! A copy of the GNU General Public License version 3 can be found in the +! file GNU_GPL_v3 in the root or at . +! +! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to +! the GNU Lesser General Public License. A copy of the GNU Lesser General Public +! License version can be found in the file GNU_LGPL in the root of EISPACK +! (/SRC/EISPACK ) or at . +! +! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk +! +! SUBROUTINE solve_real_symm(n, b, x,tol) +! SUBROUTINE dot ( n, x, y ,res) +! SUBROUTINE solve_complex_symm(n, b, x,tol) +! SUBROUTINE zdot ( n, x, y ,res) +! +! NAME +! solve_real_symm +! +! DESCRIPTION +! +! Solve the system Ax=b using a conjugate gradient method. +! The iterative solution halts when the norm of the residual is less than tol +! +! HISTORY +! +! started 19/6/2018 CJS +! +! COMMENTS +! +SUBROUTINE solve_real_symm(n, b, x,tol) + +USE type_specifications + +IMPLICIT NONE + +! variables passed to subroutine + +integer,intent(IN) :: n + +real(dp),intent(IN) :: b(n) +real(dp),intent(INOUT) :: x(n) ! x contains the initial guess for x +real(dp),intent(IN) :: tol + +! local variables + +real(dp) :: ak +real(dp) :: rk(n),rkb(n) +real(dp) :: pk(n),pkb(n) +real(dp) :: bk + +real(dp) :: t1(n),t2(n) +real(dp) :: res(n),res_norm +real(dp) :: num,den + +integer :: i,ii,itloop + +logical :: minres + +! START + +write(*,*)'CALLED solve_real_symm with n=',n + +! choice between normal conjugate gradient and minimum residual forms. + +minres=.TRUE. + +! set starting values + +! calculate the residual res=b-A.xk + CALL Aprod( n, x, t1 ) + res(:)=b(:)-t1(:) + +! calculate the size of the residual + res_norm=0d0 + do i=1,n + res_norm=res_norm+res(i)*res(i) + end do + res_norm=sqrt(res_norm) + + write(*,*)'iteration ',0,' norm=',res_norm + +! initial r1=b-A.x + CALL Aprod( n, x, t1 ) + rk(:)=b(:)-t1(:) + + if (minres) then +! For the minimum residual algorithm r1b=A.r1 + CALL Aprod( n, rk, rkb ) + else +! For the normal algorithm r1b=r1 + rkb(:)=rk(:) + end if + +! p1=r1 + pk=rk + +! p1b=r1b + pkb=rkb + +! iteration loop + do itloop=1,2*n + +! calculate ak=(rkb.rk)/(pkb.A.pk) + + CALL dot(n,rkb,rk,num) + + CALL Aprod( n, pk, t1 ) + CALL dot(n,pkb,t1,den) + + ak=num/den + +! calculate the next value of the solution vector, xk+1=xk+ak*pk + x(:)=x(:)+ak*pk(:) + +! calculate the residual res=b-A.xk + CALL Aprod( n, x, t2 ) + res(:)=b(:)-t2(:) + +! calculate the size of the residual + res_norm=0d0 + do i=1,n + res_norm=res_norm+res(i)*res(i) + end do + res_norm=sqrt(res_norm) + + if (res_norm.LT.tol) then + write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak + RETURN + end if + +! calculate rk+1=rk-ak*A.pk + rk(:) =rk(:) -ak*t1(:) + +! calculate rk+1b=rkb-ak*AT.pkb, note for the symmetric matrices here AT=A + CALL ATprod( n, pkb, t2 ) + rkb(:)=rkb(:)-ak*t2(:) + +! calculate bk=(rk+1b.rk+1)/(rkb.rk) + den=num + CALL dot(n,rkb,rk,num) + bk=num/den + + write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak,' bk=',bk + +! calculate pk+1=rk+1-bk*pk + pk(:)=rk(:)+bk*pk(:) + +! calculate pk+1b=rk+1b-bk*pkb + pkb(:)=rkb(:)+bk*pkb(:) + + end do ! next iteration + +! finish + + RETURN + + +END SUBROUTINE solve_real_symm +! +! _________________________________________________________________________ +! +! + +SUBROUTINE dot ( n, x, y ,res) + +! form the dot product of x and y + +USE type_specifications + +IMPLICIT NONE + +integer n +real(dp) :: x(n), y(n), res + +integer :: i + +! START + +res=0d0 + +do i=1,n + + res=res+x(i)*y(i) + +end do + +RETURN + +END SUBROUTINE dot +! +! NAME +! solve_complex_symm +! +! DESCRIPTION +! +! Solve the system Ax=b using a conjugate gradient method. +! The iterative solution halts when the norm of the residual is less than tol +! +! HISTORY +! +! started 19/6/2018 CJS +! +! COMMENTS +! + SUBROUTINE solve_complex_symm(n, b, x,tol) + +USE type_specifications + +IMPLICIT NONE + +! variables passed to subroutine + +integer,intent(IN) :: n + +complex(dp),intent(IN) :: b(n) +complex(dp),intent(INOUT) :: x(n) ! x contains the initial guess for x +real(dp),intent(IN) :: tol + +! local variables + +complex(dp) :: ak +complex(dp) :: rk(n),rkb(n) +complex(dp) :: pk(n),pkb(n) +complex(dp) :: bk + +complex(dp) :: t1(n),t2(n) +complex(dp) :: res(n) +real(dp) :: res_norm +complex(dp) :: num,den + +integer :: i,ii,itloop + +logical :: minres + +! START + +write(*,*)'CALLED solve_real_symm with n=',n + +! choice between normal conjugate gradient and minimum residual forms. + +minres=.TRUE. + +! set starting values + +! calculate the residual res=b-A.xk + CALL zAprod( n, x, t1 ) + res(:)=b(:)-t1(:) + +! calculate the size of the residual + res_norm=0d0 + do i=1,n + res_norm=res_norm+res(i)*res(i) + end do + res_norm=sqrt(res_norm) + + write(*,*)'iteration ',0,' norm=',res_norm + +! initial r1=b-A.x + CALL zAprod( n, x, t1 ) + rk(:)=b(:)-t1(:) + + if (minres) then +! For the minimum residual algorithm r1b=A.r1 + CALL zAprod( n, rk, rkb ) + else +! For the normal algorithm r1b=r1 + rkb(:)=rk(:) + end if + +! p1=r1 + pk=rk + +! p1b=r1b + pkb=rkb + +! iteration loop + do itloop=1,2*n + +! calculate ak=(rkb.rk)/(pkb.A.pk) + + CALL zdot(n,rkb,rk,num) + + CALL zAprod( n, pk, t1 ) + CALL zdot(n,pkb,t1,den) + + ak=num/den + +! calculate the next value of the solution vector, xk+1=xk+ak*pk + x(:)=x(:)+ak*pk(:) + +! calculate the residual res=b-A.xk + CALL zAprod( n, x, t2 ) + res(:)=b(:)-t2(:) + +! calculate the size of the residual + res_norm=0d0 + do i=1,n + res_norm=res_norm+abs(res(i))**2 + end do + res_norm=sqrt(res_norm) + + if (res_norm.LT.tol) then + write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak + RETURN + end if + +! calculate rk+1=rk-ak*A.pk + rk(:) =rk(:) -ak*t1(:) + +! calculate rk+1b=rkb-ak*AT.pkb, note for the symmetric matrices here AT=A + CALL zATprod( n, pkb, t2 ) + rkb(:)=rkb(:)-ak*t2(:) + +! calculate bk=(rk+1b.rk+1)/(rkb.rk) + den=num + CALL zdot(n,rkb,rk,num) + bk=num/den + + write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak,' bk=',bk + +! calculate pk+1=rk+1-bk*pk + pk(:)=rk(:)+bk*pk(:) + +! calculate pk+1b=rk+1b-bk*pkb + pkb(:)=rkb(:)+bk*pkb(:) + + end do ! next iteration + +! finish + + RETURN + + +END SUBROUTINE solve_complex_symm +! +! _________________________________________________________________________ +! +! + +SUBROUTINE zdot ( n, x, y ,res) + +! form the dot product of x and y + +USE type_specifications + +IMPLICIT NONE + +integer n +complex(dp) :: x(n), y(n), res + +integer :: i + +! START + +res=0d0 + +do i=1,n + + res=res+x(i)*y(i) ! note no complex conjugate in the inner product here. This gives the biconjugate gradient method. + +end do + +RETURN + +END SUBROUTINE zdot diff --git a/SRC/PUL_PARAMETER_CALCULATION/Laplace.F90 b/SRC/PUL_PARAMETER_CALCULATION/Laplace.F90 index 85e1f98..efdbdd9 100644 --- a/SRC/PUL_PARAMETER_CALCULATION/Laplace.F90 +++ b/SRC/PUL_PARAMETER_CALCULATION/Laplace.F90 @@ -106,6 +106,7 @@ ! started 5/7/2016 CJS. This subroutineis based on software from NLR and has been ! translated from Matlab to Fortran. ! 8/5/2017 CJS: Include references to Theory_Manual +! 19/6/2018 CJS: Include iterative solver based on the biconjugate gradient method ! SUBROUTINE Laplace(mesh_filename,dim,first_surface_is_free_space_boundary, & n_dielectric_regions,dielectric_region_epsr,frequency,Lmat,Cmat,Gmat,ox,oy) @@ -169,7 +170,7 @@ integer :: N_nodes_in ! number of nodes in t integer :: N_elements_in ! number of elements in the gmsh file (this is not necessarily the number of triangular elements in the FE mesh) integer :: n_boundaries ! number of boundaries, not including dielectric (internal) boundaries -integer :: N_boundaries_max ! maximum boundary number generated by PUL_LC_Laplace and found in the mesh file +integer :: N_boundaries_max ! maximum boundary number generated by PUL_LC_Laplace and found insolve_real_symm(n, b, x,tol) the mesh file integer :: N_boundary ! number of viable boundaries in the mesh i.e. boundaries with two or more points on integer :: boundary_number integer :: N_elements_boundary_temp @@ -202,7 +203,7 @@ integer :: N_unknown ! number of known node voltages i.e integer :: N_known ! number of unknown node voltages integer :: jmax ! maximum boundary number (i.e. number of conductors) - ! this is used to determine the number of finite element solutions (right hand sides to solve for) to + ! this is used to determine the number of finite element solutions (right hand solve_real_symm(n, b, x,tol)sides to solve for) to ! fill the capacitance/conductance matrix integer :: total_n_boundary_nodes ! total number of boundary nodes i.e. the number of known node values @@ -216,10 +217,12 @@ complex(dp),allocatable :: c(:,:) ! element based array of constants re complex(dp),allocatable :: delta(:) ! element based array with a value related to the element geometry complex(dp),allocatable :: eps_r(:) ! element based array with the complex relative permittivity of the element -! 1D arrays used in the construction of the K matrix ( K(i_K(:),j_K(:))=K(i_K(:),j_K(:))+s_K(:) ) -integer,allocatable :: i_K(:) -integer,allocatable :: j_K(:) -complex(dp),allocatable :: s_K(:) +! +!integer :: n_entrysolve_real_symm(n, b, x,tol) +!! 1D arrays used in the construction of the K matrix ( K(i_K(:),j_K(:))=K(i_K(:),j_K(:))+s_K(:) ) +!integer,allocatable :: i_K(:) +!integer,allocatable :: j_K(:) +!complex(dp),allocatable :: s_K(:) ! 1D arrays used in the construction of the right hand side vectors ( K_rhs(i_K_rhs(:),j_K_rhs(:))=K_rhs(i_K_rhs(:),j_K_rhs(:))+s_K_rhs(:) ) integer,allocatable :: i_K_rhs(:) ! row number @@ -248,7 +251,6 @@ real(dp),allocatable :: Conductance_energy(:,:) integer :: node,new_node integer :: element -integer :: n_entry integer :: n_entry_K_rhs integer :: n1,n2,n3 complex(dp) :: x1,x2,x3,y1,y2,y3 @@ -280,6 +282,37 @@ integer :: nr integer :: ierr ! error code for matrix inversion +! variables for iterative solver + +logical checkA +integer itnlim +real(dp) rtol +integer nout +logical goodb +logical precon +real(dp) shift + +integer istop, itn +real(dp) anorm, acond, rnorm, ynorm + +real(dp),allocatable :: r1(:) +real(dp),allocatable :: r2(:) +real(dp),allocatable :: vt(:) +real(dp),allocatable :: wt(:) +real(dp),allocatable :: yt(:) +real(dp),allocatable :: bt(:) +real(dp),allocatable :: xt(:) + +logical :: lossy_dielectric ! flag to indicate lossy dielectric i.e. we must solve a complex problem + +logical :: wantse +integer :: n,m +real(dp) :: atol, btol, conlim, damp +real(dp) :: Arnorm, xnorm +real(dp),allocatable :: se(:) + +real(dp) :: Vout + ! START if (verbose) write(*,*)'CALLED: Laplace' @@ -542,12 +575,16 @@ if (verbose) write(*,*)' Number of materials=',N_materials ALLOCATE ( Mat_prop(1:N_materials,1:4) ) +lossy_dielectric=.FALSE. + do i=1,N_materials Mat_prop(i,1)=i ! material number Mat_prop(i,2)=dble(dielectric_region_epsr(i-1)) ! Re{epsr} Mat_prop(i,3)=aimag(dielectric_region_epsr(i-1)) ! Im{epsr} Mat_prop(i,4)=frequency + if ( abs(Mat_prop(i,3)).GT.small ) lossy_dielectric=.TRUE. + if (verbose) then count=0 do element=1,N_elements @@ -942,7 +979,16 @@ do loop=1,2 ls=sqrt((x2-x1)**2+(y2-y1)**2) ! l2= boundary element edge length - gamma=-eps_r(el)/(log(rho)*rho) ! gamma (see equation 4.3 with equation 4.93) + if (use_ABC) then +! This implements an asymptotic boundary condition. This apears to cause convergence +! issues for the Capacitance matrix calculation so should be used with caution. +! The default is the Neumann boundary condition (gamma=0) + gamma=-eps_r(el)/(log(rho)*rho) ! gamma (see equation 4.3 with equation 4.93) + else +! Neumann boundary condition on the outer boundary. This effectively sets the charge to zero on the outer +! boundary and the convergence of the capacitance matrix is improved compared with the ABC + gamma=0d0 + end if ! K11 contributions according to equation 4.51 (note delta_ij=1 if i=j, 0 otherwise) n_entry=n_entry+1 @@ -1002,8 +1048,6 @@ ALLOCATE ( v_tmp(1:N_known) ) ! solution based on a full matrix inverse -ALLOCATE ( K(1:N_unknown,1:N_unknown) ) -ALLOCATE ( KI(1:N_unknown,1:N_unknown) ) ALLOCATE ( K_rhs(1:N_unknown,1:N_known) ) if(verbose) then @@ -1011,12 +1055,6 @@ if(verbose) then write(*,*)'Number of entries in K_rhs',n_entry_K_rhs end if -! STAGE 11a: fill the K matrix -K(1:N_unknown,1:N_unknown)=0d0 -do i=1,n_entry - K(i_K(i),j_K(i))=K(i_K(i),j_K(i))+s_K(i) -end do - ! STAGE 11b: fill the K_rhs matrix K_rhs(1:N_unknown,1:N_known)=0d0 do i=1,n_entry_K_rhs @@ -1030,20 +1068,108 @@ else write(*,*)'Dimension of K in Laplace is',N_unknown,N_unknown end if +if (direct_solver) then + ALLOCATE ( K(1:N_unknown,1:N_unknown) ) + ALLOCATE ( KI(1:N_unknown,1:N_unknown) ) + +! STAGE 11a: fill the K matrix + K(1:N_unknown,1:N_unknown)=0d0 + do i=1,n_entry + K(i_K(i),j_K(i))=K(i_K(i),j_K(i))+s_K(i) + end do + ! STAGE 11c: Invert the K matrix -if(verbose) write(*,*)'Invert the K matrix' -ierr=0 ! set ierr=0 to cause an error within cinvert_Gauss_Jordan if there is a problem calculating the inverse -CALL cinvert_Gauss_Jordan(K,N_unknown,KI,N_unknown,ierr) + if(verbose) write(*,*)'Invert the K matrix' + ierr=0 ! set ierr=0 to cause an error within cinvert_Gauss_Jordan if there is a problem calculating the inverse + CALL cinvert_Gauss_Jordan(K,N_unknown,KI,N_unknown,ierr) ! STAGE 11d loop over all the RHS vectors solving the matrix equation -do j1=1,jmax-1 - do j2=j1,jmax-1 + do j1=1,jmax-1 + do j2=j1,jmax-1 + v_tmp(1:n_known)=V(j1,j2,1:n_known) + b_tmp(1:N_unknown)=-matmul(K_rhs,v_tmp) + x_tmp=matmul(KI,b_tmp) + x(j1,j2,1:N_unknown)=x_tmp(1:N_unknown) + end do + end do + +else if(.NOT.lossy_dielectric) then + + checkA = .true. + itnlim = N_unknown * 2 + rtol = 1.0D-12 + nout=6 + goodb=.FALSE. + precon = .FALSE. + shift=0d0 + + allocate(r1(1:N_unknown)) + allocate(r2(1:N_unknown)) + allocate(vt(1:N_unknown)) + allocate(wt(1:N_unknown)) + allocate(yt(1:N_unknown)) + allocate(xt(1:N_unknown)) + allocate(bt(1:N_unknown)) + + ALLOCATE( s_K_re(1:n_entry) ) + s_K_re(1:n_entry)=dble(s_K(1:n_entry)) + +! Iterative solution + do j1=1,jmax-1 + do j2=j1,jmax-1 + v_tmp(1:n_known)=V(j1,j2,1:n_known) b_tmp(1:N_unknown)=-matmul(K_rhs,v_tmp) - x_tmp=matmul(KI,b_tmp) + bt(1:N_unknown)=b_tmp(1:N_unknown) + +! UoN conjugate gradient solution + CALL solve_real_symm(N_unknown, bt, xt, rtol) + + x(j1,j2,1:N_unknown)=xt(1:N_unknown) + + end do + end do + + deallocate(r1) + deallocate(r2) + deallocate(vt) + deallocate(wt) + deallocate(yt) + deallocate(xt) + deallocate(bt) + +else if (lossy_dielectric) then + + itnlim=4*N_unknown + nout=6 + wantse=.FALSE. + atol=1D-8 + btol=1d-8 + conlim=1d10 + damp=0d0 + allocate(se(1:N_unknown)) + +! Iterative solution + do j1=1,jmax-1 + do j2=j1,jmax-1 + + v_tmp(1:n_known)=V(j1,j2,1:n_known) + b_tmp(1:N_unknown)=-matmul(K_rhs,v_tmp) + n=N_unknown + m=N_unknown + +! UoN conjugate gradient solution + rtol = 1.0D-12 + CALL solve_complex_symm(N_unknown, b_tmp, x_tmp, rtol) + x(j1,j2,1:N_unknown)=x_tmp(1:N_unknown) - end do -end do + + end do + end do + + deallocate(se) + +end if ! STAGE 12 Determine the voltage phi in each node of the mesh @@ -1321,6 +1447,7 @@ DEALLOCATE( eps_r ) DEALLOCATE( i_K ) DEALLOCATE( j_K ) DEALLOCATE( s_K ) +if ( ALLOCATED(s_K_re) ) DEALLOCATE( s_K_re ) DEALLOCATE( i_K_rhs ) DEALLOCATE( j_K_rhs ) @@ -1330,8 +1457,8 @@ DEALLOCATE ( x ) DEALLOCATE ( x_tmp ) DEALLOCATE ( b_tmp ) DEALLOCATE ( v_tmp ) -DEALLOCATE ( K ) -DEALLOCATE ( KI ) +if ( ALLOCATED(K) ) DEALLOCATE ( K ) +if ( ALLOCATED(KI) ) DEALLOCATE ( KI ) DEALLOCATE ( K_rhs ) DEALLOCATE( phi ) diff --git a/SRC/PUL_PARAMETER_CALCULATION/Makefile b/SRC/PUL_PARAMETER_CALCULATION/Makefile index 2853426..f546b2c 100644 --- a/SRC/PUL_PARAMETER_CALCULATION/Makefile +++ b/SRC/PUL_PARAMETER_CALCULATION/Makefile @@ -1,5 +1,5 @@ default: $(PUL_PARAMETER_CALCULATION_OBJS) # $(OBJ_MOD_DIR)/%.o: %.F90 $(TYPE_SPEC_MODULE) $(CONSTANTS_MODULE) $(GENERAL_MODULE) \ - PUL_analytic.F90 PUL_LC_Laplace.F90 Laplace.F90 + PUL_analytic.F90 PUL_LC_Laplace.F90 Laplace.F90 Aprod.F90 CG_solve.F90 $(FC) $(FLAGS) -c -o $@ $< diff --git a/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 b/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 index 510dc4d..ec7e4b2 100644 --- a/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 +++ b/SRC/PUL_PARAMETER_CALCULATION/PUL_LC_Laplace.F90 @@ -57,8 +57,9 @@ ! 14/10/2017 CJS make the filter fitting minimum aorder=1, border=0 and ! ensure that border=aorder-1 to make the choice of model order more sensible ! 16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions +! 4/3/2018 CJS Start to include the infinite ground plane option ! -! +! 16/3/2018 CJS add offset to x and y for shapes. Relates to ML_flex cables ! SUBROUTINE PUL_LC_Laplace(PUL,name,fit_order,freq_spec,domain) ! @@ -89,6 +90,8 @@ IMPLICIT NONE ! flags to indicate the configuration logical :: ground_plane + logical :: finite_ground_plane + logical :: infinite_ground_plane logical :: overshield logical :: open_boundary @@ -115,16 +118,20 @@ IMPLICIT NONE ! variables required for generating the gmsh geometry input file: we may be able to simplify ! things here as there is some overlap with the PUL structure integer,allocatable :: shape_type(:) ! holds a nuber which indicates the shape type + integer,allocatable :: loop_number(:) ! holds the loop number for this shape real(dp),allocatable :: x(:) ! x coordinate of the centre of the cable to which this shape belongs real(dp),allocatable :: y(:) ! y coordinate of the centre of the cable to which this shape belongs real(dp),allocatable :: r(:) ! radius of a circular shape or curve radius for a Dshape real(dp),allocatable :: rw(:) ! width1 (x dimension) of rectangular/ Dshape real(dp),allocatable :: rw2(:) ! width2 (x dimension) of rectangular/ Dshape real(dp),allocatable :: rh(:) ! height (y dimension) of rectangular shape - real(dp),allocatable :: ro(:) ! offset in the x direction of this shape from the centre (x():),y(:) above) + real(dp),allocatable :: rox(:) ! offset in the x direction of this shape from the centre (x():),y(:) above) + real(dp),allocatable :: roy(:) ! offset in the y direction of this shape from the centre (x():),y(:) above) real(dp),allocatable :: rtheta(:) ! rotation angle of conductor real(dp),allocatable :: dl(:) ! edge length for the mesh on this shape logical,allocatable :: conductor_has_dielectric(:) ! flag to indicate that a conductor has a dielectric around it + + real(dp) :: dl_local ! edge length for the mesh at rectangle edge centre complex(dp) :: epsr ! relative permittivity value @@ -151,6 +158,7 @@ IMPLICIT NONE real(dp) :: rl,rt,rdl,rx1,rx2,ry1,ry2 real(dp) :: xp,yp,xt,yt integer :: p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12 + real(dp) :: edge_length character(LEN=filename_length) :: command ! string used for the command which runs gmsh from the shell integer :: gmsh_exit_stat ! exit status from the shell command which runs gmsh. @@ -171,6 +179,8 @@ IMPLICIT NONE logical :: dielectric_defined ! flag to indicate whether there is a dielectric logical :: frequency_dependent_dielectric ! flag to indicate whether there is a frequency dependent dielectric + integer :: gpline1,gpline2 ! line segments for infinite ground plane + integer :: row,col ! loop variables for matrix elements integer :: i,ii ! temporary loop variables @@ -181,9 +191,11 @@ IMPLICIT NONE ! STAGE 1: work out the configuration for the calculation i.e. is there a ground plane, is the outer boundary free space or a conductor (overshield) ground_plane=.FALSE. + infinite_ground_plane=.FALSE. + finite_ground_plane=.FALSE. overshield=.FALSE. open_boundary=.FALSE. - + if ((.NOT.PUL%overshield_present).AND.(.NOT.PUL%ground_plane_present)) then ! SOLUTION TYPE 1: NO OVERSHIELD, NO GROUND PLANE @@ -200,10 +212,16 @@ IMPLICIT NONE ! SOLUTION TYPE 2: NO OVERSHIELD, WITH GROUND PLANE ground_plane=.TRUE. + if (inf_gnd) then + infinite_ground_plane=.TRUE. + else + finite_ground_plane=.TRUE. + end if + + tot_n_boundaries_without_dielectric=PUL%n_conductors+1 ! add a free space boundary first_surface_is_free_space_boundary=.TRUE. end_conductor_loop=PUL%n_conductors-1 ! the ground plane is not included here - added separately - end_conductor_to_dielectric_loop=PUL%n_conductors - tot_n_boundaries_without_dielectric=PUL%n_conductors+1 ! add a free space boundary + end_conductor_to_dielectric_loop=PUL%n_conductors first_conductor_shape=2 ! shape 1 is for the outer boundary else @@ -221,6 +239,8 @@ IMPLICIT NONE if (verbose) then write(*,*)'ground_plane:',ground_plane + write(*,*)'infinite_ground_plane:',infinite_ground_plane + write(*,*)'finite_ground_plane:',finite_ground_plane write(*,*)'open_boundary:',open_boundary write(*,*)'overshield:',overshield write(*,*)'first_surface_is_free_space_boundary:',first_surface_is_free_space_boundary @@ -304,13 +324,16 @@ IMPLICIT NONE n_shapes=tot_n_boundaries ALLOCATE( shape_type(1:n_shapes) ) shape_type(1:n_shapes)=circle ! set all boudaries to type circle initially + ALLOCATE( loop_number(1:n_shapes) ) + loop_number(1:n_shapes)=-1 ! set all boudaries to -1 i.e. not set ALLOCATE( x(1:n_shapes) ) ALLOCATE( y(1:n_shapes) ) ALLOCATE( r(1:n_shapes) ) ALLOCATE( rw(1:n_shapes) ) ALLOCATE( rw2(1:n_shapes) ) ALLOCATE( rh(1:n_shapes) ) - ALLOCATE( ro(1:n_shapes) ) + ALLOCATE( rox(1:n_shapes) ) + ALLOCATE( roy(1:n_shapes) ) ALLOCATE( rtheta(1:n_shapes) ) ALLOCATE( dl(1:n_shapes) ) ALLOCATE( shape_to_region(1:n_shapes,1:2) ) @@ -332,17 +355,17 @@ IMPLICIT NONE do i=1,end_conductor_loop ! if the last conductor is the ground plane it is included in x/y,max/min calculation later if (PUL%shape(i).EQ.circle) then - xmin=min(xmin,PUL%x(i)-PUL%r(i)+PUL%o(i)) - xmax=max(xmax,PUL%x(i)+PUL%r(i)+PUL%o(i)) - ymin=min(ymin,PUL%y(i)-PUL%r(i)) - ymax=max(ymax,PUL%y(i)+PUL%r(i)) + xmin=min(xmin,PUL%x(i)-PUL%r(i)+PUL%ox(i)) + xmax=max(xmax,PUL%x(i)+PUL%r(i)+PUL%ox(i)) + ymin=min(ymin,PUL%y(i)-PUL%r(i)+PUL%oy(i)) + ymax=max(ymax,PUL%y(i)+PUL%r(i)+PUL%oy(i)) else if (PUL%shape(i).EQ.rectangle) then - xmin=min(xmin,PUL%x(i)-PUL%rw(i)+PUL%o(i)) - xmax=max(xmax,PUL%x(i)+PUL%rw(i)+PUL%o(i)) - ymin=min(ymin,PUL%y(i)-PUL%rh(i)) - ymax=max(ymax,PUL%y(i)+PUL%rh(i)) + xmin=min(xmin,PUL%x(i)-PUL%rw(i)+PUL%ox(i)) + xmax=max(xmax,PUL%x(i)+PUL%rw(i)+PUL%ox(i)) + ymin=min(ymin,PUL%y(i)-PUL%rh(i)+PUL%oy(i)) + ymax=max(ymax,PUL%y(i)+PUL%rh(i)+PUL%oy(i)) else if (PUL%shape(i).EQ.Dshape) then @@ -352,13 +375,17 @@ IMPLICIT NONE ymin=min(ymin,PUL%y(i)-PUL%rh(i)) ymax=max(ymax,PUL%y(i)+PUL%rh(i)) - else +! note don't check semicircle as this is only for infinite ground only + + else if (PUL%shape(i).NE.semicircle) then + write(*,*)'Unknown shape type',PUL%shape(i) + end if end do ! next conductor if (ground_plane) then -! include the ground plane in the bundle sizing +! include the ground plane in the bundle sizing by adding the origin point ! here we assume that the ground plane is along the x axis. xmin=min(xmin,0d0) @@ -382,21 +409,48 @@ IMPLICIT NONE write(*,*)'boundary radius=',rboundary end if -! set the first shape information to relate to the outer boundary, this is a circle, centred on the bundle centre +! set the first shape information to relate to the outer boundary +! +! If we have a finite gound plane of free space boundary then this is a circle, centred on the bundle centre ! note that for the Laplace calculation the geometry will be shifted so that the origin is at xc,yc so the ! free space outer boundary is centred on the origin. This is required for the free space boundary condition in Laplace to work correctly. - x(1)=xc - y(1)=yc - r(1)=rboundary - rw(1)=0d0 - rw2(1)=0d0 - rh(1)=0d0 - ro(1)=0d0 - rtheta(1)=0d0 - dl(1)=r(1)/(2*Laplace_surface_mesh_constant) - shape_to_region(1,inside) =0 ! inside is the dielectric region - shape_to_region(1,outside)=-1 ! no outside region + if (.NOT.infinite_ground_plane) then + + x(1)=xc + y(1)=yc + r(1)=rboundary + rw(1)=0d0 + rw2(1)=0d0 + rh(1)=0d0 + rox(1)=0d0 + roy(1)=0d0 + rtheta(1)=0d0 + dl(1)=r(1)/(2*Laplace_surface_mesh_constant) + shape_to_region(1,inside) =0 ! inside is the dielectric region + shape_to_region(1,outside)=-1 ! no outside region + + else +! +! If we have an infinite ground plane then we have a semicircle centred on the origin + + shape_type(1)=semicircle + xc=0d0 + yc=0d0 + x(1)=xc + y(1)=yc + r(1)=rboundary + rw(1)=0d0 + rw2(1)=0d0 + rh(1)=0d0 + rox(1)=0d0 + roy(1)=0d0 + rtheta(1)=0d0 + dl(1)=r(1)/(2*Laplace_surface_mesh_constant) + shape_to_region(1,inside) =0 ! inside is the dielectric region + shape_to_region(1,outside)=-1 ! no outside region + + end if ! infinite ground plane else ! there is an overshield specified so we do not need to offset the bundle to define a free space outer boundary @@ -418,7 +472,8 @@ IMPLICIT NONE rw(shape)=PUL%rw(i) rw2(shape)=PUL%rw2(i) rh(shape)=PUL%rh(i) - ro(shape)=PUL%o(i) + rox(shape)=PUL%ox(i) + roy(shape)=PUL%oy(i) rtheta(shape)=PUL%rtheta(i) if (shape_type(shape).EQ.circle) then dl(shape)=r(shape)/Laplace_surface_mesh_constant @@ -426,6 +481,8 @@ IMPLICIT NONE dl(shape)=min(rw(shape),rh(shape))/Laplace_surface_mesh_constant else if (shape_type(shape).EQ.Dshape) then dl(shape)=r(shape)/(2D0*Laplace_surface_mesh_constant) + else if (shape_type(shape).EQ.semicircle) then + dl(shape)=r(shape)/Laplace_surface_mesh_constant else write(*,*)'Unknown shape type',shape_type(shape),' i=',i,' shape=',shape end if @@ -434,22 +491,46 @@ IMPLICIT NONE if (ground_plane) then ! Add the ground plane information shape=PUL%n_conductors+1 ! ground plane - shape_type(shape)=rectangle - - rl=rboundary*Laplace_ground_plane_ratio ! width of the ground plane (x dimension) see constants.F90 for parameter - rt=rl*Laplace_ground_plane_thickness_ratio ! height of the ground plane (y dimension) see constants.F90 for parameter - rdl=rt/2d0 - x(shape)=xc ! x centre of ground plane rectangle - y(shape)=-rt ! y centre of ground plane rectangle - r(shape)=0d0 ! not used - rw(shape)=rl*2d0 ! x extent of the ground plane - rw2(shape)=rw(shape) ! x extent of the ground plane - rh(shape)=rt*2d0 ! y extent of the ground plane - ro(shape)=0d0 ! no offset from centre - rtheta(shape)=0d0 ! always zero for ground plane - dl(shape)=rdl ! this is the thickness of the ground plane - shape_to_region(shape,inside) =-1 ! no internal region - shape_to_region(shape,outside) = 0 ! no dielectric so a free space region + + if (finite_ground_plane) then + + shape_type(shape)=rectangle + rl=rboundary*Laplace_ground_plane_ratio ! width of the ground plane (x dimension) see constants.F90 for parameter + rt=rl*Laplace_ground_plane_thickness_ratio ! height of the ground plane (y dimension) see constants.F90 for parameter + rdl=rt/2d0 + x(shape)=xc ! x centre of ground plane rectangle + y(shape)=-rt ! y centre of ground plane rectangle + r(shape)=0d0 ! not used + rw(shape)=rl*2d0 ! x extent of the ground plane + rw2(shape)=rw(shape) ! x extent of the ground plane + rh(shape)=rt*2d0 ! y extent of the ground plane + rox(shape)=0d0 ! no offset from centre + roy(shape)=0d0 ! no offset from centre + rtheta(shape)=0d0 ! always zero for ground plane + dl(shape)=rdl ! this is the thickness of the ground plane + shape_to_region(shape,inside) =-1 ! no internal region + shape_to_region(shape,outside) = 0 ! no dielectric so a free space region + + else + + shape_type(shape)=semicircle_diameter + rl=0d0 ! not used + rt=0d0 ! not used + rdl=0d0 ! not used + x(shape)=0d0 ! x centre of ground plane + y(shape)=0d0 ! y centre of ground plane + r(shape)=rboundary ! semicircle radius + rw(shape)=0d0 ! not used + rw2(shape)=0d0 ! not used + rh(shape)=0d0 ! not used + rox(shape)=0d0 ! no offset from centre + roy(shape)=0d0 ! no offset from centre + rtheta(shape)=0d0 ! always zero for infinite ground plane + dl(shape)=r(1)/(2*Laplace_surface_mesh_constant) + shape_to_region(shape,inside) =-1 ! no internal region + shape_to_region(shape,outside) = 0 ! no dielectric so a free space region + + end if ! infinite ground plane else if (overshield) then ! Add the overshield conductor information @@ -461,7 +542,8 @@ IMPLICIT NONE rw(shape)=PUL%overshield_w rw2(shape)=PUL%overshield_w2 rh(shape)=PUL%overshield_h - ro(shape)=0d0 + rox(shape)=0d0 + roy(shape)=0d0 rtheta(shape)=0d0 dl(shape)=r(shape)/(2d0*Laplace_surface_mesh_constant) shape_to_region(shape,inside) = 0 ! inside the overshield is the free space region @@ -489,7 +571,8 @@ IMPLICIT NONE rw(shape)=PUL%rdw(i) rw2(shape)=PUL%rw2(i) rh(shape)=PUL%rdh(i) - ro(shape)=PUL%rdo(i) + rox(shape)=PUL%rdox(i) + roy(shape)=PUL%rdoy(i) rtheta(shape)=PUL%rtheta(i) if (shape_type(shape).EQ.circle) then @@ -525,8 +608,8 @@ IMPLICIT NONE ! calculate the centre of the conductor (the test point) ! xp,yp is the coordinate of the centre of the conductor relative to the centre of the cable - xp=ro(conductor_shape) - yp=0d0 + xp=rox(conductor_shape) + yp=roy(conductor_shape) ! rotate the cable about it's origin then translate to the cable position to give the final conductor position xt=x(conductor_shape)+xp*cos(rtheta(conductor_shape))-yp*sin(rtheta(conductor_shape)) @@ -537,7 +620,8 @@ IMPLICIT NONE dielectric_region=dielectric_region+1 ! check whether the conductor is inside the dielectric - if (point_is_inside(xt,yt,shape_type(shape),x(shape),y(shape),r(shape),rh(shape),rw(shape),ro(shape),rtheta(shape))) then + if (point_is_inside(xt,yt,shape_type(shape),x(shape),y(shape),r(shape),rh(shape),rw(shape), & + rox(shape),roy(shape),rtheta(shape))) then shape_to_region(conductor_shape,outside) =dielectric_region @@ -581,56 +665,121 @@ IMPLICIT NONE write(gmsh_geometry_file_unit,*)' // circle ',i point=point+1 - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+ro(i),',' ,y(i)-yc,',',0.0,',' ,dl(i),' };' + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' & + ,y(i)-yc+roy(i),',',0.0,',' ,dl(i),' };' point=point+1 - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+ro(i),',' ,y(i)-yc-r(i),',',0.0,',',dl(i),' };' + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' & + ,y(i)-yc-r(i)+roy(i),',',0.0,',',dl(i),' };' point=point+1 - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+ro(i)+r(i),',',y(i)-yc,',',0.0,',' ,dl(i),' };' + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)+r(i),',' & + ,y(i)-yc+roy(i),',',0.0,',' ,dl(i),' };' point=point+1 - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+ro(i),',' ,y(i)-yc+r(i),',',0.0,',',dl(i),' };' + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' & + ,y(i)-yc+r(i)+roy(i),',',0.0,',',dl(i),' };' point=point+1 - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+ro(i)-r(i),',',y(i)-yc,',',0.0,',' ,dl(i),' };' + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)-r(i),',' & + ,y(i)-yc+roy(i),',',0.0,',' ,dl(i),' };' else if (shape_type(i).EQ.rectangle) then write(gmsh_geometry_file_unit,*)' // rectangle ',i - point=point+1 - - xt=rw(i)/2d0+ro(i) - yt=rh(i)/2d0 + +! corner point + point=point+1 + xt=rw(i)/2d0+rox(i) + yt=rh(i)/2d0+roy(i) xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(i),' };' - point=point+1 - - xt=-rw(i)/2d0+ro(i) - yt=rh(i)/2d0 + edge_length=min(max_mesh_edge_length,dl(i)) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };' + +! edge centre point + point=point+1 + xt=rox(i) + yt=rh(i)/2d0+roy(i) xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(i),' };' - point=point+1 - - xt=-rw(i)/2d0+ro(i) - yt=-rh(i)/2d0 + dl_local=min(max_mesh_edge_length,rw(i)/Laplace_surface_mesh_constant) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };' + +! corner point + point=point+1 + xt=-rw(i)/2d0+rox(i) + yt=rh(i)/2d0+roy(i) xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(i),' };' - point=point+1 - - xt=rw(i)/2d0+ro(i) - yt=-rh(i)/2d0 + edge_length=min(max_mesh_edge_length,dl(i)) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };' + +! edge centre point + point=point+1 + xt=-rw(i)/2d0+rox(i) + yt=roy(i) + xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) + yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) + dl_local=min(max_mesh_edge_length,rh(i)/Laplace_surface_mesh_constant) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };' + +! corner point + point=point+1 + xt=-rw(i)/2d0+rox(i) + yt=-rh(i)/2d0+roy(i) + xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) + yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) + edge_length=min(max_mesh_edge_length,dl(i)) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };' + +! edge centre point + point=point+1 + xt=rox(i) + yt=-rh(i)/2d0+roy(i) + xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) + yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) + dl_local=min(max_mesh_edge_length,rw(i)/Laplace_surface_mesh_constant) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };' + +! corner point + point=point+1 + xt=rw(i)/2d0+rox(i) + yt=-rh(i)/2d0+roy(i) xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) - write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl(i),' };' + edge_length=min(max_mesh_edge_length,dl(i)) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',edge_length,' };' + +! edge centre point + point=point+1 + xt=rw(i)/2d0+rox(i) + yt=roy(i) + xp=x(i)-xc+xt*cos(rtheta(i))-yt*sin(rtheta(i)) + yp=y(i)-yc+xt*sin(rtheta(i))+yt*cos(rtheta(i)) + dl_local=min(max_mesh_edge_length,rh(i)/Laplace_surface_mesh_constant) + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl_local,' };' else if (shape_type(i).EQ.Dshape) then - CALL write_Dshape_gmsh(x(i),y(i),rw(i),rw2(i),rh(i),r(i),ro(i),rtheta(i),dl(i),point,i,gmsh_geometry_file_unit) + CALL write_Dshape_gmsh(x(i),y(i),rw(i),rw2(i),rh(i),r(i),rox(i),roy(i),rtheta(i),dl(i),point,i,gmsh_geometry_file_unit) + + else if (shape_type(i).EQ.semicircle) then - end if + write(gmsh_geometry_file_unit,*)' // semicircle ',i + point=point+1 + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i),',' & + ,y(i)-yc+roy(i),',',0.0,',' ,dl(i),' };' + point=point+1 + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)+r(i),',' & + ,y(i)-yc+roy(i),',',0.0,',',dl(i),' };' + point=point+1 + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i) ,',' & + ,y(i)-yc+r(i)+roy(i),',',0.0,',' ,dl(i),' };' + point=point+1 + write(gmsh_geometry_file_unit,*)'Point(',point,' ) = {',x(i)-xc+rox(i)-r(i),',' & + ,y(i)-yc+roy(i),',',0.0,',',dl(i),' };' + end if + end do ! next shape - + ! write the boudary line segment definitions for each shape, these are constructed from the previously defined points write(gmsh_geometry_file_unit,*)' ' @@ -666,6 +815,10 @@ IMPLICIT NONE p2=point+2 p3=point+3 p4=point+4 + p5=point+5 + p6=point+6 + p7=point+7 + p8=point+8 write(gmsh_geometry_file_unit,*)' // rectangle line segments ',i line=line+1 @@ -675,9 +828,17 @@ IMPLICIT NONE line=line+1 write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p3,',',p4,' };' line=line+1 - write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p4,',',p1,' };' + write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p4,',',p5,' };' + line=line+1 + write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p5,',',p6,' };' + line=line+1 + write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p6,',',p7,' };' + line=line+1 + write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p7,',',p8,' };' + line=line+1 + write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p8,',',p1,' };' - point=point+4 + point=point+8 else if (shape_type(i).EQ.Dshape) then @@ -713,7 +874,36 @@ IMPLICIT NONE write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p11,',',p12,',',p1,' };' point=point+12 - + + else if (shape_type(i).EQ.semicircle) then + + p1=point+1 + p2=point+2 + p3=point+3 + p4=point+4 + + write(gmsh_geometry_file_unit,*)' // semicircle ',i + line=line+1 + write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p2,',',p1,',',p3,' };' + line=line+1 + write(gmsh_geometry_file_unit,*)'Circle(',line,' ) = {',p3,',',p1,',',p4,' };' + + point=point+4 + + else if (shape_type(i).EQ.semicircle_diameter) then +! this is the infinite ground plane line so pick points from the outer boundary circle + p1=4 ! semicircle end point 2 + p2=1 ! centre + p3=2 ! semicircle end point 1 + + write(gmsh_geometry_file_unit,*)' // semicircle_diameter ',i + line=line+1 + write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p1,',',p2,' };' + gpline1=line + line=line+1 + write(gmsh_geometry_file_unit,*)'Line(',line,' ) = {',p2,',',p3,' };' + gpline2=line + end if end do ! next shape @@ -728,10 +918,11 @@ IMPLICIT NONE line_loop=0 ! initialise line_loop counter line=0 ! reset the line counter to the first line do i=1,n_shapes - - line_loop=line_loop+1 + + if ( (shape_type(i).EQ.circle) ) then - if ( (shape_type(i).EQ.circle).OR.(shape_type(i).EQ.rectangle) ) then + line_loop=line_loop+1 + loop_number(i)=line_loop write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',line+3,',',line+4,' };' @@ -739,11 +930,46 @@ IMPLICIT NONE do ii=1,4 write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number' end do - + line=line+4 + + else if ( shape_type(i).EQ.rectangle ) then + + line_loop=line_loop+1 + loop_number(i)=line_loop + + write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',line+3,',',line+4,',', & + line+5,',',line+6,',',line+7,',',line+8,' };' + +! write the boundary segment and the boundary number to the boundary file + do ii=1,8 + write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number' + end do + + line=line+8 + + else if (shape_type(i).EQ.semicircle) then + +! this is the outer boundary with infinite ground plane. +! the line loop consists of 4 lines, 2 for the semicircle and 2 for the ground plane + + line_loop=line_loop+1 + loop_number(i)=line_loop + + write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',gpline1,',',gpline2,' };' + +! write the boundary segment and the boundary number to the boundary file + do ii=1,2 + write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number' + end do + + line=line+2 else if (shape_type(i).EQ.Dshape) then + line_loop=line_loop+1 + loop_number(i)=line_loop + write(gmsh_geometry_file_unit,*)'Line Loop(',line_loop,' ) = {',line+1,',',line+2,',',line+3,',',line+4,',', & line+5,',',line+6,',',line+7,',',line+8,' };' @@ -753,6 +979,18 @@ IMPLICIT NONE end do line=line+8 + + + else if (shape_type(i).EQ.semicircle_diameter) then +! note there is no line loop for the semicircle_diameter shape as this is dealt with along with the semicircle +! just skip these lines + +! write the boundary segment and the boundary number to the boundary file + do ii=1,2 + write(boundary_file_unit,*)line+ii,i,' # boundary segment number and boundary number' + end do + + line=line+2 end if @@ -766,6 +1004,7 @@ IMPLICIT NONE surface=0 ! initialise surface counter ! loop over the dielectric regions (including free space) as each one contributes a surface to be meshed. + do dielectric_region=0,n_dielectric_regions ! 0 is the free space region which always exists surface=surface+1 @@ -785,32 +1024,39 @@ IMPLICIT NONE ! The overshield reference conductor should be the last conductor i=PUL%n_conductors + line_loop=loop_number(i) if( (shape_to_region(i,inside).EQ.dielectric_region).OR. & (shape_to_region(i,outside).EQ.dielectric_region) ) then - CALL add_integer_to_string(string1,i,string2) + CALL add_integer_to_string(string1,line_loop,string2) string1=trim(string2)//', ' end if end if ! overshield - - do i=1,n_shapes + + do i=1,n_shapes ! don't include ground plane in check + line_loop=loop_number(i) if ((.NOT.overshield).OR.(i.NE.PUL%n_conductors)) then if( (shape_to_region(i,inside).EQ.dielectric_region).OR. & - (shape_to_region(i,outside).EQ.dielectric_region) ) then - CALL add_integer_to_string(string1,i,string2) - string1=trim(string2)//', ' + (shape_to_region(i,outside).EQ.dielectric_region) ) then + if (line_loop.GT.0) then + CALL add_integer_to_string(string1,line_loop,string2) + string1=trim(string2)//', ' + end if end if end if - end do ! next shape + end do ! next shape - else + else ! not the fisrt surface ! dielectric boundaries are defined second so reverse the order of the loop for dielectric regions do i=n_shapes,1,-1 + line_loop=loop_number(i) if( (shape_to_region(i,inside).EQ.dielectric_region).OR. & - (shape_to_region(i,outside).EQ.dielectric_region) ) then - CALL add_integer_to_string(string1,i,string2) - string1=trim(string2)//', ' + (shape_to_region(i,outside).EQ.dielectric_region) ) then + if (line_loop.GT.0) then + CALL add_integer_to_string(string1,line_loop,string2) + string1=trim(string2)//', ' + end if end if end do @@ -834,13 +1080,15 @@ IMPLICIT NONE ! Dealllocate the local shape geometry data. if (allocated( shape_type )) DEALLOCATE( shape_type ) + if (allocated( loop_number )) DEALLOCATE( loop_number ) if (allocated( x )) DEALLOCATE( x ) if (allocated( y )) DEALLOCATE( y ) if (allocated( r )) DEALLOCATE( r ) if (allocated( rw )) DEALLOCATE( rw ) if (allocated( rw2 )) DEALLOCATE( rw2 ) if (allocated( rh )) DEALLOCATE( rh ) - if (allocated( ro )) DEALLOCATE( ro ) + if (allocated( rox )) DEALLOCATE( rox ) + if (allocated( roy )) DEALLOCATE( roy ) if (allocated( rtheta )) DEALLOCATE( rtheta ) if (allocated( dl )) DEALLOCATE( dl ) if (allocated( shape_to_region )) DEALLOCATE( shape_to_region ) diff --git a/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 b/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 index 10c088c..f1401c0 100644 --- a/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 +++ b/SRC/PUL_PARAMETER_CALCULATION/PUL_parameter_module.F90 @@ -40,6 +40,14 @@ !PUL_analytic.F90: SUBROUTINE PUL_LC_calc_overshield_wide_separation_approximation !PUL_analytic.F90: SUBROUTINE calculate_height_over_ground_plane !PUL_LC_Laplace.F90: SUBROUTINE PUL_LC_Laplace +!Aprod.F90: SUBROUTINE Aprod +!Aprod.F90: SUBROUTINE ATprod +!Aprod.F90: SUBROUTINE zAprod +!Aprod.F90: SUBROUTINE zATprod +!CG_solve.F90: SUBROUTINE solve_real_symm(n, b, x,tol) +!CG_solve.F90: SUBROUTINE solve_complex_symm(n, b, x,tol) +! +! ! ! NAME ! MODULE PUL_parameters @@ -47,13 +55,18 @@ ! Data relating to the calculation of PUL_parameters ! ! COMMENTS -! +! The conjugate gradient solutions are included in this module rather than the maths module +! as they are intimately linked with the matrix-vector product subroutine which in turn relies +! on the sparse matrix storage structure in this module ! ! HISTORY ! started 25/11/2015 CJS ! 30/3/2016 CJS fix factor of 2 in mutual inductance of wires over ground plane ! 18_4_2016 CJS include calculation for conductors within a cylindrical shield for oversheild domains ! 5_7_2016 CJS include numerical Laplace solver for L,C,G in inhomogeneous regions +! 13/3/2018CJS include iterative solver stuff based on Stanford University symmlq.f +! 16/3/2018 CJS add y offset for ML_flex_cable (PUL%ox, and PUL%oy)(PUL%rdox, and PUL%rdoy) +! 19/6/2018 CJS include iterative sparse matrix solver based on the biconjugate gradient method ! MODULE PUL_parameter_module @@ -74,14 +87,16 @@ TYPE::PUL_type real(dp),allocatable :: r(:) ! radius of a circular conductor real(dp),allocatable :: x(:) ! x coordinate of the centre of the cable to which this conductor belongs real(dp),allocatable :: y(:) ! y coordinate of the centre of the cable to which this conductor belongs - real(dp),allocatable :: o(:) ! offset in the x direction of this conductor from the cable centre (x():),y(:) above) + real(dp),allocatable :: ox(:) ! offset in the x direction of this conductor from the cable centre (x():),y(:) above) + real(dp),allocatable :: oy(:) ! offset in the y direction of this conductor from the cable centre (x():),y(:) above) real(dp),allocatable :: rh(:) ! height (y dimension) of rectangular conductor real(dp),allocatable :: rw(:) ! width1 (x dimension) of rectangular conductor/ Dshape real(dp),allocatable :: rw2(:) ! width2 (x dimension) of rectangular conductor/ Dshape real(dp),allocatable :: rd(:) ! radius of dielectric surrounding a circular conductor real(dp),allocatable :: rdh(:) ! height (y dimension) of rectangular dielectric around conductor real(dp),allocatable :: rdw(:) ! width (x dimension) of rectangular dielectric around conductor - real(dp),allocatable :: rdo(:) ! offset of dielectric in the x direction of this conductor from the cable centre + real(dp),allocatable :: rdox(:) ! offset of dielectric in the x direction of this conductor from the cable centre + real(dp),allocatable :: rdoy(:) ! offset of dielectric in the y direction of this conductor from the cable centre real(dp),allocatable :: rtheta(:) ! rotation angle of conductor type(Sfilter),allocatable :: epsr(:) ! relative permittivity of dielecrtric surrounding the conductor @@ -112,6 +127,15 @@ TYPE::PUL_type END TYPE PUL_type +integer,public :: n_entry + +! 1D arrays used in the construction of the K matrix ( K(i_K(:),j_K(:))=K(i_K(:),j_K(:))+s_K(:) ) +integer,allocatable,public :: i_K(:) +integer,allocatable,public :: j_K(:) +real(dp),allocatable,public :: s_K_re(:) +complex(dp),allocatable,public :: s_K(:) + + CONTAINS include "PUL_analytic.F90" @@ -120,6 +144,10 @@ include "PUL_LC_Laplace.F90" include "Laplace.F90" +include "Aprod.F90" + +include "CG_solve.F90" + ! NAME ! SUBROUTINE allocate_and_reset_PUL_data ! @@ -156,11 +184,13 @@ include "Laplace.F90" ALLOCATE( PUL%x( PUL%n_conductors) ) ALLOCATE( PUL%y( PUL%n_conductors) ) ALLOCATE( PUL%r( PUL%n_conductors) ) - ALLOCATE( PUL%o( PUL%n_conductors) ) + ALLOCATE( PUL%ox( PUL%n_conductors) ) + ALLOCATE( PUL%oy( PUL%n_conductors) ) ALLOCATE( PUL%rd( PUL%n_conductors) ) ALLOCATE( PUL%rdw( PUL%n_conductors) ) ALLOCATE( PUL%rdh( PUL%n_conductors) ) - ALLOCATE( PUL%rdo( PUL%n_conductors) ) + ALLOCATE( PUL%rdox( PUL%n_conductors) ) + ALLOCATE( PUL%rdoy( PUL%n_conductors) ) ALLOCATE( PUL%epsr( PUL%n_conductors) ) ALLOCATE( PUL%rh( PUL%n_conductors) ) ALLOCATE( PUL%rw( PUL%n_conductors) ) @@ -174,14 +204,16 @@ include "Laplace.F90" PUL%x(i)=0d0 PUL%y(i)=0d0 PUL%r(i)=0d0 - PUL%o(i)=0d0 + PUL%ox(i)=0d0 + PUL%oy(i)=0d0 PUL%rh(i)=0d0 PUL%rw(i)=0d0 PUL%rw2(i)=0d0 PUL%rd(i)=0d0 PUL%rdw(i)=0d0 PUL%rdh(i)=0d0 - PUL%rdo(i)=0d0 + PUL%rdox(i)=0d0 + PUL%rdoy(i)=0d0 PUL%rtheta(i)=0d0 PUL%epsr(i)=1d0 @@ -240,14 +272,16 @@ include "Laplace.F90" if (allocated(PUL%x ) ) DEALLOCATE( PUL%x ) if (allocated(PUL%y ) ) DEALLOCATE( PUL%y ) if (allocated(PUL%r ) ) DEALLOCATE( PUL%r ) - if (allocated(PUL%o ) ) DEALLOCATE( PUL%o ) + if (allocated(PUL%ox ) ) DEALLOCATE( PUL%ox ) + if (allocated(PUL%oy ) ) DEALLOCATE( PUL%oy ) if (allocated(PUL%rh ) ) DEALLOCATE( PUL%rh ) if (allocated(PUL%rw ) ) DEALLOCATE( PUL%rw ) if (allocated(PUL%rw2 ) ) DEALLOCATE( PUL%rw2 ) if (allocated(PUL%rd ) ) DEALLOCATE( PUL%rd ) if (allocated(PUL%rdw ) ) DEALLOCATE( PUL%rdw ) if (allocated(PUL%rdh ) ) DEALLOCATE( PUL%rdh ) - if (allocated(PUL%rdo ) ) DEALLOCATE( PUL%rdo ) + if (allocated(PUL%rdox ) ) DEALLOCATE( PUL%rdox ) + if (allocated(PUL%rdoy ) ) DEALLOCATE( PUL%rdoy ) if (allocated(PUL%rtheta ) ) DEALLOCATE( PUL%rtheta ) if (allocated(PUL%epsr ) ) then @@ -281,9 +315,10 @@ include "Laplace.F90" ! ! HISTORY ! started 6/10/2016 CJS +! add rox, roy instead of only the x offset, ro 16/3/2018 CJS ! - logical FUNCTION point_is_inside(xt,yt,shape_type,x,y,r,rh,rw,ro,rtheta) + logical FUNCTION point_is_inside(xt,yt,shape_type,x,y,r,rh,rw,rox,roy,rtheta) real(dp),intent(IN) :: xt ! test point x real(dp),intent(IN) :: yt ! test point y @@ -293,7 +328,8 @@ include "Laplace.F90" real(dp),intent(IN) :: r ! radius of cylindrical test object real(dp),intent(IN) :: rh ! height (y dimension) of rectangular test object real(dp),intent(IN) :: rw ! width (x dimension) of test object - real(dp),intent(IN) :: ro ! offset in the x direction of this test object from the centre + real(dp),intent(IN) :: rox ! offset in the x direction of this test object from the centre + real(dp),intent(IN) :: roy ! offset in the y direction of this test object from the centre real(dp),intent(IN) :: rtheta ! rotation angle of test object ! local variables @@ -335,8 +371,8 @@ include "Laplace.F90" write(*,*)'Inverse rotation on test point',xt_r,yt_r ! apply the inverse of the offset to the shape i.e. move it into the shifted, un-rotated coordinate system of the rectangle - xt_ro=xt_r-ro - yt_ro=yt_r + xt_ro=xt_r-rox + yt_ro=yt_r-roy write(*,*)'Offset Inverse rotation on test point',xt_ro,yt_ro @@ -345,9 +381,9 @@ include "Laplace.F90" dx=xt_ro dy=yt_ro - write(*,*)'dielectric centre point',x,y,' offset',ro,' angle',rtheta - write(*,*)'dx',dx,' rw/2',rw/2d0 - write(*,*)'dy',dy,' rh/2',rh/2d0 +! write(*,*)'dielectric centre point',x,y,' offset',rox,roy,' angle',rtheta +! write(*,*)'dx',dx,' rw/2',rw/2d0 +! write(*,*)'dy',dy,' rh/2',rh/2d0 if ( (abs(dx).LT.rw/2d0).AND.(abs(dy).LT.rh/2d0) ) then point_is_inside=.TRUE. @@ -442,7 +478,7 @@ include "Laplace.F90" ! started 15/11/2016 CJS ! - SUBROUTINE write_Dshape_gmsh(x,y,w_in,w2_in,h_in,r,ox,theta,dl,point,number,unit) + SUBROUTINE write_Dshape_gmsh(x,y,w_in,w2_in,h_in,r,ox,oy,theta,dl,point,number,unit) USE type_specifications USE general_module @@ -458,6 +494,7 @@ include "Laplace.F90" real(dp),intent(IN) :: h_in ! height (x dimension) of the Dshape real(dp),intent(IN) :: r ! radius of curves in Dshape real(dp),intent(IN) :: ox ! x offset + real(dp),intent(IN) :: oy ! y offset real(dp),intent(IN) :: theta ! rotation angle of Dshape real(dp),intent(IN) :: dl ! edge length for the line segments making up the Dshape @@ -498,7 +535,7 @@ include "Laplace.F90" ! POINT 1 ! top right xt=w1+ox - yt=h+r + yt=h+r+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -506,7 +543,7 @@ include "Laplace.F90" ! POINT 2 ! top left xt=-w1+ox - yt=h+r + yt=h+r+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -514,7 +551,7 @@ include "Laplace.F90" ! POINT 3 ! top left centre xt=-w1+ox - yt=h + yt=h+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -522,7 +559,7 @@ include "Laplace.F90" ! POINT 4 ! top left edge xt=-w1+ox+voxl - yt=h+voyl + yt=h+voyl+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -530,7 +567,7 @@ include "Laplace.F90" ! POINT 5 ! bottom left edge xt=-w2+ox+voxl - yt=-h+voyl + yt=-h+voyl+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -538,7 +575,7 @@ include "Laplace.F90" ! POINT 6 ! bottom left centre xt=-w2+ox - yt=-h + yt=-h+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -546,7 +583,7 @@ include "Laplace.F90" ! POINT 7 ! bottom left xt=-w2+ox - yt=-(h+r) + yt=-(h+r)+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -554,7 +591,7 @@ include "Laplace.F90" ! POINT 8 ! bottom right xt=w2+ox - yt=-(h+r) + yt=-(h+r)+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -562,7 +599,7 @@ include "Laplace.F90" ! POINT 9 ! bottom right centre xt=w2+ox - yt=-h + yt=-h+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -570,7 +607,7 @@ include "Laplace.F90" ! POINT 10 ! bottom right edge xt=w2+ox-voxl - yt=-h+voyl + yt=-h+voyl+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -578,7 +615,7 @@ include "Laplace.F90" ! POINT 11 ! top right edge xt=w1+ox-voxl - yt=h+voyl + yt=h+voyl+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' @@ -586,7 +623,7 @@ include "Laplace.F90" ! POINT 12 ! top right centre xt=w1+ox - yt=h + yt=h+oy xp=x+xt*cos(theta)-yt*sin(theta) yp=y+xt*sin(theta)+yt*cos(theta) write(unit,*)'Point(',point,' ) = {',xp,',',yp,',',0.0,',',dl,' };' diff --git a/SRC/cable_bundle_model_builder.F90 b/SRC/cable_bundle_model_builder.F90 index cc8f544..6a75b36 100644 --- a/SRC/cable_bundle_model_builder.F90 +++ b/SRC/cable_bundle_model_builder.F90 @@ -74,6 +74,8 @@ ! 29/6/2016 CJS: Check that the cables are all on one side of the ground plane if it is present ! December 2016 CJS Version 2: Rationalise cable types so that there is only a single version of each type of cable ! 24/2/2017 CJS Allow the input name to include a path i.e. the _spec file does not need to be local. +! 13/3/2018 CJS Add flag for direct/ iterative matrix solver in Laplace solution and inf/finite ground plane +! 19/6/2018 CJS Add flag for Neumann/ Asymptotic boundary condition in Laplace solver. Default is Neumann ! PROGRAM cable_bundle_model_builder @@ -333,6 +335,15 @@ integer :: dim if (INDEX(line,'no_plot_potential').NE.0) plot_potential=.FALSE. if (INDEX(line,'plot_mesh').NE.0) plot_mesh=.TRUE. if (INDEX(line,'no_plot_mesh').NE.0) plot_mesh=.FALSE. + + if (INDEX(line,'direct_solver').NE.0) direct_solver=.TRUE. + if (INDEX(line,'iterative_solver').NE.0) direct_solver=.FALSE. + + if (INDEX(line,'inf_gnd').NE.0) inf_gnd=.TRUE. + if (INDEX(line,'finite_gnd').NE.0) inf_gnd=.FALSE. + + if (INDEX(line,'abc').NE.0) use_ABC=.TRUE. + if (INDEX(line,'neumann').NE.0) use_ABC=.FALSE. ! redefine mesh generation parameters if required if (INDEX(line,'laplace_boundary_constant').NE.0) then @@ -346,6 +357,10 @@ integer :: dim if (INDEX(line,'twisted_pair_equivalent_radius').NE.0) then read(bundle_spec_file_unit,*,END=100,ERR=100)Twisted_pair_equivalent_radius end if + + if (INDEX(line,'max_mesh_edge_length').NE.0) then + read(bundle_spec_file_unit,*,END=100,ERR=100)max_mesh_edge_length + end if end do ! continue until all flags are read - indicated by an end of file. diff --git a/SRC/cable_model_builder.F90 b/SRC/cable_model_builder.F90 index a6336dc..ae06c14 100644 --- a/SRC/cable_model_builder.F90 +++ b/SRC/cable_model_builder.F90 @@ -84,6 +84,7 @@ ! 24/2/2017 CJS Allow the input name to include a path i.e. the _spec file does not need to be local. ! 23/10/2017 CJS Check the pole stability of the dielectric and transfer impedance filters. ! 16/11/2017 CJS Include network synthesis process to replace s-domain transfer functions +! 13/3/2018 CJS Add flag for direct/ iterative matrix solver in Laplace solution ! ! PROGRAM cable_model_builder @@ -214,10 +215,14 @@ integer :: i CALL spacewire_set_parameters(cable_spec) - else if (cable_spec%cable_type_string.eq.'flex_cable') then + else if (cable_spec%cable_type_string.eq.'original_flex_cable') then CALL flex_cable_set_parameters(cable_spec) + else if (cable_spec%cable_type_string.eq.'flex_cable') then + + CALL ML_flex_cable_set_parameters(cable_spec) + else if (cable_spec%cable_type_string.eq.'dconnector') then CALL Dconnector_set_parameters(cable_spec) @@ -274,11 +279,20 @@ integer :: i STOP 1 end if + if (cable_spec%n_parameters.NE.0) then + ! Check that the number of parameters is consistent with this cable type - if (np_in.NE.cable_spec%n_parameters) then - run_status='ERROR, Wrong number of parameters' - CALL write_program_status() - STOP 1 + if (np_in.NE.cable_spec%n_parameters) then + run_status='ERROR, Wrong number of parameters' + CALL write_program_status() + STOP 1 + end if + + else + +! for a ML_flex cable the number of parameters depends on the cable specification + cable_spec%n_parameters=np_in + end if ! allocate and read parameters @@ -410,6 +424,9 @@ integer :: i if (INDEX(line,'no_plot_potential').NE.0) plot_potential=.FALSE. if (INDEX(line,'plot_mesh').NE.0) plot_mesh=.TRUE. if (INDEX(line,'no_plot_mesh').NE.0) plot_mesh=.FALSE. + + if (INDEX(line,'direct_solver').NE.0) direct_solver=.TRUE. + if (INDEX(line,'iterative_solver').NE.0) direct_solver=.FALSE. ! redefine mesh generation parameters if required if (INDEX(line,'laplace_boundary_constant').NE.0) then @@ -423,6 +440,10 @@ integer :: i if (INDEX(line,'twisted_pair_equivalent_radius').NE.0) then read(cable_spec_file_unit,*,END=100,ERR=100)Twisted_pair_equivalent_radius end if + + if (INDEX(line,'max_mesh_edge_length').NE.0) then + read(cable_spec_file_unit,*,END=100,ERR=100)max_mesh_edge_length + end if end do ! continue until all flags are read - indicated by an end of file. @@ -482,6 +503,10 @@ integer :: i CALL flex_cable_set_internal_domain_information(cable_spec) + else if (cable_spec%cable_type.EQ.cable_geometry_type_ML_flex_cable) then + + CALL ML_flex_cable_set_internal_domain_information(cable_spec) + else if (cable_spec%cable_type.EQ.cable_geometry_type_dconnector) then CALL dconnector_set_internal_domain_information(cable_spec) diff --git a/TEST_CASES/FLEX_CABLE_3_AND_SINGLE_WIRE_EINC/flex_cable_3.cable_spec b/TEST_CASES/FLEX_CABLE_3_AND_SINGLE_WIRE_EINC/flex_cable_3.cable_spec index ed53828..6c5b699 100644 --- a/TEST_CASES/FLEX_CABLE_3_AND_SINGLE_WIRE_EINC/flex_cable_3.cable_spec +++ b/TEST_CASES/FLEX_CABLE_3_AND_SINGLE_WIRE_EINC/flex_cable_3.cable_spec @@ -2,13 +2,17 @@ . flex_cable 3 # number of conductors -6 # number of parameters -1.0e-3 # parameter 1: conductor width (x dimension) -0.25e-3 # parameter 2: conductor height (y dimension) -0.5e-3 # parameter 3: conductor separation (x dimension) -0.5e-3 # parameter 4: dielectric offset x -0.5e-3 # parameter 5: dielectric offset y -0.0 # parameter 6: conductivity +10 # number of parameters +5.0e-3 # parameter 1: dielectric width (x dimension) +1.25e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +3 # parameter 9: row 1 number of conductors +0.0 # parameter 10: conductivity 1 # number of frequency dependent parameters # dielectric relative permittivity model follows 1E9 # w normalisation constant diff --git a/TEST_CASES/FLEX_CABLE_3_CONDUCTOR/flex_cable_3.cable_spec b/TEST_CASES/FLEX_CABLE_3_CONDUCTOR/flex_cable_3.cable_spec index 6e77b95..8d82406 100644 --- a/TEST_CASES/FLEX_CABLE_3_CONDUCTOR/flex_cable_3.cable_spec +++ b/TEST_CASES/FLEX_CABLE_3_CONDUCTOR/flex_cable_3.cable_spec @@ -2,6 +2,33 @@ . flex_cable 3 # number of conductors +10 # number of parameters +4.4e-3 # parameter 1: dielectric width (x dimension) +0.45e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +3 # parameter 9: row 1 number of conductors +0.0 # parameter 10: conductivity +1 # number of frequency dependent parameters +# dielectric relative permittivity model follows + 1E9 # w normalisation constant + 0 # a order, a coefficients follow below: + 1.0 + 0 # b order, b coefficients follow below: + 1.0 + + + +# OLD model + +#MOD_cable_lib_dir +. +flex_cable +3 # number of conductors 6 # number of parameters 1.0e-3 # parameter 1: conductor width (x dimension) 0.25e-3 # parameter 2: conductor height (y dimension) @@ -16,11 +43,3 @@ flex_cable 1.0 0 # b order, b coefficients follow below: 1.0 - - -# dielectric relative permittivity model follows - 1E9 # w normalisation constant - 1 # a order, a coefficients follow below: - 4.0 2.0 - 1 # b order, b coefficients follow below: - 1.0 1.0 diff --git a/TEST_CASES/FLEX_CABLE_3_EINC/flex_cable_3.cable_spec b/TEST_CASES/FLEX_CABLE_3_EINC/flex_cable_3.cable_spec index ed53828..63bb0ba 100644 --- a/TEST_CASES/FLEX_CABLE_3_EINC/flex_cable_3.cable_spec +++ b/TEST_CASES/FLEX_CABLE_3_EINC/flex_cable_3.cable_spec @@ -2,6 +2,33 @@ . flex_cable 3 # number of conductors +10 # number of parameters +5.0e-3 # parameter 1: dielectric width (x dimension) +1.25e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +3 # parameter 9: row 1 number of conductors +0.0 # parameter 10: conductivity +1 # number of frequency dependent parameters +# dielectric relative permittivity model follows + 1E9 # w normalisation constant + 0 # a order, a coefficients follow below: + 2.2 + 0 # b order, b coefficients follow below: + 1.0 + + + + + +#MOD_cable_lib_dir +. +flex_cable +3 # number of conductors 6 # number of parameters 1.0e-3 # parameter 1: conductor width (x dimension) 0.25e-3 # parameter 2: conductor height (y dimension) diff --git a/TEST_CASES/FLEX_CABLE_3_GROUND_PLANE/flex_cable_3.cable_spec b/TEST_CASES/FLEX_CABLE_3_GROUND_PLANE/flex_cable_3.cable_spec index ed53828..c34411e 100644 --- a/TEST_CASES/FLEX_CABLE_3_GROUND_PLANE/flex_cable_3.cable_spec +++ b/TEST_CASES/FLEX_CABLE_3_GROUND_PLANE/flex_cable_3.cable_spec @@ -2,6 +2,32 @@ . flex_cable 3 # number of conductors +10 # number of parameters +5.0e-3 # parameter 1: dielectric width (x dimension) +1.25e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +3 # parameter 9: row 1 number of conductors +0.0 # parameter 10: conductivity +1 # number of frequency dependent parameters +# dielectric relative permittivity model follows + 1E9 # w normalisation constant + 0 # a order, a coefficients follow below: + 2.2 + 0 # b order, b coefficients follow below: + 1.0 + + + + +#MOD_cable_lib_dir +. +flex_cable +3 # number of conductors 6 # number of parameters 1.0e-3 # parameter 1: conductor width (x dimension) 0.25e-3 # parameter 2: conductor height (y dimension) diff --git a/TEST_CASES/FLEX_CABLE_3_OVERSHIELD/flex_cable_3.cable_spec b/TEST_CASES/FLEX_CABLE_3_OVERSHIELD/flex_cable_3.cable_spec index ed53828..c34411e 100644 --- a/TEST_CASES/FLEX_CABLE_3_OVERSHIELD/flex_cable_3.cable_spec +++ b/TEST_CASES/FLEX_CABLE_3_OVERSHIELD/flex_cable_3.cable_spec @@ -2,6 +2,32 @@ . flex_cable 3 # number of conductors +10 # number of parameters +5.0e-3 # parameter 1: dielectric width (x dimension) +1.25e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +3 # parameter 9: row 1 number of conductors +0.0 # parameter 10: conductivity +1 # number of frequency dependent parameters +# dielectric relative permittivity model follows + 1E9 # w normalisation constant + 0 # a order, a coefficients follow below: + 2.2 + 0 # b order, b coefficients follow below: + 1.0 + + + + +#MOD_cable_lib_dir +. +flex_cable +3 # number of conductors 6 # number of parameters 1.0e-3 # parameter 1: conductor width (x dimension) 0.25e-3 # parameter 2: conductor height (y dimension) diff --git a/TEST_CASES/FLEX_CABLE_DIELECTRIC_3_CONDUCTOR/flex_cable_3.cable_spec b/TEST_CASES/FLEX_CABLE_DIELECTRIC_3_CONDUCTOR/flex_cable_3.cable_spec index ed53828..b6390d0 100644 --- a/TEST_CASES/FLEX_CABLE_DIELECTRIC_3_CONDUCTOR/flex_cable_3.cable_spec +++ b/TEST_CASES/FLEX_CABLE_DIELECTRIC_3_CONDUCTOR/flex_cable_3.cable_spec @@ -2,6 +2,32 @@ . flex_cable 3 # number of conductors +10 # number of parameters +5.0e-3 # parameter 1: dielectric width (x dimension) +1.25e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +3 # parameter 9: row 1 number of conductors +0.0 # parameter 10: conductivity +1 # number of frequency dependent parameters +# dielectric relative permittivity model follows + 1E9 # w normalisation constant + 0 # a order, a coefficients follow below: + 2.2 + 0 # b order, b coefficients follow below: + 1.0 + + + + +#MOD_cable_lib_dir +. +flex_cable +3 # number of conductors 6 # number of parameters 1.0e-3 # parameter 1: conductor width (x dimension) 0.25e-3 # parameter 2: conductor height (y dimension) @@ -16,3 +42,29 @@ flex_cable 2.2 0 # b order, b coefficients follow below: 1.0 +#MOD_cable_lib_dir +. +ML_flex_cable +3 # number of conductors +10 # number of parameters +5.0e-3 # parameter 1: dielectric width (x dimension) +1.25e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +3 # parameter 9: row 1 number of conductors +0.0 # parameter 10: conductivity +1 # number of frequency dependent parameters +# dielectric relative permittivity model follows + 1E9 # w normalisation constant + 0 # a order, a coefficients follow below: + 2.2 + 0 # b order, b coefficients follow below: + 1.0 + + + + diff --git a/TEST_CASES/FLEX_CABLE_LOSSY_2_CONDUCTOR/flex_cable_2.cable_spec b/TEST_CASES/FLEX_CABLE_LOSSY_2_CONDUCTOR/flex_cable_2.cable_spec index 8b9f3a7..91e80cb 100644 --- a/TEST_CASES/FLEX_CABLE_LOSSY_2_CONDUCTOR/flex_cable_2.cable_spec +++ b/TEST_CASES/FLEX_CABLE_LOSSY_2_CONDUCTOR/flex_cable_2.cable_spec @@ -2,6 +2,32 @@ . flex_cable 2 # number of conductors +10 # number of parameters +2.9e-3 # parameter 1: dielectric width (x dimension) +0.45e-3 # parameter 2: dielectric height (y dimension) +1 # parameter 3: number of rows of conductors +0.0e-3 # parameter 4: row 1 centre offset x +0.0e-3 # parameter 5: row 1 centre offset y +1.0e-3 # parameter 6: row 1 conductor width (x dimension) +0.25e-3 # parameter 7: row 1 conductor height (y dimension) +0.5e-3 # parameter 8: row 1 conductor separation +2 # parameter 9: row 1 number of conductors +5e7 # parameter 10: conductivity +1 # number of frequency dependent parameters +# dielectric relative permittivity model follows + 1E9 # w normalisation constant + 0 # a order, a coefficients follow below: + 2.25 + 0 # b order, b coefficients follow below: + 1.0 + + + + +#MOD_cable_lib_dir +. +flex_cable +2 # number of conductors 6 # number of parameters 1.0e-3 # parameter 1: conductor width (x dimension) 0.25e-3 # parameter 2: conductor height (y dimension) diff --git a/TEST_CASES/generate_spice_cable_bundle_model b/TEST_CASES/generate_spice_cable_bundle_model index 1444172..dddea4a 100755 --- a/TEST_CASES/generate_spice_cable_bundle_model +++ b/TEST_CASES/generate_spice_cable_bundle_model @@ -163,7 +163,6 @@ CM_TO_DM_TWISTED_PAIR_OVER_GROUND_PLANE LOW_MASS_SPACEWIRE " - ERROR_STRING="\ Run using the following command: diff --git a/clean_project b/clean_project new file mode 100755 index 0000000..4a11dfa --- /dev/null +++ b/clean_project @@ -0,0 +1,16 @@ +# Script to clean object files, mod files and +# all temporary test case files before updating +# the git repository. + +cd SRC +make clean +rm -f compilation_date.inc +cd .. +# +cd TEST_CASES +./generate_spice_cable_bundle_model clean_all +cd .. +# +cd DOCUMENTATION +make clean +cd .. -- libgit2 0.21.2