diff --git a/.github/workflows/ubuntu-gnu.yml b/.github/workflows/ubuntu-gnu.yml index aa7d5d57..7f36f152 100644 --- a/.github/workflows/ubuntu-gnu.yml +++ b/.github/workflows/ubuntu-gnu.yml @@ -1,45 +1,45 @@ -name: ubuntu-gnu - -permissions: - actions: write - -on: - push: - branches: - - main - - dev - pull_request: - branches: - - main - - dev - -jobs: - builds-and-test: - runs-on: ubuntu-latest - steps: - - name: Checkout - uses: actions/checkout@v4 - with: - submodules: 'recursive' - - - name: Setup ninja - uses: seanmiddleditch/gha-setup-ninja@master - - - name: Build application - run: | - cmake -S . -B build -G Ninja -DCMAKE_BUILD_TYPE=Release - cmake --build build -j - - - name: Run fdtd unit tests - run: | - build/bin/fdtd_tests - - - name: Install python wrapper requirements - run: | - python -m pip install -r requirements.txt - - - name: Run wrapper tests - run: | - python -m pytest test - - \ No newline at end of file +name: ubuntu-gnu + +permissions: + actions: write + +on: + push: + branches: + - main + - dev + pull_request: + branches: + - main + - dev + +jobs: + builds-and-test: + runs-on: ubuntu-latest + steps: + - name: Checkout + uses: actions/checkout@v4 + with: + submodules: 'recursive' + + - name: Setup ninja + uses: seanmiddleditch/gha-setup-ninja@master + + - name: Build application + run: | + cmake -S . -B build -G Ninja -DCMAKE_BUILD_TYPE=Release + cmake --build build -j + + - name: Run fdtd unit tests + run: | + build/bin/fdtd_tests + + - name: Install python wrapper requirements + run: | + python -m pip install -r requirements.txt + + - name: Run wrapper tests + run: | + python -m pytest test + + diff --git a/.github/workflows/ubuntu-intelLLVM.yml b/.github/workflows/ubuntu-intelLLVM.yml index 38738da1..407ec60e 100644 --- a/.github/workflows/ubuntu-intelLLVM.yml +++ b/.github/workflows/ubuntu-intelLLVM.yml @@ -81,4 +81,4 @@ jobs: - name: exclude unused files from cache if: steps.cache-install.outputs.cache-hit != 'true' - run: sh -c .github/workflows/oneapi_cache_exclude_linux.sh \ No newline at end of file + run: sh -c .github/workflows/oneapi_cache_exclude_linux.sh diff --git a/.github/workflows/windows-intelLLVM.yml b/.github/workflows/windows-intelLLVM.yml index 1e70dddc..918b5f92 100644 --- a/.github/workflows/windows-intelLLVM.yml +++ b/.github/workflows/windows-intelLLVM.yml @@ -54,7 +54,7 @@ jobs: shell: bash run: | .github/workflows/build_windows.bat - + - name: Install python wrapper requirements run: | python -m pip install -r requirements.txt @@ -63,9 +63,9 @@ jobs: shell: bash run: | .github/workflows/run_tests_windows.bat - + - name: exclude unused files from cache if: steps.cache-install.outputs.cache-hit != 'true' shell: bash - run: .github/workflows/oneapi_cache_exclude_windows.sh \ No newline at end of file + run: .github/workflows/oneapi_cache_exclude_windows.sh diff --git a/CMakeLists.txt b/CMakeLists.txt index 767bb183..36bda8be 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required (VERSION 3.18) -ENABLE_LANGUAGE(Fortran) +enable_language (CXX Fortran) project(semba-fdtd Fortran) set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) @@ -49,6 +49,7 @@ if (CompileExecutable) "src_main_pub/magneticdispersive.F90" "src_main_pub/maloney_nostoch.F90" "src_main_pub/mpicomm.F90" + "src_main_pub/nfde_rotate.F90" "src_main_pub/nodalsources.F90" "src_main_pub/observation.F90" "src_main_pub/planewaves.F90" @@ -63,19 +64,20 @@ if (CompileExecutable) "src_main_pub/version.F90" "src_wires_pub/wires_types.F90" "src_wires_pub/wires.F90" + "src_wires_pub/wires_mtln.F90" ) - target_link_libraries(semba-fdtd semba-types smbjson) + target_link_libraries(semba-fdtd semba-types smbjson mtlnsolver ngspice_interface ${NGSPICE_LIB}) include_directories(${CMAKE_BINARY_DIR}/mod) endif() if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") message(STATUS "Gfortran found...") - set(CMAKE_Fortran_FLAGS "-ffree-form -ffree-line-length-none -fdec -fcheck=bounds") + set(CMAKE_Fortran_FLAGS "-ffree-form -ffree-line-length-none -fdec -fcheck=bounds -fopenmp") add_definitions(-DCompileWithIncludeMpifh) endif() if(CMAKE_Fortran_COMPILER_ID MATCHES "^Intel") message(STATUS "INTEL FORTRAN found...") - set(CMAKE_Fortran_FLAGS "-fpp -check all,nouninit -debug full -traceback ") + set(CMAKE_Fortran_FLAGS "-fpp -check all,nouninit -debug full -traceback") endif() if (CompileWithHDF) diff --git a/README.md b/README.md index 830b7121..099b98d1 100755 --- a/README.md +++ b/README.md @@ -49,7 +49,9 @@ IEEE Transactions on Electromagnetic Compatibility. 2016. https://doi.org/10.110 # Usage -## Compilation +## Compilation and testing + +Tests must be run from the root folder. `python` wrapper test assumes that `semba-fdtd` has been compiled successfully and is located in folder `build/bin/`. For intel compilation it also assumes that the intel runtime libraries are accessible. ## Running cases diff --git a/doc/smbjson.md b/doc/smbjson.md index e1e6c1b3..1ad5a9d9 100644 --- a/doc/smbjson.md +++ b/doc/smbjson.md @@ -592,12 +592,16 @@ If not `magnitudeFile` is specified and only one `source` is defined, the `magni #### `movie` Probes of type `movie` record a vector field in a volume region indicated by `elementIds`. `[field]` can be `electric`, `magnetic`, or `currentDensity`; defaults to `electric`. +`currentDensity` will store only the surface density currents on `pec` or lossy surfaces. +The stored values can be selected using `[components]`, which stores an array of the following labels `x`, `y`, `z`, or `magnitude`; if no components are specified, defaults to `magnitude`. +An example follows: ```json { "name": "electric_field_movie", "type": "movie", "field": "electric", + "components": ["x"], "elementIds": [4] } ``` diff --git a/requirements.txt b/requirements.txt index a6c4edea..39affa7e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ pytest numpy matplotlib +pandas \ No newline at end of file diff --git a/src_json_parser/labels_mod.F90 b/src_json_parser/labels_mod.F90 index c4466660..90ab3326 100644 --- a/src_json_parser/labels_mod.F90 +++ b/src_json_parser/labels_mod.F90 @@ -9,6 +9,7 @@ module labels_mod character (len=*), parameter :: J_DIR_X = "x" character (len=*), parameter :: J_DIR_Y = "y" character (len=*), parameter :: J_DIR_Z = "z" + character (len=*), parameter :: J_DIR_M = 'magnitude' character (len=*), parameter :: J_FIELD = "field" character (len=*), parameter :: J_FIELD_ELECTRIC = "electric" @@ -139,6 +140,8 @@ module labels_mod character (len=*), parameter :: J_PR_POINT_DIRECTIONS = "directions" + character (len=*), parameter :: J_PR_MOVIE_COMPONENTS = "components" + character (len=*), parameter :: J_PR_FAR_FIELD_THETA = "theta" character (len=*), parameter :: J_PR_FAR_FIELD_PHI = "phi" character (len=*), parameter :: J_PR_FAR_FIELD_DIR_INITIAL = "initial" diff --git a/src_json_parser/nfdetypes_extension.F90 b/src_json_parser/nfdetypes_extension.F90 index 5dabedf3..6c41d713 100644 --- a/src_json_parser/nfdetypes_extension.F90 +++ b/src_json_parser/nfdetypes_extension.F90 @@ -128,7 +128,9 @@ subroutine initializeProblemDescription(pD) allocate(pD%BloquePrb) allocate(pD%BloquePrb%bp(0)) + allocate(pD%VolPrb) + allocate(pD%VolPrb%collection(0)) ! allocate(pD%tWires) allocate(pD%tWires%tw(0)) @@ -766,28 +768,28 @@ elemental logical function abstractSonda_eq(a, b) result(res) res = .false. if (a%n_FarField /= b%n_FarField) return - if (a%n_Electric /= b%n_Electric) return - if (a%n_Magnetic /= b%n_Magnetic) return - if (a%n_NormalElectric /= b%n_NormalElectric) return - if (a%n_NormalMagnetic /= b%n_NormalMagnetic) return - if (a%n_SurfaceElectricCurrent /= b%n_SurfaceElectricCurrent) return - if (a%n_SurfaceMagneticCurrent /= b%n_SurfaceMagneticCurrent) return + ! if (a%n_Electric /= b%n_Electric) return + ! if (a%n_Magnetic /= b%n_Magnetic) return + ! if (a%n_NormalElectric /= b%n_NormalElectric) return + ! if (a%n_NormalMagnetic /= b%n_NormalMagnetic) return + ! if (a%n_SurfaceElectricCurrent /= b%n_SurfaceElectricCurrent) return + ! if (a%n_SurfaceMagneticCurrent /= b%n_SurfaceMagneticCurrent) return if (.not. associated(a%FarField) .or. .not. associated(b%FarField)) return - if (.not. associated(a%Electric) .or. .not. associated(b%Electric)) return - if (.not. associated(a%Magnetic) .or. .not. associated(b%Magnetic)) return - if (.not. associated(a%NormalElectric) .or. .not. associated(b%NormalElectric)) return - if (.not. associated(a%NormalMagnetic) .or. .not. associated(b%NormalMagnetic)) return - if (.not. associated(a%SurfaceElectricCurrent) .or. .not. associated(b%SurfaceElectricCurrent)) return - if (.not. associated(a%SurfaceMagneticCurrent) .or. .not. associated(b%SurfaceMagneticCurrent)) return + ! if (.not. associated(a%Electric) .or. .not. associated(b%Electric)) return + ! if (.not. associated(a%Magnetic) .or. .not. associated(b%Magnetic)) return + ! if (.not. associated(a%NormalElectric) .or. .not. associated(b%NormalElectric)) return + ! if (.not. associated(a%NormalMagnetic) .or. .not. associated(b%NormalMagnetic)) return + ! if (.not. associated(a%SurfaceElectricCurrent) .or. .not. associated(b%SurfaceElectricCurrent)) return + ! if (.not. associated(a%SurfaceMagneticCurrent) .or. .not. associated(b%SurfaceMagneticCurrent)) return if (any(.not. a%FarField == b%FarField)) return - if (any(.not. a%Electric == b%Electric)) return - if (any(.not. a%Magnetic == b%Magnetic)) return - if (any(.not. a%NormalElectric == b%NormalElectric)) return - if (any(.not. a%NormalMagnetic == b%NormalMagnetic)) return - if (any(.not. a%SurfaceElectricCurrent == b%SurfaceElectricCurrent)) return - if (any(.not. a%SurfaceMagneticCurrent == b%SurfaceMagneticCurrent)) return + ! if (any(.not. a%Electric == b%Electric)) return + ! if (any(.not. a%Magnetic == b%Magnetic)) return + ! if (any(.not. a%NormalElectric == b%NormalElectric)) return + ! if (any(.not. a%NormalMagnetic == b%NormalMagnetic)) return + ! if (any(.not. a%SurfaceElectricCurrent == b%SurfaceElectricCurrent)) return + ! if (any(.not. a%SurfaceMagneticCurrent == b%SurfaceMagneticCurrent)) return res = .true. end function abstractSonda_eq @@ -915,8 +917,8 @@ end function bloqueprobes_eq elemental logical function volprobe_eq(a, b) type(VolProbe), intent(in) :: a, b - volprobe_eq = & - all(a%cordinates == b%cordinates(1:b%len_cor)) .and. & + volprobe_eq = all(a%cordinates == b%cordinates) + volprobe_eq = volprobe_eq .and. & (a%tstart == b%tstart) .and. & (a%tstop == b%tstop) .and. & (a%tstep == b%tstep) .and. & @@ -931,6 +933,10 @@ end function volprobe_eq elemental logical function volprobes_eq(a, b) type(VolProbes), intent(in) :: a, b + if (.not. associated(a%collection) .or. .not. associated(b%collection)) then + volprobes_eq = .false. + return + end if volprobes_eq = & (a%length == b%length) .and. & (a%length_max == b%length_max) .and. & diff --git a/src_json_parser/smbjson.F90 b/src_json_parser/smbjson.F90 index 0c2628f2..ce61ab8b 100644 --- a/src_json_parser/smbjson.F90 +++ b/src_json_parser/smbjson.F90 @@ -16,7 +16,6 @@ module smbjson implicit none - character (len=*), parameter :: SMBJSON_LOG_SUFFIX = "_log_" integer, private, parameter :: MAX_LINE = 256 type, public :: parser_t @@ -44,9 +43,8 @@ module smbjson procedure :: readProbes procedure :: readMoreProbes procedure :: readBlockProbes + procedure :: readVolumicProbes procedure :: readThinWires - procedure :: readSlantedWires - procedure :: readThinSlots ! procedure :: readMesh ! @@ -97,8 +95,8 @@ function parser_ctor(filename) result(res) subroutine initializeJson(this) class(parser_t) :: this - integer :: stat - + integer :: stat + allocate(this%jsonfile) call this%jsonfile%initialize() if (this%jsonfile%failed()) then @@ -120,8 +118,8 @@ subroutine initializeJson(this) function readProblemDescription(this) result (res) class(parser_t) :: this type(Parseador) :: res - integer :: stat - + integer :: stat + allocate(this%jsonfile) call this%jsonfile%initialize() if (this%jsonfile%failed()) then @@ -142,8 +140,6 @@ function readProblemDescription(this) result (res) this%mesh = this%readMesh() this%matTable = IdChildTable_t(this%core, this%root, J_MATERIALS) - - call initializeProblemDescription(res) ! Basics @@ -155,25 +151,19 @@ function readProblemDescription(this) result (res) ! Materials res%pecRegs = this%readPECRegions() res%pmcRegs = this%readPMCRegions() - ! res%DielRegs = this%readDielectricRegions() - ! res%LossyThinSurfs = this%readLossyThinSurfaces() - ! res%frqDepMats = this%readFrequencyDependentMaterials() - ! res%aniMats = this%readAnisotropicMaterials() ! Sources - ! res%boxSrc = this%readBoxSources() res%plnSrc = this%readPlanewaves() res%nodSrc = this%readNodalSources() + ! Probes res%oldSonda = this%readProbes() res%sonda = this%readMoreProbes() res%BloquePrb = this%readBlockProbes() - ! res%VolPrb = this%readVolumicProbes() + res%VolPrb = this%readVolumicProbes() + ! Thin elements res%tWires = this%readThinWires() - res%sWires = this%readSlantedWires() - res%tSlots = this%readThinSlots() - ! mtln res%mtln = this%readMTLN(res%despl) @@ -183,7 +173,6 @@ function readProblemDescription(this) result (res) call this%jsonfile%destroy() nullify(this%root) - end function function readMesh(this) result(res) @@ -195,7 +184,7 @@ function readMesh(this) result(res) type(coordinate_t) :: c integer :: stat logical :: found - + call this%core%get(this%root, J_MESH//'.'//J_COORDINATES, jcs, found=found) if (found) then call res%allocateCoordinates(10*this%core%count(jcs)) @@ -411,23 +400,6 @@ function buildPECPMCRegion(cRs) result(res) res%nVols_max = size(res%Vols) end function - function readDielectricRegions() result (res) - type(DielectricRegions) :: res - ! TODO - end function - - function readLossyThinSurfaces() result (res) - type(LossyThinSurfaces) :: res - ! TODO - end function - - function readFrequencyDependentMaterials() result (res) - type(FreqDepenMaterials) :: res - - - ! TODO - end function - function readPlanewaves(this) result (res) class(parser_t) :: this type(PlaneWaves) :: res @@ -596,6 +568,7 @@ function readFarFieldProbe(p) result (res) allocate(res%FarField(1)) ff => res%FarField(1)%probe + ff%grname = " " call this%core%get(p, J_NAME, outputName) ff%outputrequest = trim(adjustl(outputName)) @@ -603,20 +576,22 @@ function readFarFieldProbe(p) result (res) domain = this%getDomain(p, J_PR_DOMAIN) if (domain%type2 /= NP_T2_FREQ) & write(error_unit, *) "ERROR at far field probe: Only accepted domain is frequency." - + ff%tstart = 0.0 + ff%tstop = 0.0 + ff%tstep = 0.0 ff%fstart = domain%fstart ff%fstop = domain%fstop ff%fstep = domain%fstep - + block logical :: sourcesFound type(json_value), pointer :: sources, src character (len=:), allocatable :: fn - + fn = this%getStrAt(p, J_PR_DOMAIN//J_PR_DOMAIN_MAGNITUDE_FILE, found=transferFunctionFound) if (.not. transferFunctionFound) then call this%core%get(this%root, J_SOURCES, sources, sourcesFound) - if (sourcesFound) then + if (sourcesFound) then if (this%core%count(sources) == 1) then call this%core%get_child(sources, 1, src) call this%core%get(src, J_SRC_MAGNITUDE_FILE, fn, found=transferFunctionFound) @@ -624,16 +599,16 @@ function readFarFieldProbe(p) result (res) end if end if - if (transferFunctionFound) then + if (transferFunctionFound) then ff%FileNormalize = trim(adjustl(fn)) else ff%FileNormalize = " " end if - + end block - + if (domain%isLogarithmicFrequencySpacing) then - ff%outputrequest = ff%outputrequest // SMBJSON_LOG_SUFFIX + call appendLogSufix(ff%outputrequest) end if block @@ -645,6 +620,7 @@ function readFarFieldProbe(p) result (res) allocate(ff%i(2)) allocate(ff%j(2)) allocate(ff%k(2)) + allocate(ff%node(0)) ff%i(1) = nfdeCoords(1)%Xi ff%i(2) = nfdeCoords(1)%Xe ff%j(1) = nfdeCoords(1)%Yi @@ -770,7 +746,7 @@ subroutine setDomain(res, domain) res%type2 = domain%type2 if (domain%isLogarithmicFrequencySpacing) then - res%outputrequest = res%outputrequest // SMBJSON_LOG_SUFFIX + call appendLogSufix(res%outputrequest) end if end subroutine @@ -845,9 +821,9 @@ function readBlockProbe(bp) result(res) type(cell_region_t), dimension(:), allocatable :: cRs cRs = this%mesh%getCellRegions(this%getIntsAt(bp, J_ELEMENTIDS)) - if (size(cRs) /= 1) write(error_unit, *) "Block probe must be defined by a single cell region." + if (size(cRs) /= 1) write(error_unit, *) "Bulk current probe must be defined by a single cell region." - if (size(cRs(1)%intervals) /= 1) write(error_unit, *) "Block probe must be defined by a single cell interval." + if (size(cRs(1)%intervals) /= 1) write(error_unit, *) "Bulk current probe must be defined by a single cell interval." cs = cellIntervalsToCoords(cRs(1)%intervals) res%i1 = cs(1)%xi @@ -862,7 +838,7 @@ function readBlockProbe(bp) result(res) call setDomain(res, this%getDomain(bp, J_PR_DOMAIN)) res%skip = 1 - res%tag = trim(adjustl(this%getStrAt(bp, J_NAME))) + res%tag = trim(adjustl(this%getStrAt(bp, J_NAME, default=" "))) res%t = BcELECT end function @@ -881,16 +857,157 @@ subroutine setDomain(res, domain) res%type2 = domain%type2 if (domain%isLogarithmicFrequencySpacing) then - res%outputrequest = res%outputrequest // SMBJSON_LOG_SUFFIX + call appendLogSufix(res%outputrequest) end if end subroutine end function - function readVolumicProbes() result (res) + function readVolumicProbes(this) result (res) + class(parser_t) :: this type(VolProbes) :: res - ! TODO + type(json_value_ptr), dimension(:), allocatable :: ps + type(json_value), pointer :: probes + logical :: found + integer :: i + + call this%core%get(this%root, J_PROBES, probes, found) + if (.not. found) then + res = buildNoVolProbes() + return + end if + + ps = this%jsonValueFilterByKeyValues(probes, J_TYPE, [J_PR_TYPE_MOVIE]) + if (size(ps) == 0) then + res = buildNoVolProbes() + return + end if + + res%length = size(ps) + res%length_max = size(ps) + res%len_cor_max = 2*size(ps) + allocate(res%collection(size(ps))) + do i = 1, size(ps) + res%collection(i) = readVolProbe(ps(i)%p) + end do + + contains + function buildNoVolProbes() result(res) + type(VolProbes) :: res + allocate(res%collection(0)) + res%length = 0 + res%length_max = 0 + res%len_cor_max = 0 + end function + + function readVolProbe(p) result(res) + type(VolProbe) :: res + type(json_value), pointer :: p, compsPtr, compPtr + type(coords), dimension(:), allocatable :: cs + type(cell_region_t), dimension(:), allocatable :: cRs + character(len=:), allocatable :: fieldType, component + integer :: i + integer :: numberOfComponents + logical :: componentsFound + + cRs = this%mesh%getCellRegions(this%getIntsAt(p, J_ELEMENTIDS)) + if (size(cRs) /= 1) then + write(error_unit, *) "Movie probe must be defined over a single cell region." + end if + + if (size(cRs(1)%intervals) /= 1) then + write(error_unit, *) "Movie probe must be defined by a single cell interval." + end if + cs = cellIntervalsToCoords(cRs(1)%intervals) + + fieldType = this%getStrAt(p, J_FIELD, default=J_FIELD_ELECTRIC) + call this%core%get(p, J_PR_MOVIE_COMPONENTS, compsPtr, found=componentsFound) + if (componentsFound) then + numberOfComponents = this%core%count(compsPtr) + allocate(res%cordinates(numberOfComponents)) + do i = 1, numberOfComponents + call this%core%get_child(compsPtr, i, compPtr) + call this%core%get(compPtr, component) + res%cordinates(i) = cs(1) + res%cordinates(i)%Or = buildVolProbeType(fieldType, component) + end do + else + allocate(res%cordinates(1)) + res%cordinates(1) = cs(1) + component = J_DIR_M + res%cordinates(1)%Or = buildVolProbeType(fieldType, component) + endif + res%len_cor = size(res%cordinates) + + res%outputrequest = trim(adjustl(this%getStrAt(p, J_NAME, default=" "))) + call setDomain(res, this%getDomain(p, J_PR_DOMAIN)) + end function + + integer function buildVolProbeType(fieldType, component) result(res) + character(len=:), allocatable, intent(in) :: fieldType, component + select case (fieldType) + case (J_FIELD_ELECTRIC) + select case (component) + case (J_DIR_X) + res = iExC + case (J_DIR_Y) + res = iEyC + case (J_DIR_Z) + res = iEzC + case (J_DIR_M) + res = iMEC + end select + case (J_FIELD_MAGNETIC) + select case (component) + case (J_DIR_X) + res = iHxC + case (J_DIR_Y) + res = iHyC + case (J_DIR_Z) + res = iHzC + case (J_DIR_M) + res = iMHC + end select + case (J_FIELD_CURRENT_DENSITY) + select case (component) + case (J_DIR_X) + res = iCurX + case (J_DIR_Y) + res = iCurY + case (J_DIR_Z) + res = iCurZ + case (J_DIR_M) + res = iCur + end select + case default + write(error_unit,*) "ERROR Determining vol probe type: invalid field type." + end select + end function + + subroutine setDomain(res, domain) + type(VolProbe), intent(inout) :: res + type(domain_t), intent(in) :: domain + + res%tstart = domain%tstart + res%tstep = domain%tstep + res%tstop = domain%tstop + res%fstart = domain%fstart + res%fstep = domain%fstep + res%fstop = domain%fstop + res%filename = domain%filename + res%type2 = domain%type2 + + if (domain%isLogarithmicFrequencySpacing) then + call appendLogSufix(res%outputrequest) + end if + end subroutine end function + subroutine appendLogSufix(fn) + character(len=BUFSIZE), intent(inout) :: fn + character (len=*), parameter :: SMBJSON_LOG_SUFFIX = "_log_" + fn = trim(fn) // SMBJSON_LOG_SUFFIX + end subroutine + function readThinWires(this) result (res) class(parser_t) :: this type(ThinWires) :: res @@ -1039,7 +1156,7 @@ function readGeneratorOnThinWire(linels, plineElemIds) result(res) if (.not. found) then return end if - + genSrcs = this%jsonValueFilterByKeyValues(sources, J_TYPE, [J_SRC_TYPE_GEN]) if (size(genSrcs) == 0) then return @@ -1054,27 +1171,27 @@ function readGeneratorOnThinWire(linels, plineElemIds) result(res) call this%core%get(genSrcs(i)%p, J_ELEMENTIDS, sourceElemIds) srcCoord = this%mesh%getNode(sourceElemIds(1)) polylineCoords = this%mesh%getPolyline(plineElemIds(1)) - if (.not. any(polylineCoords%coordIds == srcCoord%coordIds(1))) then + if (.not. any(polylineCoords%coordIds == srcCoord%coordIds(1))) then cycle ! generator is not in this polyline end if position = findSourcePositionInLinels(sourceElemIds, linels) - - if (.not. this%existsAt(genSrcs(i)%p, J_SRC_MAGNITUDE_FILE)) then + + if (.not. this%existsAt(genSrcs(i)%p, J_SRC_MAGNITUDE_FILE)) then write(error_unit, *) 'magnitudeFile of source missing' return end if select case(this%getStrAt(genSrcs(i)%p, J_FIELD)) - case (J_FIELD_VOLTAGE) + case (J_FIELD_VOLTAGE) res(position)%srctype = "VOLT" res(position)%srcfile = this%getStrAt(genSrcs(i)%p, J_SRC_MAGNITUDE_FILE) res(position)%multiplier = 1.0 - case (J_FIELD_CURRENT) + case (J_FIELD_CURRENT) res(position)%srctype = "CURR" res(position)%srcfile = this%getStrAt(genSrcs(i)%p, J_SRC_MAGNITUDE_FILE) res(position)%multiplier = 1.0 - case default + case default write(error_unit, *) 'Field block of source of type generator must be current or voltage' end select @@ -1089,13 +1206,15 @@ function findSourcePositionInLinels(srcElemIds, linels) result(res) type(pixel_t), dimension(:), allocatable :: pixels integer :: res integer :: i + pixels = this%mesh%convertNodeToPixels(this%mesh%getNode(srcElemIds(1))) do i = 1, size(linels) - if (all(linels(i)%cell ==pixels(1)%cell)) then - res = i + if (linels(i)%tag == pixels(1)%tag) then + res = i return end if end do + write (error_unit, * ) "ERROR: Source could not be found in linels." end function @@ -1123,11 +1242,11 @@ function readThinWireTermination(terminal) result(res) end if select case(label) - case(J_MAT_TERM_TYPE_OPEN) + case(J_MAT_TERM_TYPE_OPEN) res%r = 0.0 res%l = 0.0 res%c = 0.0 - case default + case default call this%core%get(tm, J_MAT_TERM_RESISTANCE, res%r, default=0.0) call this%core%get(tm, J_MAT_TERM_INDUCTANCE, res%l, default=0.0) call this%core%get(tm, J_MAT_TERM_CAPACITANCE, res%c, default=1e22) @@ -1181,12 +1300,6 @@ logical function isThinWire(cable) end function end function - function readSlantedWires(this) result (res) - class(parser_t) :: this - type(SlantedWires) :: res - ! TODO - end function - function getDomain(this, place, path) result(res) class(parser_t) :: this type(domain_t) :: res @@ -1211,7 +1324,7 @@ function getDomain(this, place, path) result(res) endif res%type1 = NP_T1_PLAIN - + call this%core%get(domain, J_TYPE, domainType) res%type2 = getNPDomainType(domainType, transferFunctionFound) @@ -1220,7 +1333,7 @@ function getDomain(this, place, path) result(res) call this%core%get(domain, J_PR_DOMAIN_TIME_STEP, res%tstep, default=0.0) call this%core%get(domain, J_PR_DOMAIN_FREQ_START, res%fstart, default=0.0) call this%core%get(domain, J_PR_DOMAIN_FREQ_STOP, res%fstop, default=0.0) - + call this%core%get(domain, J_PR_DOMAIN_FREQ_NUMBER, numberOfFrequencies, default=0) if (numberOfFrequencies == 0) then res%fstep = 0.0 @@ -1282,12 +1395,6 @@ function getNPDomainType(typeLabel, hasTransferFunction) result(res) end function end function - function readThinSlots(this) result (res) - class(parser_t) :: this - type(ThinSlots) :: res - ! TODO - end function - function readMTLN(this, grid) result (mtln_res) class(parser_t) :: this type(Desplazamiento), intent(in) :: grid @@ -1300,17 +1407,24 @@ function readMTLN(this, grid) result (mtln_res) mtln_res%number_of_steps = this%getRealAt(this%root, J_GENERAL//'.'//J_GEN_NUMBER_OF_STEPS) cables = readCables() - mtln_res%connectors => readConnectors() - call addConnIdToConnectorMap(connIdToConnector, mtln_res%connectors) block integer :: nW, nMW nMW = countNumberOfMultiwires(cables) - if (nMW == 0) return + if (nMW == 0) then + allocate(mtln_res%cables(0)) + allocate(mtln_res%probes(0)) + allocate(mtln_res%networks(0)) + return + end if nW = countNumberOfWires(cables) nWs = nW + nMW end block + mtln_res%connectors => readConnectors() + call addConnIdToConnectorMap(connIdToConnector, mtln_res%connectors) + + allocate (mtln_res%cables(nWs)) block logical :: is_read @@ -1352,7 +1466,7 @@ function readMTLN(this, grid) result (mtln_res) end do end if end block - + mtln_res%probes = readWireProbes() mtln_res%networks = buildNetworks() @@ -1477,8 +1591,8 @@ function buildNetworks() result(res) end do allocate(res(size(networks_coordinates))) - - + + do i = 1, size(networks_coordinates) res(i) = buildNetwork(networks_coordinates(i), aux_nodes) end do @@ -1567,7 +1681,7 @@ subroutine updateListOfNetworksCoordinates(coordinates, conductor_index) found_end = .false. polyline = this%mesh%getPolyline(conductor_index) coord_ini = this%mesh%getCoordinate(polyline%coordIds(1)) - + ub = ubound(polyline%coordIds,1) coord_end = this%mesh%getCoordinate(polyline%coordIds(ub)) @@ -1631,7 +1745,7 @@ function buildNode(termination_list, label, index, id) result(res) res%node%termination%resistance = readTerminationRLC(termination, J_MAT_TERM_RESISTANCE, default = 0.0) res%node%termination%inductance = readTerminationRLC(termination, J_MAT_TERM_INDUCTANCE, default=0.0) res%node%conductor_in_cable = index - + call elemIdToCable%get(key(id), value=cable_index) res%node%belongs_to_cable => mtln_res%cables(cable_index) @@ -1912,7 +2026,7 @@ function readMTLNCable(j_cable, is_read) result(res) end if res%step_size = buildStepSize(j_cable) - res%segment_relative_positions = mapSegmentsToGridCoordinates(j_cable) + res%external_field_segments = mapSegmentsToGridCoordinates(j_cable) ! write(*,*) 'id: ', this%getIntAt(j_cable, J_MATERIAL_ID, found) material = this%matTable%getId(this%getIntAt(j_cable, J_MATERIAL_ID, found)) if (.not. found) & @@ -2075,7 +2189,7 @@ subroutine assignPULProperties(res, mat, n) function mapSegmentsToGridCoordinates(j_cable) result(res) type(json_value), pointer :: j_cable - type(segment_relative_position_t), dimension(:), allocatable :: res + type(external_field_segment_t), dimension(:), allocatable :: res integer, dimension(:), allocatable :: elemIds type(polyline_t) :: p_line @@ -2088,49 +2202,55 @@ function mapSegmentsToGridCoordinates(j_cable) result(res) type(coordinate_t) :: c1, c2 integer :: i do i = 2, size(p_line%coordIds) - c2 = this%mesh%getCoordinate(p_line%coordIds(i)) - c1 = this%mesh%getCoordinate(p_line%coordIds(i-1)) - - if (findOrientation(c2-c1) > 0) then - res = [res, mapPositiveSegment(c1,c2)] - else if (findOrientation(c2-c1) < 0) then - res = [res, mapNegativeSegment(c1,c2)] - else - write(error_unit, *) 'Error: polyline first and last coordinate are identical' - end if + c2 = this%mesh%getCoordinate(p_line%coordIds(i)) + c1 = this%mesh%getCoordinate(p_line%coordIds(i-1)) + + if (findOrientation(c2-c1) > 0) then + res = [res, mapPositiveSegment(c1,c2)] + else if (findOrientation(c2-c1) < 0) then + res = [res, mapNegativeSegment(c1,c2)] + else + write(error_unit, *) 'Error: polyline first and last coordinate are identical' + end if end do end block end function function mapNegativeSegment(c1, c2) result(res) type(coordinate_t), intent(in) :: c1, c2 - type(segment_relative_position_t) :: curr_pos + type(external_field_segment_t) :: curr_pos integer :: axis, i, n_segments - type(segment_relative_position_t), dimension(:), allocatable :: res + type(external_field_segment_t), dimension(:), allocatable :: res axis = findDirection(c2-c1) n_segments = abs(ceiling(c2%position(axis)) - floor(c1%position(axis))) allocate(res(n_segments)) curr_pos%position = [(c1%position(i), i = 1, 3)] + curr_pos%Efield_main2wire => null() + curr_pos%Efield_wire2main => null() res = [(curr_pos, i = 1, n_segments)] res(:)%position(axis) = [(res(i)%position(axis) - i, i = 1, n_segments)] - + res(:)%direction = -axis end function function mapPositiveSegment(c1, c2) result(res) type(coordinate_t), intent(in) :: c1, c2 - type(segment_relative_position_t) :: curr_pos - integer :: axis, i, n_segments - type(segment_relative_position_t), dimension(:), allocatable :: res + type(external_field_segment_t) :: curr_pos + integer :: axis, orientation, i, n_segments + type(external_field_segment_t), dimension(:), allocatable :: res axis = findDirection(c2-c1) + n_segments = abs(floor(c2%position(axis)) - ceiling(c1%position(axis))) allocate(res(n_segments)) curr_pos%position = [(c1%position(i), i = 1, 3)] + curr_pos%Efield_main2wire => null() + curr_pos%Efield_wire2main => null() res = [(curr_pos, i = 1, n_segments)] res(:)%position(axis) = [(res(i)%position(axis) + (i-1), i = 1, n_segments)] + res(:)%direction = axis end function function buildStepSize(j_cable) result(res) @@ -2154,23 +2274,23 @@ function buildStepSize(j_cable) result(res) real :: f1, f2 real, dimension(:), allocatable :: displacement do j = 2, size(p_line%coordIds) - c2 = this%mesh%getCoordinate(p_line%coordIds(j)) - c1 = this%mesh%getCoordinate(p_line%coordIds(j-1)) - axis = findDirection(c2-c1) - f1 = abs(ceiling(c1%position(axis))-c1%position(axis)) - f2 = abs(c2%position(axis)-floor(c2%position(axis))) - displacement = assignDisplacement(desp, axis) - if (f1 /= 0) then - res = [res, f1*displacement(floor(c1%position(axis)))] - end if - index_1 = ceiling(min(abs(c1%position(axis)), abs(c2%position(axis)))) - index_2 = floor(max(abs(c1%position(axis)), abs(c2%position(axis)))) - do i = 1, index_2 - index_1 - res = [res, displacement(i)] - enddo - if (f2 /= 0) then - res = [res, f2*displacement(floor(c2%position(axis)))] - end if + c2 = this%mesh%getCoordinate(p_line%coordIds(j)) + c1 = this%mesh%getCoordinate(p_line%coordIds(j-1)) + axis = findDirection(c2-c1) + f1 = abs(ceiling(c1%position(axis))-c1%position(axis)) + f2 = abs(c2%position(axis)-floor(c2%position(axis))) + displacement = assignDisplacement(desp, axis) + if (f1 /= 0) then + res = [res, f1*displacement(floor(c1%position(axis)))] + end if + index_1 = ceiling(min(abs(c1%position(axis)), abs(c2%position(axis)))) + index_2 = floor(max(abs(c1%position(axis)), abs(c2%position(axis)))) + do i = 1, index_2 - index_1 + res = [res, displacement(i)] + enddo + if (f2 /= 0) then + res = [res, f2*displacement(floor(c2%position(axis)))] + end if end do end block end function @@ -2236,11 +2356,11 @@ function findOrientation(coordDiference) result(res) integer :: res integer :: i do i = 1, 3 - if (coordDiference%position(i) /= 0) then + if (coordDiference%position(i) /= 0) then res = coordDiference%position(i)/abs(coordDiference%position(i)) end if end do - end function + end function function findDirection(coordDiference) result(res) diff --git a/src_main_pub/fdetypes.F90 b/src_main_pub/fdetypes.F90 index 851c746e..6300e176 100755 --- a/src_main_pub/fdetypes.F90 +++ b/src_main_pub/fdetypes.F90 @@ -281,7 +281,8 @@ module FDETYPES NodalE , & NodalH , & PeriodicBorders, & - MagneticMedia, PMLMagneticMedia + MagneticMedia, PMLMagneticMedia, & + MTLNbundles end type !computational limits diff --git a/src_main_pub/interpreta_switches.F90 b/src_main_pub/interpreta_switches.F90 index d7cab987..2b99fc22 100755 --- a/src_main_pub/interpreta_switches.F90 +++ b/src_main_pub/interpreta_switches.F90 @@ -151,6 +151,7 @@ module interpreta_switches_m chain, & chain2, & fichin, & + extension, & nresumeable2, & fileFDE, & fileH5, & @@ -1901,6 +1902,7 @@ subroutine buscaswitchficheroinput(l) IF ((p-4) >= 1) THEN IF (f((p-4) :(p-4)) == NFDEEXTENSION(1:1)) THEN NFDEEXTENSION= f((p-4) :p) + l%extension=NFDEEXTENSION l%fichin = f (1:p-5) ELSE l%fichin = f (1:p) @@ -1992,6 +1994,7 @@ subroutine buscaswitchficheroinput(l) IF ((p-4) >= 1) THEN IF (f((p-4) :(p-4)) == NFDEEXTENSION(1:1)) THEN NFDEEXTENSION= f((p-4) :p) + l%extension=NFDEEXTENSION l%fichin = f (1:p-5) ELSE l%fichin = f (1:p) diff --git a/src_main_pub/nfde_types.F90 b/src_main_pub/nfde_types.F90 index d51b933b..50675fb0 100755 --- a/src_main_pub/nfde_types.F90 +++ b/src_main_pub/nfde_types.F90 @@ -802,1213 +802,3 @@ MODULE NFDETypes END MODULE NFDETypes - -!!AUX ROUTINES TO CREATE THE TYPES USED BY THE PRIVATE PARSER -MODULE BoxSourceClass - - ! - USE NFDETypes - ! -CONTAINS - ! - SUBROUTINE new_box (MPIDIR,this, nombre, i1, j1, k1, i2, j2, k2) - TYPE (Box), INTENT (INOUT) :: this - CHARACTER (LEN=*), INTENT (IN) :: nombre - INTEGER (KIND=4) :: i1, j1, k1, i2, j2, k2, mpidir - INTEGER (KIND=4) :: OXI, OXE, OYI, OYE,OZI, OZE - ! - this%coor1 (1) = i1 - this%coor1 (2) = j1 - this%coor1 (3) = k1 - this%coor2 (1) = i2 - this%coor2 (2) = j2 - this%coor2 (3) = k2 - this%nombre_fichero = nombre - !MPI ROTATE BOX - IF (MPIDIR==2 ) THEN - OXI= this%coor1 (1) - OXE= this%coor2 (1) - OYI= this%coor1 (2) - OYE= this%coor2 (2) - OZI= this%coor1 (3) - OZE= this%coor2 (3) - ! - this%coor1 (1) =OZI - this%coor2 (1) =OZE - this%coor1 (2) =OXI - this%coor2 (2) =OXE - this%coor1 (3) =OYI - this%coor2 (3) =OYE - ELSEIF (MPIDIR==1 ) THEN - OXI= this%coor1 (1) - OXE= this%coor2 (1) - OYI= this%coor1 (2) - OYE= this%coor2 (2) - OZI= this%coor1 (3) - OZE= this%coor2 (3) - ! - this%coor1 (1) =OYI - this%coor2 (1) =OYE - this%coor1 (2) =OZI - this%coor2 (2) =OZE - this%coor1 (3) =OXI - this%coor2 (3) =OXE - ENDIF - !!!!! - ! - RETURN - END SUBROUTINE new_box - ! - FUNCTION incluirBox (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (Box), INTENT (IN) :: n_Volume - TYPE (Box), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (Box), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION incluirBox - ! - -END MODULE BoxSourceClass -!------------------------------------------------------------------------------ -! File : BloqueProbe.f90 -! Created on April 24, 2009, 12 : 52 AM -! This MODULE defines the class to deal with the Bloque probes in the -! nfde input files. It encodes the information provided by it and provides -! methods to generate them as outputs. -!------------------------------------------------------------------------------ -MODULE BloqueProbeClass - - ! - USE NFDETypes - ! - IMPLICIT NONE - ! -CONTAINS - !------------------------------------------------------------------------------ - ! Add probes to the collection of probes object - !------------------------------------------------------------------------------ - FUNCTION addBloqueCurrentProbes (n_BloqueProbe, o_BloqueProbes, steps) RESULT (BloqueProbes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (BloqueProbe), INTENT (INOUT) :: n_BloqueProbe - TYPE (BloqueProbe), DIMENSION (steps-1), INTENT (IN) :: o_BloqueProbes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (BloqueProbe), DIMENSION (:), POINTER :: BloqueProbes - !--> Allocate space for the new array - ALLOCATE (BloqueProbes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - BloqueProbes (n) = o_BloqueProbes (n) - END DO - !--> Add new value found - BloqueProbes (steps) = n_BloqueProbe - ! - RETURN - END FUNCTION addBloqueCurrentProbes - ! - -END MODULE BloqueProbeClass - -!------------------------------------------------------------------------------ -! File : Frontera.f90 -! Created on March 16, 2009, 13 : 03 AM -! Definicin de la clase frontera definida para trabajar con las condiciones -! frontera del fichero nfde -!------------------------------------------------------------------------------ -MODULE FronteraClass - - USE NFDETypes -CONTAINS - !------------------------------------------------------------------------------ - ! Mtodo que establece el tipo de frontera para la posicin sealada - ! introduciendo el tipo sealado - !------------------------------------------------------------------------------ - SUBROUTINE establecerFrontera (this, posicion, tipo) - TYPE (Frontera), INTENT (INOUT) :: this - INTEGER (KIND=4) :: posicion, tipo - this%tipoFrontera (posicion) = tipo - END SUBROUTINE establecerFrontera - !------------------------------------------------------------------------------ - ! Mtodo para asignar propiedades en la posicin indicadada - ! con el valor dado - !------------------------------------------------------------------------------ - SUBROUTINE asignarPropiedades (this, valor, posicion) - TYPE (Frontera), INTENT (INOUT) :: this - INTEGER (KIND=4) :: posicion, valor - this%tipoFrontera (posicion) = valor - END SUBROUTINE asignarPropiedades - !------------------------------------------------------------------------------ - ! Mtodo por el cual se establecen las propiedades PML en la posicin - ! sealada por el mtodo - !------------------------------------------------------------------------------ - SUBROUTINE asignarPropiedadesPML (this, pml, posicion) - TYPE (Frontera), INTENT (INOUT) :: this - TYPE (FronteraPML) :: pml - INTEGER (KIND=4) :: posicion - this%propiedadesPML (posicion) = pml - END SUBROUTINE asignarPropiedadesPML - !------------------------------------------------------------------------------ - ! Constructor de la clase de propiedades del PML - ! necesario introducir siempre los datos - !------------------------------------------------------------------------------ - SUBROUTINE new_pml (this, numCapas, orden, refl) - TYPE (FronteraPML), INTENT (INOUT) :: this - REAL (KIND=RK) :: orden, refl - REAL (KIND=RK) :: numCapas - this%numCapas = nint(numCapas) !por algun misterio ORIGINAL pone las capas como reales - this%orden = orden - this%refl = refl - END SUBROUTINE new_pml - ! - -END MODULE FronteraClass -!------------------------------------------------------------------------------ -! File : MasSonda.f90 -! Created on April 24, 2009, 12 : 52 AM -! This MODULE defines the class to deal with the new probes in the -! nfde input files. It encodes the information provided by it and provides -! methods to generate them as outputs. -! Despite it not being the most efficient solution, the following code uses -! dynamic dispatching -!------------------------------------------------------------------------------ -! -MODULE MasSondaClass - - ! - USE NFDETypes - IMPLICIT NONE - ! -CONTAINS - !------------------------------------------------------------------------------ - ! SUBROUTINE to generate the construction of the - ! new probe class - !------------------------------------------------------------------------------ - SUBROUTINE new_MasSonda (this, type1, type2, filename, tstart, tstop, tstep, fstart, fstop, fstep, outputrequest) - TYPE (MasSonda), INTENT (INOUT) :: this - INTEGER (KIND=4), INTENT (INOUT) :: type1, type2 - CHARACTER (LEN=*), INTENT (INOUT) :: filename - REAL (KIND=RK), INTENT (INOUT) :: tstart, tstop, tstep - REAL (KIND=RK), INTENT (INOUT) :: fstart, fstop, fstep - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - this%type1 = type1 - this%type2 = type2 - this%filename = filename - this%tstart = tstart - this%tstop = tstop - this%tstep = tstep - this%fstart = fstart - this%fstop = fstop - this%fstep = fstep - this%len_cor = 0 - NULLIFY (this%cordinates) - ! - RETURN - END SUBROUTINE new_MasSonda - ! - -END MODULE MasSondaClass -! -!------------------------------------------------------------------------------ -! File : MasSonda.f90 -! Created on April 24, 2009, 12 : 52 AM -! This MODULE defines the class to deal with the new probes in the -! nfde input files. It encodes the information provided by it and provides -! methods to generate them as outputs. -! Despite it not being the most efficient solution, the following code uses -! dynamic dispatching -!------------------------------------------------------------------------------ -! -MODULE VolProbeClass - - ! - USE NFDETypes - IMPLICIT NONE - ! -CONTAINS - !------------------------------------------------------------------------------ - ! SUBROUTINE to generate the construction of the - ! Scr probe class - !------------------------------------------------------------------------------ - SUBROUTINE New_VolProbe (this, tstart, tstop, tstep, fstart, fstop, fstep, type2, filename, outputrequest) - TYPE (VolProbe), INTENT (INOUT) :: this - REAL (KIND=RK), INTENT (INOUT) :: tstart, tstop, tstep,fstart, fstop, fstep - !sgg'10 - INTEGER (KIND=4) :: comi,type2 - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - CHARACTER (LEN=BUFSIZE) :: filename - comi = index (outputrequest, '*') + 1 - this%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - this%tstart = tstart - this%tstop = tstop - this%tstep = tstep - this%len_cor = 0 - - this%fstart=fstart - this%fstop=fstop - this%fstep=fstep - this%type2=type2 - this%filename=filename - - NULLIFY (this%cordinates) - ! - RETURN - END SUBROUTINE New_VolProbe - ! - -END MODULE VolProbeClass -!------------------------------------------------------------------------------ -! File : planeWave.f90 -! Created on March 19, 2009, 10 : 59 AM -!------------------------------------------------------------------------------ -! -MODULE PlaneWaveClass - - ! - USE NFDETypes - ! -CONTAINS - !------------------------------------------------------------------------------ - !------------------------------------------------------------------------------ - !------------------------------------------------------------------------------ - SUBROUTINE new_planewave (mpidir,this, nombre, attribute, i1, j1, k1, i2, j2, k2, theta, phi, alpha, beta,isRC,numModes,INCERTMAX) - TYPE (PlaneWave), INTENT (INOUT) :: this - CHARACTER (LEN=*), INTENT (IN) :: nombre, attribute - INTEGER (KIND=4), INTENT (IN) :: i1, j1, k1, i2, j2, k2,numModes - logical, INTENT (IN) :: isRC - REAL (KIND=RK), INTENT (IN) :: theta, phi, alpha, beta, INCERTMAX - INTEGER (KIND=4) :: OXI, OXE, OYI, OYE,OZI, OZE,mpidir - ! - this%coor1 (1) = i1 - this%coor1 (2) = j1 - this%coor1 (3) = k1 - this%coor2 (1) = i2 - this%coor2 (2) = j2 - this%coor2 (3) = k2 - this%nombre_fichero = nombre - this%atributo = attribute - this%theta = theta - this%phi = phi - this%alpha = alpha - this%beta = beta - this%isRC = isRC - this%INCERTMAX = INCERTMAX - this%numModes = numModes - - !MPI ROTATE PLANEWAVE - IF (MPIDIR==2 ) THEN - OXI= this%coor1 (1) - OXE= this%coor2 (1) - OYI= this%coor1 (2) - OYE= this%coor2 (2) - OZI= this%coor1 (3) - OZE= this%coor2 (3) - ! - this%coor1 (1) =OZI - this%coor2 (1) =OZE - this%coor1 (2) =OXI - this%coor2 (2) =OXE - this%coor1 (3) =OYI - this%coor2 (3) =OYE - - this%theta = atan2(Sqrt(Cos(theta)**2.0_RKIND+ Cos(phi)**2*Sin(theta)**2),Sin(phi)*Sin(theta)) - this%phi = atan2(Cos(phi)*Sin(theta),Cos(theta)) - this%alpha = atan2(Sqrt(Cos(alpha)**2.0_RKIND+ Cos(beta)**2*Sin(alpha)**2),Sin(beta)*Sin(alpha)) - this%beta = atan2(Cos(beta)*Sin(alpha),Cos(alpha)) - - ELSEIF (MPIDIR==1 ) THEN - OXI= this%coor1 (1) - OXE= this%coor2 (1) - OYI= this%coor1 (2) - OYE= this%coor2 (2) - OZI= this%coor1 (3) - OZE= this%coor2 (3) - ! - this%coor1 (1) =OYI - this%coor2 (1) =OYE - this%coor1 (2) =OZI - this%coor2 (2) =OZE - this%coor1 (3) =OXI - this%coor2 (3) =OXE - - this%theta = atan2(Sqrt(Cos(theta)**2.0_RKIND+ Sin(phi)**2*Sin(theta)**2),Cos(phi)*Sin(theta)) - this%phi = atan2(Cos(theta),Sin(phi)*Sin(theta)) - this%alpha = atan2(Sqrt(Cos(alpha)**2.0_RKIND+ Sin(beta)**2*Sin(alpha)**2),Cos(beta)*Sin(alpha)) - this%beta = atan2(Cos(alpha),Sin(beta)*Sin(alpha)) - ENDIF - !!!!! - RETURN - END SUBROUTINE new_planewave - ! - FUNCTION incluirPlaneWave (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (PlaneWave), INTENT (IN) :: n_Volume - TYPE (PlaneWave), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (PlaneWave), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION incluirPlaneWave - ! - -END MODULE PlaneWaveClass -!------------------------------------------------------------------------------ -! File : ThinSlots.f90 -! MODULE to cover the information provided by thin Slots in the development of -! the parsing application -!------------------------------------------------------------------------------ -MODULE ThinSlotsClass - - ! - USE NFDETypes - ! -CONTAINS - !-------------------------------------------------------------------------- - ! Adds a ThinSlotComponent inside the ThinSlot object given as a - ! parameter - !-------------------------------------------------------------------------- - FUNCTION addThinSlotComp (n_thinSlotComp, o_thinSlotsComp, steps) RESULT (thinSlotsComp) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (ThinSlotComp), INTENT (IN) :: n_thinSlotComp - TYPE (ThinSlotComp), DIMENSION (steps-1), INTENT (IN) :: o_thinSlotsComp - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (ThinSlotComp), DIMENSION (:), POINTER :: thinSlotsComp - !--> Allocate space for the new array - ALLOCATE (thinSlotsComp(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - thinSlotsComp (n) = o_thinSlotsComp (n) - END DO - !--> Add new value found - thinSlotsComp (steps) = n_thinSlotComp - ! - RETURN - END FUNCTION addThinSlotComp - !-------------------------------------------------------------------------- - ! Adds a ThinSlot inside the ThinSlots object given as a parameter - !-------------------------------------------------------------------------- - FUNCTION addThinSlot (n_thinSlot, o_thinSlots, steps) RESULT (ThinSlots) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (ThinSlot), INTENT (IN) :: n_thinSlot - TYPE (ThinSlot), DIMENSION (steps-1), INTENT (IN) :: o_thinSlots - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (ThinSlot), DIMENSION (:), POINTER :: ThinSlots - !--> Allocate space for the new array - ALLOCATE (ThinSlots(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - ThinSlots (n) = o_thinSlots (n) - END DO - !--> Add new value found - ThinSlots (steps) = n_thinSlot - ! - RETURN - END FUNCTION addThinSlot - ! - -END MODULE ThinSlotsClass -!------------------------------------------------------------------------------ -! File : ThinWires.f90 -! Created on August 10, 2009, 12 : 27 AM -! MODULE to cover the information provided by thin wires in the development of -! the parsing application -!------------------------------------------------------------------------------ -MODULE ThinWiresClass - - ! - USE NFDETypes - !---------------------------------------------------------------------------- - ! Adds a ThinWireComponent inside the ThinWire object given as a - ! parameter - !---------------------------------------------------------------------------- -CONTAINS - FUNCTION addThinWireComp (n_thinWireComp, o_thinWiresComp, steps) RESULT (thinWiresComp) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (ThinWireComp), INTENT (IN) :: n_thinWireComp - TYPE (ThinWireComp), DIMENSION (steps-1), INTENT (IN) :: o_thinWiresComp - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (ThinWireComp), DIMENSION (:), POINTER :: thinWiresComp - !--> Allocate space for the new array - ALLOCATE (thinWiresComp(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - thinWiresComp (n) = o_thinWiresComp (n) - END DO - !--> Add new value found - thinWiresComp (steps) = n_thinWireComp - ! - RETURN - END FUNCTION addThinWireComp - !---------------------------------------------------------------------------- - ! Adds a ThinWireComponent inside the ThinWire object given as a - ! parameter - !---------------------------------------------------------------------------- - FUNCTION addThinWire (n_thinWire, o_thinWires, steps) RESULT (ThinWires) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (ThinWire), INTENT (IN) :: n_thinWire - TYPE (ThinWire), DIMENSION (steps-1), INTENT (IN) :: o_thinWires - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (ThinWire), DIMENSION (:), POINTER :: ThinWires - !--> Allocate space for the new array - ALLOCATE (ThinWires(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - ThinWires (n) = o_thinWires (n) - END DO - !--> Add new value found - ThinWires (steps) = n_thinWire - ! - RETURN - END FUNCTION addThinWire - ! - -END MODULE ThinWiresClass -!------------------------------------------------------------------------------ -! File : SlantedWires.f90 -! Created on September 28, 2015, 18 : 51 AM -! MODULE to cover the information provided by slanted wires in the development of -! the parsing application -!------------------------------------------------------------------------------ -MODULE SlantedWiresClass - - ! - USE NFDETypes - !---------------------------------------------------------------------------- - ! Adds a SlantedWireComponent inside the SlantedWire object given as a - ! parameter - !---------------------------------------------------------------------------- -CONTAINS - FUNCTION addSlantedWireComp (n_slantedWireComp, o_slantedWiresComp, steps) RESULT (slantedWiresComp) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (SlantedWireComp), INTENT (IN) :: n_slantedWireComp - TYPE (SlantedWireComp), DIMENSION (steps-1), INTENT (IN) :: o_slantedWiresComp - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (SlantedWireComp), DIMENSION (:), POINTER :: slantedWiresComp - !--> Allocate space for the new array - ALLOCATE (slantedWiresComp(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - slantedWiresComp (n) = o_slantedWiresComp (n) - END DO - !--> Add new value found - slantedWiresComp (steps) = n_slantedWireComp - ! - RETURN - END FUNCTION addSlantedWireComp - !---------------------------------------------------------------------------- - ! Adds a SlantedWireComponent inside the SlantedWire object given as a - ! parameter - !---------------------------------------------------------------------------- - FUNCTION addSlantedWire (n_slantedWire, o_slantedWires, steps) RESULT (SlantedWires) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (SlantedWire), INTENT (IN) :: n_slantedWire - TYPE (SlantedWire), DIMENSION (steps-1), INTENT (IN) :: o_slantedWires - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (SlantedWire), DIMENSION (:), POINTER :: SlantedWires - !--> Allocate space for the new array - ALLOCATE (SlantedWires(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - SlantedWires (n) = o_slantedWires (n) - END DO - !--> Add new value found - SlantedWires (steps) = n_slantedWire - ! - RETURN - END FUNCTION addSlantedWire - ! - -END MODULE SlantedWiresClass -!------------------------------------------------------------------------------ -! File : SondaClass.f90 -! Created on April 24, 2009, 12 : 52 AM -! This MODULE defines the class to deal with the probes in the -! nfde input files. It encodes the information provided by it and provides -! methods to generate them as outputs. -! Despite it not being the most efficient solution, the following code uses -! dynamic dispatching -!------------------------------------------------------------------------------ -! -MODULE SondaClass - - !------------------------------------------------------------------------------ - ! External classes used - !------------------------------------------------------------------------------ - ! - USE NFDETypes - !------------------------------------------------------------------------------ - ! PRIVATE methods - !------------------------------------------------------------------------------ - PRIVATE :: FarField_Sonda_new - PRIVATE :: Electric_Sonda_new, Magnetic_Sonda_new - PRIVATE :: NormalElectric_Sonda_new, NormalMagnetic_Sonda_new - PRIVATE :: SurfaceElectricCurrent_Sonda_new, SurfaceMagneticCurrent_Sonda_new - PRIVATE :: assign_FarField,assign_Electric, assign_Magnetic, assign_NormalElectric, assign_NormalMagnetic, assign_SurfaceMagneticCurrent, assign_SurfaceElectricCurrent - PRIVATE :: addCoordinates_1, addCoordinates_2 - !------------------------------------------------------------------------------ - ! INTERFACE for the constructor of classes - !------------------------------------------------------------------------------ - INTERFACE new_SONDA - MODULE PROCEDURE FarField_Sonda_new - MODULE PROCEDURE Electric_Sonda_new - MODULE PROCEDURE Magnetic_Sonda_new - MODULE PROCEDURE NormalElectric_Sonda_new - MODULE PROCEDURE NormalMagnetic_Sonda_new - MODULE PROCEDURE SurfaceElectricCurrent_Sonda_new - MODULE PROCEDURE SurfaceMagneticCurrent_Sonda_new - MODULE PROCEDURE abstractProbe_new - END INTERFACE new_SONDA - !------------------------------------------------------------------------------ - ! INTERFACE for assigning to the abstract class a non-abstract - ! object - !------------------------------------------------------------------------------ - INTERFACE assign - MODULE PROCEDURE assign_FarField - MODULE PROCEDURE assign_Electric - MODULE PROCEDURE assign_Magnetic - MODULE PROCEDURE assign_NormalElectric - MODULE PROCEDURE assign_NormalMagnetic - MODULE PROCEDURE assign_SurfaceElectricCurrent - MODULE PROCEDURE assign_SurfaceMagneticCurrent - END INTERFACE assign - !------------------------------------------------------------------------------ - ! INTERFACE to distinguish from different add coordinates - !------------------------------------------------------------------------------ - INTERFACE addCoordinates - MODULE PROCEDURE addCoordinates_1 - MODULE PROCEDURE addCoordinates_2 - END INTERFACE addCoordinates - !------------------------------------------------------------------------------ - ! Begining of the subroutines and functions - !------------------------------------------------------------------------------ -CONTAINS - !------------------------------------------------------------------------------ - ! Initialices the object for abstract probes - !------------------------------------------------------------------------------ - SUBROUTINE abstractProbe_new (this) - TYPE (abstractSonda), INTENT (INOUT) :: this - ! - this%n_FarField = 0 - this%n_Electric = 0 - this%n_Magnetic = 0 - this%n_NormalElectric = 0 - this%n_NormalMagnetic = 0 - this%n_SurfaceElectricCurrent = 0 - this%n_SurfaceMagneticCurrent = 0 - this%n_FarField_max = 0 - this%n_Electric_max = 0 - this%n_Magnetic_max = 0 - this%n_NormalElectric_max = 0 - this%n_NormalMagnetic_max = 0 - this%n_SurfaceElectricCurrent_max = 0 - this%n_SurfaceMagneticCurrent_max = 0 - NULLIFY (this%FarField,this%Electric, this%Magnetic, this%NormalElectric, this%NormalMagnetic, this%SurfaceElectricCurrent, this%SurfaceMagneticCurrent) - ! - RETURN - END SUBROUTINE abstractProbe_new - !------------------------------------------------------------------------------ - ! Constructor for the FAR electric field probe - !------------------------------------------------------------------------------ - SUBROUTINE FarField_Sonda_new (MPIDIR,this, name, i, j, K, node, fstart, fstop, fstep, & - thetastart, thetastop, thetastep,phistart, phistop, phistep, FileNormalize,& - outputrequest) - - TYPE (FarField_Sonda), INTENT (INOUT) :: this - CHARACTER (LEN=*) :: name,FileNormalize - INTEGER (KIND=4) :: i, j, K, node,MPIDIR - REAL (KIND=RK) :: fstart, fstop, fstep - REAL (KIND=RK) :: thetastart, thetastop, thetastep,phistart, phistop, phistep - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%probe%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - ! - this%probe%grname = name - this%probe%fstart = fstart - this%probe%fstop = fstop - this%probe%fstep = fstep - this%probe%thetastart = thetastart - this%probe%thetastop = thetastop - this%probe%thetastep = thetastep - this%probe%phistart = phistart - this%probe%phistop = phistop - this%probe%phistep = phistep - - !!!mpirotate angulos farfield .... las coordenadas se rotan luego - IF (MPIDIR==2 ) THEN - this%probe%thetastart = atan2(Sqrt(Cos(thetastart)**2.0_RKIND+ Cos(phistart)**2*Sin(thetastart)**2),Sin(phistart)*Sin(thetastart)) - this%probe%phistart = atan2(Cos(phistart)*Sin(thetastart),Cos(thetastart)) - ELSEIF (MPIDIR==1 ) THEN - this%probe%thetastart = atan2(Sqrt(Cos(thetastart)**2.0_RKIND+ Sin(phistart)**2*Sin(thetastart)**2),Cos(phistart)*Sin(thetastart)) - this%probe%phistart = atan2(Cos(thetastart),Sin(phistart)*Sin(thetastart)) - ENDIF - !!!!!! - - - this%probe%FileNormalize=trim(adjustl(FileNormalize)) - - this%probe%n_cord = 1 - this%probe%n_cord_max = 1 - NULLIFY (this%probe%i, this%probe%j, this%probe%K, this%probe%node) - CALL addSondaCoordinates (MPIDIR,this%probe, i, j, K, node, this%probe%n_cord) - ! - RETURN - END SUBROUTINE FarField_Sonda_new - !------------------------------------------------------------------------------ - ! Constructor for the electric field probe - !------------------------------------------------------------------------------ - SUBROUTINE Electric_Sonda_new (MPIDIR,this, name, i, j, K, node, tstart, tstop, tstep, outputrequest) - TYPE (Electric_Sonda), INTENT (INOUT) :: this - CHARACTER (LEN=*) :: name - INTEGER (KIND=4) :: MPIDIR,i, j, K, node - REAL (KIND=RK) :: tstart, tstop, tstep - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%probe%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - ! - this%probe%grname = name - this%probe%tstart = tstart - this%probe%tstop = tstop - this%probe%tstep = tstep - this%probe%n_cord = 1 - this%probe%n_cord_max = 1 - NULLIFY (this%probe%i, this%probe%j, this%probe%K, this%probe%node) - CALL addSondaCoordinates (MPIDIR,this%probe, i, j, K, node, this%probe%n_cord) - ! - RETURN - END SUBROUTINE Electric_Sonda_new - !------------------------------------------------------------------------------ - ! Constructor for the magnetic field probe - !------------------------------------------------------------------------------ - SUBROUTINE Magnetic_Sonda_new (MPIDIR,this, name, i, j, K, node, tstart, tstop, tstep, outputrequest) - TYPE (Magnetic_Sonda), INTENT (INOUT) :: this - CHARACTER (LEN=*), INTENT (IN) :: name - INTEGER (KIND=4) :: MPIDIR,i, j, K, node - REAL (KIND=RK) :: tstart, tstop, tstep - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%probe%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - ! - this%probe%grname = name - this%probe%tstart = tstart - this%probe%tstop = tstop - this%probe%tstep = tstep - this%probe%n_cord = 1 - this%probe%n_cord_max = 1 - NULLIFY (this%probe%i, this%probe%j, this%probe%K, this%probe%node) - CALL addSondaCoordinates (MPIDIR,this%probe, i, j, K, node, this%probe%n_cord) - ! - RETURN - END SUBROUTINE Magnetic_Sonda_new - !------------------------------------------------------------------------------ - ! Constructor for the normal electric field probe - !------------------------------------------------------------------------------ - SUBROUTINE NormalElectric_Sonda_new (MPIDIR,this, name, i, j, K, node, nml, tstart, tstop, tstep, outputrequest) - TYPE (NormalElectric_Sonda), INTENT (INOUT) :: this - CHARACTER (LEN=*), INTENT (IN) :: name - INTEGER (KIND=4) :: MPIDIR,i, j, K, node, nml - REAL (KIND=RK) :: tstart, tstop, tstep - INTEGER (KIND=4), DIMENSION (:), POINTER :: aux_f - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%probe%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - ! - this%n_nml = 1 - NULLIFY (this%nml) - this%probe%grname = name - this%probe%tstart = tstart - this%probe%tstop = tstop - this%probe%tstep = tstep - this%probe%n_cord = 1 - this%probe%n_cord_max = 1 - NULLIFY (this%probe%i, this%probe%j, this%probe%K, this%probe%node) - CALL addSondaCoordinates (MPIDIR,this%probe, i, j, K, node, this%probe%n_cord) - aux_f => addNml (nml, this%nml, this%n_nml) - IF (ASSOCIATED(this%nml)) DEallocate (this%nml) - this%nml => aux_f - ! - RETURN - END SUBROUTINE NormalElectric_Sonda_new - !------------------------------------------------------------------------------ - ! Constructor for the magnetic current surface probe object - !------------------------------------------------------------------------------ - SUBROUTINE SurfaceMagneticCurrent_Sonda_new (MPIDIR,this, name, i, j, K, node, nml, tstart, tstop, tstep, outputrequest) - TYPE (SurfaceMagneticCurrent_Sonda), INTENT (INOUT) :: this - CHARACTER (LEN=*) :: name - INTEGER (KIND=4) :: MPIDIR,i, j, K, node, nml - REAL (KIND=RK) :: tstart, tstop, tstep - INTEGER (KIND=4), DIMENSION (:), POINTER :: aux_f - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%probe%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - ! - this%n_nml = 1 - NULLIFY (this%nml) - this%probe%grname = name - this%probe%tstart = tstart - this%probe%tstop = tstop - this%probe%tstep = tstep - this%probe%n_cord = 1 - this%probe%n_cord_max = 1 - NULLIFY (this%probe%i, this%probe%j, this%probe%K, this%probe%node) - CALL addSondaCoordinates (MPIDIR,this%probe, i, j, K, node, this%probe%n_cord) - aux_f => addNml (nml, this%nml, this%n_nml) - IF (ASSOCIATED(this%nml)) DEallocate (this%nml) - this%nml => aux_f - ! - RETURN - END SUBROUTINE SurfaceMagneticCurrent_Sonda_new - !------------------------------------------------------------------------------ - ! Constructor for the normal magnetic field probe - !------------------------------------------------------------------------------ - SUBROUTINE NormalMagnetic_Sonda_new (MPIDIR,this, name, i, j, K, node, nml, tstart, tstop, tstep, outputrequest) - TYPE (NormalMagnetic_Sonda), INTENT (INOUT) :: this - CHARACTER (LEN=*) :: name - INTEGER (KIND=4) :: MPIDIR,i, j, K, node, nml - REAL (KIND=RK) :: tstart, tstop, tstep - INTEGER (KIND=4), DIMENSION (:), POINTER :: aux_f - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%probe%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - ! - this%n_nml = 1 - NULLIFY (this%nml) - this%probe%grname = name - this%probe%tstart = tstart - this%probe%tstop = tstop - this%probe%tstep = tstep - this%probe%n_cord = 1 - this%probe%n_cord_max = 1 - NULLIFY (this%probe%i, this%probe%j, this%probe%K, this%probe%node) - CALL addSondaCoordinates (MPIDIR,this%probe, i, j, K, node, this%probe%n_cord) - aux_f => addNml (nml, this%nml, this%n_nml) - IF (ASSOCIATED(this%nml)) DEallocate (this%nml) - this%nml => aux_f - ! - RETURN - END SUBROUTINE NormalMagnetic_Sonda_new - !------------------------------------------------------------------------------ - ! Constructor for the electric current surface probe object - !------------------------------------------------------------------------------ - SUBROUTINE SurfaceElectricCurrent_Sonda_new (MPIDIR,this, name, i, j, K, node, nml, tstart, tstop, tstep, outputrequest) - TYPE (SurfaceElectricCurrent_Sonda), INTENT (INOUT) :: this - CHARACTER (LEN=*) :: name - INTEGER (KIND=4) :: MPIDIR,i, j, K, node, nml - REAL (KIND=RK) :: tstart, tstop, tstep - INTEGER (KIND=4), DIMENSION (:), POINTER :: aux_f - !sgg'10 - INTEGER (KIND=4) :: comi - CHARACTER (LEN=*), INTENT (INOUT) :: outputrequest - comi = index (outputrequest, '*') + 1 - this%probe%outputrequest = trim (adjustl(outputrequest(comi:))) - !! fin sgg'10 - ! - this%n_nml = 1 - NULLIFY (this%nml) - this%probe%grname = name - this%probe%tstart = tstart - this%probe%tstop = tstop - this%probe%tstep = tstep - this%probe%n_cord = 1 - this%probe%n_cord_max = 1 - NULLIFY (this%probe%i, this%probe%j, this%probe%K, this%probe%node) - CALL addSondaCoordinates (MPIDIR,this%probe, i, j, K, node, this%probe%n_cord) - aux_f => addNml (nml, this%nml, this%n_nml) - IF (ASSOCIATED(this%nml)) DEallocate (this%nml) - this%nml => aux_f - ! - RETURN - END SUBROUTINE SurfaceElectricCurrent_Sonda_new - !------------------------------------------------------------------------------ - ! Add probes to the list of probes - !------------------------------------------------------------------------------ - FUNCTION addSONDA (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (abstractSonda), INTENT (IN) :: n_Volume - TYPE (abstractSonda), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (abstractSonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - !--> Assigning new values - ! IF ( ASSOCIATED( o_Volumes ) ) DEallocate (o_Volumes) - ! - RETURN - END FUNCTION addSONDA - !------------------------------------------------------------------------------ - ! Selector to add coordinates - !------------------------------------------------------------------------------ - SUBROUTINE addCoordinates_1 (MPIDIR,this, i, j, K, node, nml, precounting) - TYPE (abstractSonda), INTENT (INOUT) :: this - INTEGER (KIND=4), INTENT (INOUT) :: MPIDIR,i, j, K, node, nml - INTEGER (KIND=4) :: precounting - ! - IF (ASSOCIATED(this%NormalElectric)) THEN - this%NormalElectric(1)%probe%n_cord_max = this%NormalElectric(1)%probe%n_cord_max + 1 - this%NormalElectric(1)%probe%n_cord = this%NormalElectric(1)%probe%n_cord * precounting + 1 - CALL addSondaCoordinates (MPIDIR,this%NormalElectric(1)%probe, i, j, K, node, this%NormalElectric(1)%probe%n_cord) - this%NormalElectric(1)%n_nml_max = this%NormalElectric(1)%n_nml_max + 1 - this%NormalElectric(1)%n_nml = this%NormalElectric(1)%n_nml * precounting + 1 - this%NormalElectric(1)%nml => addNml (nml, this%NormalElectric(1)%nml, this%NormalElectric(1)%n_nml) - ELSE IF (ASSOCIATED(this%NormalMagnetic)) THEN - this%NormalMagnetic(1)%probe%n_cord_max = this%NormalMagnetic(1)%probe%n_cord_max + 1 - this%NormalMagnetic(1)%probe%n_cord = this%NormalMagnetic(1)%probe%n_cord * precounting + 1 - CALL addSondaCoordinates (MPIDIR,this%NormalMagnetic(1)%probe, i, j, K, node, this%NormalMagnetic(1)%probe%n_cord) - this%NormalMagnetic(1)%n_nml_max = this%NormalMagnetic(1)%n_nml_max + 1 - this%NormalMagnetic(1)%n_nml = this%NormalMagnetic(1)%n_nml * precounting + 1 - this%NormalMagnetic(1)%nml => addNml (nml, this%NormalMagnetic(1)%nml, this%NormalMagnetic(1)%n_nml) - ELSE IF (ASSOCIATED(this%SurfaceElectricCurrent)) THEN - this%SurfaceElectricCurrent(1)%probe%n_cord_max = this%SurfaceElectricCurrent(1)%probe%n_cord_max + 1 - this%SurfaceElectricCurrent(1)%probe%n_cord = this%SurfaceElectricCurrent(1)%probe%n_cord * precounting + 1 - CALL addSondaCoordinates (MPIDIR,this%SurfaceElectricCurrent(1)%probe, i, j, K, node, this%SurfaceElectricCurrent(1)%probe%n_cord) - this%SurfaceElectricCurrent(1)%n_nml_max = this%SurfaceElectricCurrent(1)%n_nml_max + 1 - this%SurfaceElectricCurrent(1)%n_nml = this%SurfaceElectricCurrent(1)%n_nml * precounting + 1 - this%SurfaceElectricCurrent(1)%nml => addNml (nml, this%SurfaceElectricCurrent(1)%nml, this%SurfaceElectricCurrent(1)%n_nml) - ELSE IF (ASSOCIATED(this%SurfaceMagneticCurrent)) THEN - this%SurfaceMagneticCurrent(1)%probe%n_cord_max = this%SurfaceMagneticCurrent(1)%probe%n_cord_max + 1 - this%SurfaceMagneticCurrent(1)%probe%n_cord = this%SurfaceMagneticCurrent(1)%probe%n_cord * precounting + 1 - CALL addSondaCoordinates (MPIDIR,this%SurfaceMagneticCurrent(1)%probe, i, j, K, node, this%SurfaceMagneticCurrent(1)%probe%n_cord) - this%SurfaceMagneticCurrent(1)%n_nml_max = this%SurfaceMagneticCurrent(1)%n_nml_max + 1 - this%SurfaceMagneticCurrent(1)%n_nml = this%SurfaceMagneticCurrent(1)%n_nml * precounting + 1 - this%SurfaceMagneticCurrent(1)%nml => addNml (nml, this%SurfaceMagneticCurrent(1)%nml, this%SurfaceMagneticCurrent(1)%n_nml) - END IF - ! - RETURN - END SUBROUTINE addCoordinates_1 - !------------------------------------------------------------------------------ - ! Selector to add coordinates - !------------------------------------------------------------------------------ - SUBROUTINE addCoordinates_2 (MPIDIR,this, i, j, K, node, precounting) - TYPE (abstractSonda), INTENT (INOUT) :: this - INTEGER (KIND=4) :: MPIDIR,i, j, K, node - INTEGER (KIND=4) :: precounting - ! - IF (ASSOCIATED(this%FarField)) THEN - this%FarField(1)%probe%n_cord_max = this%FarField(1)%probe%n_cord_max + 1 - this%FarField(1)%probe%n_cord = this%FarField(1)%probe%n_cord * precounting + 1 - CALL addSondaCoordinates (MPIDIR,this%FarField(1)%probe, i, j, K, node, this%FarField(1)%probe%n_cord) - ELSEIF (ASSOCIATED(this%Electric)) THEN - this%Electric(1)%probe%n_cord_max = this%Electric(1)%probe%n_cord_max + 1 - this%Electric(1)%probe%n_cord = this%Electric(1)%probe%n_cord * precounting + 1 - CALL addSondaCoordinates (MPIDIR,this%Electric(1)%probe, i, j, K, node, this%Electric(1)%probe%n_cord) - ELSE IF (ASSOCIATED(this%Magnetic)) THEN - this%Magnetic(1)%probe%n_cord_max = this%Magnetic(1)%probe%n_cord_max + 1 - this%Magnetic(1)%probe%n_cord = this%Magnetic(1)%probe%n_cord * precounting + 1 - CALL addSondaCoordinates (MPIDIR,this%Magnetic(1)%probe, i, j, K, node, this%Magnetic(1)%probe%n_cord) - END IF - ! - RETURN - END SUBROUTINE addCoordinates_2 - !------------------------------------------------------------------------------ - ! Method to include the coordinates introduced as parameters - ! into the Sonda object included as a parameter - ! as well. - !------------------------------------------------------------------------------ - SUBROUTINE addSondaCoordinates (mpidir,this, n_i, n_j, n_k, n_node, steps) - !-------------> Parameters - TYPE (Sonda), INTENT (INOUT) :: this - INTEGER (KIND=4), INTENT (INOUT) :: n_i, n_j, n_k, n_node,mpidir - INTEGER (KIND=4), INTENT (INOUT) :: steps - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - INTEGER (KIND=4), DIMENSION (:), POINTER :: i, j, K, node - !--> Allocate space for the new array - ALLOCATE (i(steps), j(steps), K(steps), node(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - i (n) = this%i(n) - J (n) = this%j(n) - K (n) = this%K(n) - node (n) = this%node(n) - END DO - !--> Add new value found - i (steps) = n_i - j (steps) = n_j - K (steps) = n_k - node (steps) = n_node - !--> Assigning new values - NULLIFY (this%i, this%j, this%K, this%node) - ! - this%i => i - this%j => j - this%K => K - this%node => node - !!!ROTATE MPI - IF (MPIDIR==2 ) THEN - this%i => K - this%j => I - this%K => J - ELSEIF (MPIDIR==1 ) THEN - this%i => J - this%j => K - this%K => I - ENDIF - !!!FIN - - ! - RETURN - END SUBROUTINE addSondaCoordinates - !------------------------------------------------------------------------------ - ! Method to include the nml attribute into the normal electric and magnetic - ! field probe. - !------------------------------------------------------------------------------ - FUNCTION addNml (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - INTEGER (KIND=4), INTENT (IN) :: n_Volume - INTEGER (KIND=4), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - INTEGER (KIND=4), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION addNml - !------------------------------------------------------------------------------ - ! Routine to assign the abstract class to the non-abstract Electric - !------------------------------------------------------------------------------ - FUNCTION assign_FarField (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (FarField_Sonda), INTENT (IN) :: n_Volume - TYPE (FarField_Sonda), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (FarField_Sonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION assign_FarField - !------------------------------------------------------------------------------ - ! Routine to assign the abstract class to the non-abstract Electric - !------------------------------------------------------------------------------ - FUNCTION assign_Electric (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (Electric_Sonda), INTENT (IN) :: n_Volume - TYPE (Electric_Sonda), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (Electric_Sonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION assign_Electric - !------------------------------------------------------------------------------ - ! Routine to assign the abstract class to the non-abstract Electric - !------------------------------------------------------------------------------ - FUNCTION assign_Magnetic (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (Magnetic_Sonda), INTENT (IN) :: n_Volume - TYPE (Magnetic_Sonda), DIMENSION (:), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (Magnetic_Sonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION assign_Magnetic - !------------------------------------------------------------------------------ - ! Routine to assign the abstract class to the non-abstract Electric - !------------------------------------------------------------------------------ - FUNCTION assign_NormalElectric (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (NormalElectric_Sonda), INTENT (IN) :: n_Volume - TYPE (NormalElectric_Sonda), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (NormalElectric_Sonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION assign_NormalElectric - !------------------------------------------------------------------------------ - ! Routine to assign the abstract class to the non-abstract Electric - !------------------------------------------------------------------------------ - FUNCTION assign_NormalMagnetic (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (NormalMagnetic_Sonda), INTENT (IN) :: n_Volume - TYPE (NormalMagnetic_Sonda), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (NormalMagnetic_Sonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION assign_NormalMagnetic - !------------------------------------------------------------------------------ - ! Routine to assign the abstract class to the non-abstract Electric - !------------------------------------------------------------------------------ - FUNCTION assign_SurfaceElectricCurrent (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (SurfaceElectricCurrent_Sonda), INTENT (IN) :: n_Volume - TYPE (SurfaceElectricCurrent_Sonda), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (SurfaceElectricCurrent_Sonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION assign_SurfaceElectricCurrent - !------------------------------------------------------------------------------ - ! Routine to assign the abstract class to the non-abstract Electric - !------------------------------------------------------------------------------ - FUNCTION assign_SurfaceMagneticCurrent (n_Volume, o_Volumes, steps) RESULT (Volumes) - !-------------> Parameters - INTEGER (KIND=4), INTENT (IN) :: steps - TYPE (SurfaceMagneticCurrent_Sonda), INTENT (IN) :: n_Volume - TYPE (SurfaceMagneticCurrent_Sonda), DIMENSION (steps-1), INTENT (IN) :: o_Volumes - !-------------> Used inside - INTEGER (KIND=4) :: n - !-------------> Output - TYPE (SurfaceMagneticCurrent_Sonda), DIMENSION (:), POINTER :: Volumes - !--> Allocate space for the new array - ALLOCATE (Volumes(steps)) - !--> Copy of the old array into the new one - DO n = 1, steps - 1 - Volumes (n) = o_Volumes (n) - END DO - !--> Add new value found - Volumes (steps) = n_Volume - ! - RETURN - END FUNCTION assign_SurfaceMagneticCurrent - ! - -END MODULE SondaClass diff --git a/src_main_pub/semba_fdtd.F90 b/src_main_pub/semba_fdtd.F90 index 25724eeb..0dd4de80 100755 --- a/src_main_pub/semba_fdtd.F90 +++ b/src_main_pub/semba_fdtd.F90 @@ -38,17 +38,21 @@ PROGRAM SEMBA_FDTD_launcher USE Getargs ! USE fdetypes - USE Solver + USE Solver !!!el timestepping.F90 USE Resuming !nfde parser stuff - USE NFDETypes -#ifdef CompilePrivateVersion + USE NFDETypes + use nfde_rotate_m + + +#ifdef CompilePrivateVersion USE ParseadorClass #endif - USE smbjson, only: fdtdjson_parser_t => parser_t + + USE smbjson, only: fdtdjson_parser_t => parser_t USE Preprocess_m USE storeData - + #ifdef CompileWithXDMF USE xdmf_h5 #endif @@ -60,7 +64,6 @@ PROGRAM SEMBA_FDTD_launcher use MPI_stochastic #endif #endif - ! !************************************************* !***[conformal] ****************** @@ -77,6 +80,9 @@ PROGRAM SEMBA_FDTD_launcher !************************************************* !************************************************* ! +!!! + use mtln_solver_mod, mtln_solver_t => mtln_t +!!! use interpreta_switches_m IMPLICIT NONE ! @@ -85,8 +91,7 @@ PROGRAM SEMBA_FDTD_launcher !!!241018 fin pscaling integer (KIND=IKINDMTAG) , allocatable , dimension(:,:,:) :: sggMtag integer (KIND=INTEGERSIZEOFMEDIAMATRICES) , allocatable , dimension(:,:,:) :: sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx,sggMiHy,sggMiHz - - + LOGICAL :: dummylog,finishedwithsuccess,l_auxinput, l_auxoutput, ThereArethinslots integer (KIND=4) :: myunit,jmed REAL (KIND=RKIND) :: maxSourceValue @@ -135,11 +140,16 @@ PROGRAM SEMBA_FDTD_launcher !**************************************************************************** type (entrada_t) :: l - +!!! + type (mtln_solver_t) :: mtln_solver +!!! logical :: lexis integer (kind=4) :: my_iostat -!!!!!!!!!!!!!!!!comienzo instrucciones - + + INTEGER (KIND=4) :: verdadero_mpidir + logical :: newrotate !300124 tiramos con el rotador antiguo + + newrotate=.false. !!ojo tocar luego !!200918 !!!si se lanza con -pscal se overridea esto Eps0= 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12 Mu0 = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6 @@ -343,11 +353,21 @@ PROGRAM SEMBA_FDTD_launcher !!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!! call print_credits(l) + if (trim(adjustl(l%extension))=='.nfde') then #ifdef CompilePrivateVersion - call cargaNFDE + call cargaNFDE(l%filefde,parser) #else - call cargaFDTDJSON(l%fichin, parser) + print *,'Not compiled with cargaNFDEINDEX' + stop #endif + elseif (trim(adjustl(l%extension))=='.json') then + call cargaFDTDJSON(l%fichin, parser) + else + print *, 'Neither .nfde nor .json files used as input after -i' + stop + endif + + !!!!!!!!!!!!!!!!!!!!!!! sgg%extraswitches=parser%switches !!!da preferencia a los switches por linea de comando @@ -471,12 +491,14 @@ PROGRAM SEMBA_FDTD_launcher !NOTE: md: lo necesito despues del conformal init (antes se borraba mas arriba) !REVIEW: sgg -#ifdef CompilePrivateVersion - CALL Destroy_Parser (parser) +#ifdef CompilePrivateVersion + if (trim(adjustl(l%extension))=='.nfde') then + CALL Destroy_Parser (parser) + DEALLOCATE (NFDE_FILE%lineas) + DEALLOCATE (NFDE_FILE) + nullify (NFDE_FILE) + endif #endif - DEALLOCATE (NFDE_FILE%lineas) - DEALLOCATE (NFDE_FILE) - nullify (NFDE_FILE) !!!!!!!!!!!!!!!!!!!!!!!!!!!!! #ifdef CompileWithMPI @@ -842,7 +864,8 @@ PROGRAM SEMBA_FDTD_launcher l%opcionestotales,l%sgbcfreq,l%sgbcresol,l%sgbccrank,l%sgbcdepth,l%fatalerror,l%fieldtotl,l%permitscaling, & l%EpsMuTimeScale_input_parameters, & l%stochastic,l%mpidir,l%verbose,l%precision,l%hopf,l%ficherohopf,l%niapapostprocess,l%planewavecorr, & - l%dontwritevtk,l%experimentalVideal,l%forceresampled,l%factorradius,l%factordelta,l%noconformalmapvtk ) + l%dontwritevtk,l%experimentalVideal,l%forceresampled,l%factorradius,l%factordelta,l%noconformalmapvtk, & + mtln_solver) deallocate (sggMiEx, sggMiEy, sggMiEz,sggMiHx, sggMiHy, sggMiHz,sggMiNo,sggMtag) else @@ -1047,17 +1070,19 @@ PROGRAM SEMBA_FDTD_launcher #ifdef CompilePrivateVersion -subroutine cargaNFDE +subroutine cargaNFDE(local_nfde,local_parser) + CHARACTER (LEN=BUFSIZE) :: local_nfde + TYPE (Parseador), POINTER :: local_parser INTEGER (KIND=8) :: numero,i8,troncho,longitud integer (kind=4) :: mpi_t_linea_t,longitud4 IF (l%existeNFDE) THEN - WRITE (dubuf,*) 'INIT Reading file '//trim (adjustl(whoami))//' ', trim (adjustl(l%fileFDE)) + WRITE (dubuf,*) 'INIT Reading file '//trim (adjustl(whoami))//' ', trim (adjustl(local_nfde)) CALL print11 (l%layoutnumber, dubuf) !!!!!!!!!!!!!!!!!!!!!!! #ifdef CompileWithMPI CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) if (l%layoutnumber==0) then - NFDE_FILE => cargar_NFDE_FILE (l%fileFDE) + NFDE_FILE => cargar_NFDE_FILE (local_nfde) !!!ya se allocatea dentro else ALLOCATE (NFDE_FILE) @@ -1114,21 +1139,39 @@ subroutine cargaNFDE !!!close (6729) !!!!!! #else - NFDE_FILE => cargar_NFDE_FILE (l%fileFDE) + NFDE_FILE => cargar_NFDE_FILE (local_nfde) #endif write(dubuf,*) '[OK]'; call print11(l%layoutnumber,dubuf) !---> - END IF + END IF NFDE_FILE%mpidir=l%mpidir !!!!!!!!!!!!!!!!!!! - WRITE (dubuf,*) 'INIT interpreting geometrical data from ', trim (adjustl(l%fileFDE)) + WRITE (dubuf,*) 'INIT interpreting geometrical data from ', trim (adjustl(local_nfde)) CALL print11 (l%layoutnumber, dubuf) !!!!!!!!!! - parser => newparser (NFDE_FILE) + if(newrotate) then + verdadero_mpidir=NFDE_FILE%mpidir + NFDE_FILE%mpidir=3 !no lo rota el parseador antiguo + endif + local_parser => newparser (NFDE_FILE) +#ifdef CompileWithMPI + CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) +#endif + if(newrotate) then + NFDE_FILE%mpidir=verdadero_mpidir !restorealo + call nfde_rotate (local_parser,NFDE_FILE%mpidir) !lo rota el parseador nuevo +#ifdef CompileWithMPI + CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) +#endif + endif l%thereare_stoch=NFDE_FILE%thereare_stoch l%mpidir=NFDE_FILE%mpidir !bug 100419 -!!!!!!!!!!! - write(dubuf,*) '[OK]'; call print11(l%layoutnumber,dubuf) +!!!!!!!!!!! + ! write(dubuf,*) '[OK]'; call print11(l%layoutnumber,dubuf) + write(dubuf,*) '[OK] '//trim(adjustl(whoami))//' newparser (NFDE_FILE)'; call print11(0,dubuf) +#ifdef CompileWithMPI + CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) +#endif return end subroutine cargaNFDE @@ -1161,8 +1204,7 @@ subroutine NFDE2sgg CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) #endif CALL read_limits_nogeom (l%layoutnumber,l%size, sgg, fullsize, SINPML_fullsize, parser,l%MurAfterPML,l%mur_exist) - - + dtantesdecorregir=sgg%dt !!!!!corrige el delta de t si es necesario !sgg15 310715 bug distintos sgg%dt !!!!!!!!!! @@ -1271,6 +1313,9 @@ subroutine NFDE2sgg CALL read_geomData (sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx,sggMiHy,sggMiHz, l%fichin, l%layoutnumber, l%size, SINPML_fullsize, fullsize, parser, & l%groundwires,l%attfactorc,l%mibc,l%sgbc,l%sgbcDispersive,l%MEDIOEXTRA,maxSourceValue,l%skindepthpre,l%createmapvtk,l%input_conformal_flag,l%CLIPREGION,l%boundwireradius,l%maxwireradius,l%updateshared,l%run_with_dmma, & eps0,mu0,.false.,l%hay_slanted_wires,l%verbose,l%ignoresamplingerrors,tagtype,l%wiresflavor) +!!!!mtln constructor 100424 + mtln_solver = mtlnCtor(parser%mtln) +!!!! WRITE (dubuf,*) '[OK] ENDED NFDE --------> GEOM' CALL print11 (l%layoutnumber, dubuf) !writing @@ -1295,7 +1340,7 @@ subroutine NFDE2sgg endif #endif #endif - ELSE !del l%size==1 + ELSE !del l%size==1 #ifdef CompileWithMPI CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) #ifdef CompileWithStochastic @@ -1303,15 +1348,17 @@ subroutine NFDE2sgg call HalvesStochasticMPI(l%layoutnumber,l%size,l%simu_devia) endif #endif + + CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) !!!ahora divide el espacio computacional CALL MPIdivide (sgg, fullsize, SINPML_fullsize, l%layoutnumber, l%size, l%forcing, l%forced, l%slicesoriginales, l%resume,l%fatalerror) ! - CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) + CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) if (l%fatalerror) then !intenta recuperarte return endif - + ! if the layout is pure PML then take at least a line of non PML to build the PML data insider read_geomDAta ! Uses extra memory but later matrix sggm is deallocated in favor of smaller sggMIEX, etc DO field = iEx, iHz @@ -1320,8 +1367,8 @@ subroutine NFDE2sgg sgg%Alloc(field)%ZE = Max (sgg%Alloc(field)%ZE, SINPML_fullsize(field)%ZI+1) sgg%Alloc(field)%ZI = Min (sgg%Alloc(field)%ZI, SINPML_fullsize(field)%ZE-1) END DO - ! - CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) + ! + CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) !!incluido aqui pq se precisa para clip 16/07/15 DO field = iEx, iHz sgg%SINPMLSweep(field)%XI = Max (SINPML_fullsize(field)%XI, sgg%Sweep(field)%XI) @@ -1333,10 +1380,13 @@ subroutine NFDE2sgg END DO !!fin 16/07/15 WRITE (dubuf,*) 'INIT NFDE --------> GEOM' - CALL print11 (l%layoutnumber, dubuf) + CALL print11 (l%layoutnumber, dubuf) + CALL read_geomData (sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx,sggMiHy,sggMiHz, l%fichin, l%layoutnumber, l%size, SINPML_fullsize, fullsize, parser, & l%groundwires,l%attfactorc,l%mibc,l%sgbc,l%sgbcDispersive,l%MEDIOEXTRA,maxSourceValue,l%skindepthpre,l%createmapvtk,l%input_conformal_flag,l%CLIPREGION,l%boundwireradius,l%maxwireradius,l%updateshared,l%run_with_dmma, & eps0,mu0,l%simu_devia,l%hay_slanted_wires,l%verbose,l%ignoresamplingerrors,tagtype,l%wiresflavor) + + #ifdef CompileWithMPI !wait until everything comes out CALL MPI_Barrier (SUBCOMM_MPI, l%ierr) diff --git a/src_main_pub/timestepping.F90 b/src_main_pub/timestepping.F90 index 5070d134..b24074a4 100755 --- a/src_main_pub/timestepping.F90 +++ b/src_main_pub/timestepping.F90 @@ -91,7 +91,7 @@ module Solver use HollandWires #endif #ifdef CompileWithWires_mtln - use HollandWires_mtln + use Wire_bundles_mtln_mod #endif #ifdef CompileWithBerengerWires use WiresBerenger @@ -113,7 +113,9 @@ module Solver USE CALC_CONSTANTS #ifdef CompileWithPrescale USE P_rescale -#endif +#endif + use mtln_solver_mod, mtln_solver_t => mtln_t + use Wire_bundles_mtln_mod !! #ifdef CompileWithProfiling use nvtx @@ -140,8 +142,12 @@ subroutine launch_simulation(sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx opcionestotales,sgbcFreq,sgbcresol,sgbccrank,sgbcDepth,fatalerror,fieldtotl,permitscaling, & EpsMuTimeScale_input_parameters, & stochastic,mpidir,verbose,precision,hopf,ficherohopf,niapapostprocess,planewavecorr, & - dontwritevtk,experimentalVideal,forceresampled,factorradius,factordelta,noconformalmapvtk ) - + dontwritevtk,experimentalVideal,forceresampled,factorradius,factordelta,noconformalmapvtk, & + mtln_solver) + +!!! + type (mtln_solver_t) :: mtln_solver +!!! logical :: noconformalmapvtk logical :: hopf,experimentalVideal,forceresampled character (LEN=BUFSIZE) :: ficherohopf @@ -705,11 +711,6 @@ subroutine launch_simulation(sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx inductance_model,wirethickness,groundwires,strictOLD,TAPARRABOS,g2,wiresflavor,SINPML_fullsize,fullsize,wirecrank,dtcritico,eps0,mu0,simu_devia,stochastic,verbose,factorradius,factordelta) l_auxinput=thereare%Wires l_auxoutput=l_auxinput -!! -#ifdef CompileWithWires_mtln - call InitWires_mtln() -#endif -!! #ifdef CompileWithMPI call MPI_Barrier(SUBCOMM_MPI,ierr) call MPI_AllReduce( l_auxinput, l_auxoutput, 1_4, MPI_LOGICAL, MPI_LOR, MPI_COMM_WORLD, ierr) @@ -808,6 +809,13 @@ subroutine launch_simulation(sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx call MPI_Barrier(SUBCOMM_MPI,ierr) #endif + +#ifdef CompileWithWires_mtln + call InitWires_mtln(sgg,Ex,Ey,Ez,mtln_solver,thereAre%MTLNbundles) +#endif + + + #ifdef CompileWithAnisotropic !Anisotropic #ifdef CompileWithMPI @@ -1341,10 +1349,7 @@ subroutine launch_simulation(sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx if (wirecrank) then call AdvanceWiresEcrank(sgg,n, layoutnumber,wiresflavor,simu_devia,stochastic) else - call AdvanceWiresE(sgg,n, layoutnumber,wiresflavor,simu_devia,stochastic,experimentalVideal,wirethickness,eps0,mu0) -#ifdef CompileWithWires_mtln - call AdvanceWiresE_mtln() -#endif + call AdvanceWiresE(sgg,n, layoutnumber,wiresflavor,simu_devia,stochastic,experimentalVideal,wirethickness,eps0,mu0) endif endif endif @@ -1361,6 +1366,9 @@ subroutine launch_simulation(sgg,sggMtag,sggMiNo,sggMiEx,sggMiEy,sggMiEz,sggMiHx call AdvanceWiresE_Slanted(sgg,n) endif #endif +#ifdef CompileWithWires_mtln + if (thereAre%MTLNbundles) call AdvanceWiresE_mtln(sgg, Idxe,Idye,Idze,Idxh,Idyh,Idzh,eps0,mu0,mtln_solver) +#endif If (Thereare%PMLbodies) then !waveport absorbers call AdvancePMLbodyE endif diff --git a/src_mtln/CMakeLists.txt b/src_mtln/CMakeLists.txt index 8f8b1431..00de2052 100644 --- a/src_mtln/CMakeLists.txt +++ b/src_mtln/CMakeLists.txt @@ -27,9 +27,11 @@ add_library(ngspice_interface "ngspice_interface.c" ) -target_link_libraries(ngspice_interface PRIVATE - ${NGSPICE_LIB} -) +target_link_libraries(ngspice_interface PRIVATE ${NGSPICE_LIB}) + +if (${CMAKE_SYSTEM_NAME} STREQUAL "Linux") + target_link_libraries(ngspice_interface PRIVATE stdc++) +endif () target_link_libraries(mtlnsolver PRIVATE ${LAPACK_LIB} diff --git a/src_mtln/mtl.F90 b/src_mtln/mtl.F90 index e1daa8f1..c4c461d5 100644 --- a/src_mtln/mtl.F90 +++ b/src_mtln/mtl.F90 @@ -3,7 +3,7 @@ module mtl_mod ! use NFDETypes use utils_mod use dispersive_mod - use mtln_types_mod, only: segment_relative_position_t + use mtln_types_mod, only: external_field_segment_t implicit none @@ -20,7 +20,7 @@ module mtl_mod integer :: conductor_in_parent type(transfer_impedance_per_meter_t) :: transfer_impedance - type(segment_relative_position_t), allocatable, dimension(:) :: segment_relative_positions + type(external_field_segment_t), allocatable, dimension(:) :: external_field_segments contains procedure :: setTimeStep @@ -84,7 +84,7 @@ function mtlHomogeneous(lpul, cpul, rpul, gpul, & step_size, name, & dt, parent_name, conductor_in_parent, & transfer_impedance, & - segment_relative_positions) result(res) + external_field_segments) result(res) type(mtl_t) :: res real, intent(in), dimension(:,:) :: lpul, cpul, rpul, gpul real, intent(in), dimension(:) :: step_size @@ -95,7 +95,7 @@ function mtlHomogeneous(lpul, cpul, rpul, gpul, & character(len=*), intent(in), optional :: parent_name integer, intent(in), optional :: conductor_in_parent type(transfer_impedance_per_meter_t), intent(in), optional :: transfer_impedance - type(segment_relative_position_t), intent(in), dimension(:), optional :: segment_relative_positions + type(external_field_segment_t), intent(in), dimension(:), optional :: external_field_segments res%name = name res%step_size = step_size @@ -130,15 +130,16 @@ function mtlHomogeneous(lpul, cpul, rpul, gpul, & res%transfer_impedance = transfer_impedance end if - if (present(segment_relative_positions)) then - res%segment_relative_positions = segment_relative_positions + if (present(external_field_segments)) then + res%external_field_segments = external_field_segments end if end function function mtlInhomogeneous(lpul, cpul, rpul, gpul, & step_size, name, & dt, parent_name, conductor_in_parent, & - transfer_impedance) result(res) + transfer_impedance, & + external_field_segments) result(res) type(mtl_t) :: res real, intent(in), dimension(:,:,:) :: lpul, cpul, rpul, gpul real, intent(in), dimension(:) :: step_size @@ -148,6 +149,7 @@ function mtlInhomogeneous(lpul, cpul, rpul, gpul, & character(len=*), intent(in), optional :: parent_name integer, intent(in), optional :: conductor_in_parent type(transfer_impedance_per_meter_t), intent(in), optional :: transfer_impedance + type(external_field_segment_t), intent(in), dimension(:), optional :: external_field_segments res%name = name res%step_size = step_size @@ -181,6 +183,10 @@ function mtlInhomogeneous(lpul, cpul, rpul, gpul, & res%transfer_impedance = transfer_impedance end if + if (present(external_field_segments)) then + res%external_field_segments = external_field_segments + end if + end function subroutine initDirections(this) diff --git a/src_mtln/mtl_bundle.F90 b/src_mtln/mtl_bundle.F90 index 93b12a7b..d23f7bb9 100644 --- a/src_mtln/mtl_bundle.F90 +++ b/src_mtln/mtl_bundle.F90 @@ -21,7 +21,7 @@ module mtl_bundle_mod real, dimension(:,:,:), allocatable :: v_term, i_term real, dimension(:,:,:), allocatable :: v_diff, i_diff - type(segment_relative_position_t), dimension(:), allocatable :: segment_relative_positions + type(external_field_segment_t), dimension(:), allocatable :: external_field_segments contains @@ -37,7 +37,6 @@ module mtl_bundle_mod procedure :: advanceCurrent => bundle_advanceCurrent procedure :: addTransferImpedance => bundle_addTransferImpedance ! procedure :: setConnectorTransferImpedance - procedure :: setExternalCurrent => bundle_setExternalCurrent procedure :: setExternalVoltage => bundle_setExternalVoltage procedure :: updateExternalCurrent => bundle_updateExternalCurrent @@ -64,7 +63,7 @@ function mtldCtor(levels, name) result(res) res%dt = levels(1)%lines(1)%dt res%step_size = levels(1)%lines(1)%step_size res%number_of_divisions = size(res%step_size,1) - res%segment_relative_positions = levels(1)%lines(1)%segment_relative_positions + res%external_field_segments = levels(1)%lines(1)%external_field_segments call res%initialAllocation() call res%mergePULMatrices(levels) call res%mergeDispersiveMatrices(levels) @@ -309,22 +308,25 @@ subroutine bundle_advanceCurrent(this) ! call this%transfer_impedance%updatePhi(i_prev, i_now) end subroutine - subroutine bundle_setExternalCurrent(this, current) + subroutine bundle_setExternalVoltage(this) class(mtl_bundle_t) :: this - real, dimension(:), intent(in) :: current - this%i(1,:) = current(:) - end subroutine - - subroutine bundle_setExternalVoltage(this, voltage) - class(mtl_bundle_t) :: this - real, dimension(:), intent(in) :: voltage - this%v(1,:) = voltage(:) + integer :: i + do i = 1, size(this%v,2) + this%v(1, i) = this%external_field_segments(i)%Efield_main2wire * this%step_size(i) + end do end subroutine subroutine bundle_updateExternalCurrent(this, current) class(mtl_bundle_t) :: this - real, dimension(:), intent(inout) :: current - current(:) = this%i(1,:) + real, dimension(:,:,:), intent(inout) :: current + ! real, dimension(:), intent(inout) :: current + integer, dimension(:), allocatable :: position + integer :: i + do i = 1, size(this%i,2) + position = this%external_field_segments(i)%position + current(position(1), position(2), position(3)) = this%i(1,i) + end do + ! current(:) = this%i(1,:) end subroutine end module mtl_bundle_mod \ No newline at end of file diff --git a/src_mtln/mtln_solver.F90 b/src_mtln/mtln_solver.F90 index 9e016912..391d8846 100644 --- a/src_mtln/mtln_solver.F90 +++ b/src_mtln/mtln_solver.F90 @@ -47,6 +47,11 @@ function mtlnCtor(parsed) result(res) type(preprocess_t) :: pre pre = preprocess(parsed) + if (size(pre%bundles) == 0) then + res%number_of_bundles = 0 + return + end if + res%dt = pre%dt res%time = 0.0 @@ -70,13 +75,11 @@ subroutine initNodes(this) end do end subroutine - subroutine mtln_step(this, currents, voltages) + subroutine mtln_step(this) class(mtln_t) :: this - real, dimension(:,:), intent(in out) :: currents - real, dimension(:,:), intent(in) :: voltages integer :: i - call this%setExternalVoltage(voltages) + call this%setExternalVoltage() call this%advanceBundlesVoltage() call this%advanceNWVoltage() @@ -85,7 +88,7 @@ subroutine mtln_step(this, currents, voltages) call this%advanceTime() call this%updateProbes() - call this%updateExternalCurrent(currents) + ! call this%updateExternalCurrent() end subroutine subroutine step_alone(this) @@ -102,22 +105,23 @@ subroutine step_alone(this) end subroutine - subroutine setExternalVoltage(this, voltages) + subroutine setExternalVoltage(this) class(mtln_t) :: this - real, dimension(:,:), intent(in) :: voltages integer :: i do i = 1, this%number_of_bundles - call this%bundles(i)%setExternalVoltage(voltages(i,:)) + call this%bundles(i)%setExternalVoltage() end do end subroutine subroutine updateExternalCurrent(this, currents) class(mtln_t) :: this - real, dimension(:,:), intent(inout) :: currents + ! real, dimension(:,:), intent(inout) :: currents + real, dimension(:,:,:), intent(inout) :: currents integer :: i do i = 1, this%number_of_bundles - call this%bundles(i)%updateExternalCurrent(currents(i,:)) + call this%bundles(i)%updateExternalCurrent(currents) + ! call this%bundles(i)%updateExternalCurrent(currents(i,:)) end do diff --git a/src_mtln/mtln_types.F90 b/src_mtln/mtln_types.F90 index 41ec3bc5..4d515258 100644 --- a/src_mtln/mtln_types.F90 +++ b/src_mtln/mtln_types.F90 @@ -1,5 +1,5 @@ module mtln_types_mod - + use fdetypes, ONLY: RKIND !sggmtln implicit none integer, parameter :: TERMINATION_UNDEFINED = -1 @@ -21,12 +21,21 @@ module mtln_types_mod integer, parameter :: PROBE_TYPE_VOLTAGE = 1 integer, parameter :: PROBE_TYPE_CURRENT = 2 - type :: segment_relative_position_t + integer, parameter :: DIRECTION_X_POS = 1 + integer, parameter :: DIRECTION_X_NEG = -1 + integer, parameter :: DIRECTION_Y_POS = 2 + integer, parameter :: DIRECTION_Y_NEG = -2 + integer, parameter :: DIRECTION_Z_POS = 3 + integer, parameter :: DIRECTION_Z_NEG = -3 + + type :: external_field_segment_t integer, dimension(3) ::position + integer :: direction + real (kind=rkind) , pointer :: Efield_wire2main, Efield_main2wire contains private - procedure :: segment_relative_positions_eq - generic, public :: operator(==) => segment_relative_positions_eq + procedure :: external_field_segments_eq + generic, public :: operator(==) => external_field_segments_eq end type type :: termination_t @@ -111,7 +120,7 @@ module mtln_types_mod integer :: conductor_in_parent = -1 type(connector_t), pointer :: initial_connector => null() type(connector_t), pointer :: end_connector => null() - type(segment_relative_position_t), allocatable, dimension(:) :: segment_relative_positions + type(external_field_segment_t), allocatable, dimension(:) :: external_field_segments contains private @@ -175,7 +184,7 @@ recursive elemental logical function cable_eq(a,b) cable_eq = cable_eq .and. all(a%step_size == b%step_size) cable_eq = cable_eq .and. (a%transfer_impedance == b%transfer_impedance) cable_eq = cable_eq .and. (a%conductor_in_parent == b%conductor_in_parent) - cable_eq = cable_eq .and. all(a%segment_relative_positions == b%segment_relative_positions) + cable_eq = cable_eq .and. all(a%external_field_segments == b%external_field_segments) if (.not. cable_eq) then @@ -298,10 +307,31 @@ elemental logical function terminal_network_eq(a,b) all(a%connections == b%connections) end function - elemental logical function segment_relative_positions_eq(a,b) - class(segment_relative_position_t), intent(in) :: a,b - segment_relative_positions_eq = & - all(a%position == b%position) + elemental logical function external_field_segments_eq(a,b) + class(external_field_segment_t), intent(in) :: a,b + external_field_segments_eq = & + all(a%position == b%position) .and. & + a%direction == b%direction + + if (.not. associated(a%Efield_main2wire) .and. .not. associated(b%Efield_main2wire)) then + external_field_segments_eq = external_field_segments_eq .and. .true. + else if ((associated(a%Efield_main2wire) .and. .not. associated(b%Efield_main2wire)) .or. & + (.not. associated(a%Efield_main2wire) .and. associated(b%Efield_main2wire))) then + external_field_segments_eq = external_field_segments_eq .and. .false. + else + external_field_segments_eq = external_field_segments_eq .and. (a%Efield_main2wire == b%Efield_main2wire) + end if + + if (.not. associated(a%Efield_wire2main) .and. .not. associated(b%Efield_wire2main)) then + external_field_segments_eq = external_field_segments_eq .and. .true. + else if ((associated(a%Efield_wire2main) .and. .not. associated(b%Efield_wire2main)) .or. & + (.not. associated(a%Efield_wire2main) .and. associated(b%Efield_wire2main))) then + external_field_segments_eq = external_field_segments_eq .and. .false. + else + external_field_segments_eq = external_field_segments_eq .and. (a%Efield_wire2main == b%Efield_wire2main) + end if + + end function subroutine terminal_connection_add_node(this, node) @@ -313,14 +343,14 @@ subroutine terminal_connection_add_node(this, node) subroutine terminal_network_add_connection(this, connection) class(terminal_network_t) :: this - type(terminal_connection_t) :: connection - type(terminal_connection_t), dimension(:), allocatable :: newConnections - integer :: newConnectionsSize - if (.not. allocated(this%connections)) allocate(this%connections(0)) - - allocate(newConnections( size(this%connections) + 1 ) ) + type(terminal_connection_t) :: connection + type(terminal_connection_t), dimension(:), allocatable :: newConnections + integer :: newConnectionsSize + if (.not. allocated(this%connections)) allocate(this%connections(0)) + + allocate(newConnections( size(this%connections) + 1 ) ) newConnectionsSize = size(newConnections) - newConnections(1:newConnectionsSize-1) = this%connections + newConnections(1:newConnectionsSize-1) = this%connections newConnections(newConnectionsSize) = connection call MOVE_ALLOC(from=newConnections, to=this%connections) end subroutine diff --git a/src_mtln/ngspice_interface.c b/src_mtln/ngspice_interface.c index 552572da..d4df02a6 100644 --- a/src_mtln/ngspice_interface.c +++ b/src_mtln/ngspice_interface.c @@ -1,3 +1,4 @@ + #include #include #include diff --git a/src_mtln/preprocess.F90 b/src_mtln/preprocess.F90 index c7b5e3dc..b84852af 100644 --- a/src_mtln/preprocess.F90 +++ b/src_mtln/preprocess.F90 @@ -211,7 +211,7 @@ function buildLineFromCable(cable) result(res) parent_name = parent_name, & conductor_in_parent = conductor_in_parent, & transfer_impedance = cable%transfer_impedance, & - segment_relative_positions = cable%segment_relative_positions) + external_field_segments = cable%external_field_segments) if (associated(cable%initial_connector)) call addConnector(res, cable%initial_connector, 0) if (associated(cable%end_connector)) call addConnector(res, cable%initial_connector, size(res%rpul,1)) diff --git a/src_pyWrapper/pyWrapper.py b/src_pyWrapper/pyWrapper.py index db495e8f..b0c8e432 100644 --- a/src_pyWrapper/pyWrapper.py +++ b/src_pyWrapper/pyWrapper.py @@ -2,23 +2,82 @@ import json import os import glob, re +import pandas as pd +import numpy as np + +def positionStrToCell(pos_str): + pos = pos_str.split('_') + return np.array([int(pos[0]), int(pos[1]), int(pos[2])]) + +class Probe(): + + def __init__(self, probe_filename): + self.filename = probe_filename + + current_probe_tags = ['_Wx_', '_Wy_', '_Wz_'] + far_field_tag = ['_FF_'] + movie_tags = ['_ExC_', '_EyC_', '_EzC_', '_HxC_', '_HyC_', '_HzC_'] + all_tags = current_probe_tags + far_field_tag + movie_tags + + + basename = os.path.basename(self.filename) + self.case_name, basename_with_no_case_name = basename.split('.fdtd_') + basename_with_no_case_name = os.path.splitext(basename_with_no_case_name)[0] + + + for tag in all_tags: + ids = [m.start() for m in re.finditer(tag, basename_with_no_case_name)] + if len(ids) == 0: + continue + elif len(ids) == 1: + if tag in current_probe_tags: + self.type = 'wire' + self.name, position_str = basename_with_no_case_name.split(tag) + self.cell = positionStrToCell(position_str) + self.segment_tag = int(position_str.split('_s')[1]) + self.df = pd.read_csv(self.filename, sep='\s+') + self.df = self.df.rename(columns={'t': 'time', basename: 'current'}) + elif tag in far_field_tag: + self.type = 'farField' + self.name, positions_str = basename_with_no_case_name.split(tag) + init_str, end_str = pos = positions_str.split('__') + self.cell_init = positionStrToCell(init_str) + self.cell_end = positionStrToCell(end_str) + self.df = pd.read_csv(self.filename, sep='\s+') + elif tag in movie_tags: + self.type = 'movie' + self.name, positions_str = basename_with_no_case_name.split(tag) + init_str, end_str = pos = positions_str.split('__') + self.cell_init = positionStrToCell(init_str) + self.cell_end = positionStrToCell(end_str) + else: + raise ValueError("Unable to determine probe name or type for a probe with name:" + basename) + try: + self.type + self.name + except: + raise ValueError('Unable to determine type for probe' + basename) + + + def __getitem__(self, key): + return self.df[key] class FDTD(): def __init__(self, input_filename, path_to_exe): - self.file_name = input_filename + self.filename = input_filename self.path_to_exe = path_to_exe - self.folder = os.path.dirname(self.file_name) - self.case = os.path.basename(self.file_name).split('.json')[0] + self.folder = os.path.dirname(self.filename) + self.case = os.path.basename(self.filename).split('.json')[0] self.hasRun = False def run(self): os.chdir(self.folder) - self.output = subprocess.run([self.path_to_exe, "-i",self.file_name]) + self.output = subprocess.run([self.path_to_exe, "-i",self.filename]) self.hasRun = True def readJsonDict(self): - with open(self.file_name) as input_file: + with open(self.filename) as input_file: return json.load(input_file) def getSolvedProbeFilenames(self, probe_name): @@ -26,9 +85,13 @@ def getSolvedProbeFilenames(self, probe_name): if not "probes" in input_json: raise ValueError('Solver does not contain probes.') - for probe in input_json["probes"]: - probeFiles = [x for x in glob.glob('*dat') if re.match(self.case + '_' + probe_name + '.*dat',x)] - return probeFiles + file_extensions = ('*.dat', '*.bin') + probeFiles = [] + for ext in file_extensions: + newProbes = [x for x in glob.glob(ext) if re.match(self.case + '_' + probe_name, x)] + probeFiles.extend(newProbes) + + return probeFiles def hasFinishedSuccessfully(self): if self.hasRun: diff --git a/src_wires_pub/wires_mtln.F90 b/src_wires_pub/wires_mtln.F90 index e29423d5..79061094 100644 --- a/src_wires_pub/wires_mtln.F90 +++ b/src_wires_pub/wires_mtln.F90 @@ -3,45 +3,112 @@ ! Module thin wires from Wires paper !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -module HollandWires_mtln +module Wire_bundles_mtln_mod use report use fdetypes - use wiresHolland_constants -!!! - use mtln_solver_mod -!!! - implicit none + use mtln_solver_mod , mtln_solver_t => mtln_t - type(Thinwires_t) , target ,save :: HWires - - private - public InitWires_mtln,AdvanceWiresE_mtln + implicit none + real (kind=rkind_wires) :: eps0, mu0 + private + public InitWires_mtln, AdvanceWiresE_mtln + contains - subroutine InitWires_mtln() - - return - endsubroutine InitWires_mtln - - subroutine AdvanceWiresE_mtln() - class(mtln_t), pointer :: mtln_wir - real, dimension(:,:), allocatable :: currents - real, dimension(:,:), allocatable :: voltages -! do n=1,HWires%NumCurrentSegments - !segmento => Hwires%segment(n) - !segmento%efield_wire2main = real(segmento%efield_wire2main,kind=rkind_wires) - segmento%cte5 * segmento%current -!!!... - !Segmento%Current = Segmento%Current + Segmento%cte2*real(Segmento%Efield_main2wire,KIND=RKIND_wires) -! end do - - !!!call mtln_wir%mtln_step(currents, voltages) -! call mtln_step(mtln_wir, currents, voltages) - - return + + subroutine InitWires_mtln(sgg,Ex,Ey,Ez,mtln_solver,thereAreMTLNbundles) + type (SGGFDTDINFO), intent(IN), target :: sgg + REAL (KIND=RKIND), intent(inout), target :: & + Ex(sgg%Alloc(iEx)%XI : sgg%Alloc(iEx)%XE, & + sgg%Alloc(iEx)%YI : sgg%Alloc(iEx)%YE, & + sgg%Alloc(iEx)%ZI : sgg%Alloc(iEx)%ZE), & + Ey(sgg%Alloc(iEy)%XI : sgg%Alloc(iEy)%XE, & + sgg%Alloc(iEy)%YI : sgg%Alloc(iEy)%YE, & + sgg%Alloc(iEy)%ZI : sgg%Alloc(iEy)%ZE), & + Ez(sgg%Alloc(iEz)%XI : sgg%Alloc(iEz)%XE, & + sgg%Alloc(iEz)%YI : sgg%Alloc(iEz)%YE, & + sgg%Alloc(iEz)%ZI : sgg%Alloc(iEz)%ZE) + class(mtln_solver_t) :: mtln_solver + logical :: thereAreMTLNbundles + + if (mtln_solver%number_of_bundles>=1) then + thereAreMTLNbundles=.true. + else + thereAreMTLNbundles=.false. + return + endif + + block + integer (kind=4) :: i, j, k, m, n + do m = 1, mtln_solver%NUMBER_OF_BUNDLES + do n = 1, ubound(mtln_solver%bundles(m)%external_field_segments,1) + call readGridIndices(i, j, k, mtln_solver%bundles(m)%external_field_segments(n)) + select case (abs(mtln_solver%bundles(m)%external_field_segments(n)%direction)) + case(1) + mtln_solver%bundles(m)%external_field_segments(n)%Efield_wire2main => Ex(i, j, k) + mtln_solver%bundles(m)%external_field_segments(n)%Efield_main2wire => Ex(i, j, k) + case(2) + mtln_solver%bundles(m)%external_field_segments(n)%Efield_wire2main => Ey(i, j, k) + mtln_solver%bundles(m)%external_field_segments(n)%Efield_main2wire => Ey(i, j, k) + case(3) + mtln_solver%bundles(m)%external_field_segments(n)%Efield_wire2main => Ez(i, j, k) + mtln_solver%bundles(m)%external_field_segments(n)%Efield_main2wire => Ez(i, j, k) + end select + + end do + end do + end block + + end subroutine InitWires_mtln + + subroutine AdvanceWiresE_mtln(sgg,Idxh, Idyh, Idzh, eps00,mu00,mtln_solver) + class(mtln_solver_t) :: mtln_solver + type (SGGFDTDINFO), intent(IN), target :: sgg + real (kind=RKIND), dimension (:), intent(in) :: & + Idxh(sgg%ALLOC(iEx)%XI : sgg%ALLOC(iEx)%XE),& + Idyh(sgg%ALLOC(iEy)%YI : sgg%ALLOC(iEy)%YE),& + Idzh(sgg%ALLOC(iEz)%ZI : sgg%ALLOC(iEz)%ZE) + real(KIND=RKIND) :: cte, eps00, mu00 + integer (kind=4) :: i, j, k, direction, m, n + real (kind= rkind) :: current + real(kind = rkind), pointer :: punt + + eps0 = eps00 + mu0 = mu00 + + call mtln_solver%step() + + do m = 1, mtln_solver%number_of_bundles + do n = 1, ubound(mtln_solver%bundles(m)%external_field_segments,1) + call readGridIndices(i, j, k, mtln_solver%bundles(m)%external_field_segments(n)) + direction = mtln_solver%bundles(m)%external_field_segments(n)%direction + punt => mtln_solver%bundles(m)%external_field_segments(n)%Efield_wire2main + current = mtln_solver%bundles(m)%i(1, n) + select case (abs(direction)) + case(1) + cte = sgg%dt /eps0 *Idyh(j)*idzh(k) + case(2) + cte = sgg%dt /eps0 *Idzh(k)*idxh(i) + case(3) + cte = sgg%dt /eps0 *Idxh(i)*idyh(j) + end select + punt = punt - sign(cte,real(direction,kind=rkind)) * current + end do + end do end subroutine AdvanceWiresE_mtln -end module HollandWires_mtln + + subroutine readGridIndices(i, j, k, field_segment) + type(external_field_segment_t), intent(in) :: field_segment + integer, intent(inout) :: i, j, k + i = field_segment%position(1) + j = field_segment%position(2) + k = field_segment%position(3) + end subroutine + + +end module Wire_bundles_mtln_mod diff --git a/test/mtln/mtln_tests.h b/test/mtln/mtln_tests.h index 81bc0e4c..db752443 100644 --- a/test/mtln/mtln_tests.h +++ b/test/mtln/mtln_tests.h @@ -43,20 +43,26 @@ TEST(mtln, mtl_inhomogeneous) { EXPECT_EQ(0, test_mtl_init_inhomogeneous()); } TEST(mtln, mtl_time_step) { EXPECT_EQ(0, test_mtl_time_step()); } TEST(mtln, mtl_wrong_dt) { EXPECT_EQ(0, test_mtl_wrong_dt()); } TEST(mtln, mtl_bundle_init) { EXPECT_EQ(0, test_mtl_bundle_init()); } + TEST(mtln, mtln_types) {EXPECT_EQ(0, test_mtln_types()); } + TEST(mtln, fhash_arrays) { EXPECT_EQ(0, test_fhash_arrays()); } TEST(mtln, fhash_cables) { EXPECT_EQ(0, test_fhash_cables()); } TEST(mtln, fhash) { EXPECT_EQ(0, test_fhash()); } + TEST(mtln, preprocess_conductors_before_cable) { EXPECT_EQ(0, test_preprocess_conductors_before_cable()); } TEST(mtln, preprocess_conductors_in_level) { EXPECT_EQ(0, test_preprocess_conductors_in_level()); } TEST(mtln, preprocess_zt_conductor_ranges_2) { EXPECT_EQ(0, test_preprocess_zt_conductor_ranges_2()); } TEST(mtln, preprocess_zt_conductor_ranges) { EXPECT_EQ(0, test_preprocess_zt_conductor_ranges()); } + TEST(mtln, math_eigvals) { EXPECT_EQ(0, test_math_eigvals()); } TEST(mtln, math_matmul_broadcast) { EXPECT_EQ(0, test_math_matmul_broadcast()); } + TEST(mtln, dispersive_init_1_pole) { EXPECT_EQ(0, test_dispersive_init_1_pole()); } TEST(mtln, dispersive_init_2_poles) { EXPECT_EQ(0, test_dispersive_init_2_poles()); } TEST(mtln, dispersive_init_1_pole_3_levels) { EXPECT_EQ(0, test_dispersive_init_1_pole_3_levels()); } TEST(mtln, dispersive_init_1_pole_lines_with_lumped) { EXPECT_EQ(0, test_dispersive_init_1_pole_lines_with_lumped()); } + TEST(mtln, spice_tran) { EXPECT_EQ(0, test_spice_tran()); } TEST(mtln, spice_tran_2) { EXPECT_EQ(0, test_spice_tran_2()); } TEST(mtln, spice_multiple) { EXPECT_EQ(0, test_spice_multiple()); } diff --git a/test/mtln/test_Paul.F90 b/test/mtln/test_Paul.F90 index 5027cb0c..05c46952 100644 --- a/test/mtln/test_Paul.F90 +++ b/test/mtln/test_Paul.F90 @@ -749,7 +749,7 @@ integer function test_coaxial_line_paul_8_6_square() bind(C) result(error_cnt) character(20) :: charR, charL, charC, lineC integer :: i - type(segment_relative_position_t), dimension(100) :: segment_positions + type(external_field_segment_t), dimension(100) :: external_field_segments error_cnt = 0 cable%name = "wire0" @@ -758,9 +758,12 @@ integer function test_coaxial_line_paul_8_6_square() bind(C) result(error_cnt) cable%resistance_per_meter = rpul cable%capacitance_per_meter = cpul do i = 1, 100 - segment_positions(i)%position = (/i, 1, 1/) + external_field_segments(i)%position = (/i, 1, 1/) + external_field_segments(i)%direction = 1 + external_field_segments(i)%Efield_main2wire => null() + external_field_segments(i)%Efield_wire2main => null() end do - cable%segment_relative_positions = segment_positions + cable%external_field_segments = external_field_segments cable%step_size = [(4.0, i = 1, 100)] parsed%time_step = 2e-8 diff --git a/test/mtln/test_mtl_bundle.F90 b/test/mtln/test_mtl_bundle.F90 index 46134fd0..8e66ea58 100644 --- a/test/mtln/test_mtl_bundle.F90 +++ b/test/mtln/test_mtl_bundle.F90 @@ -14,21 +14,24 @@ integer function test_mtl_bundle_init() bind(C) result(error_cnt) real, dimension(5) :: step_size = [20.0, 20.0, 20.0, 20.0, 20.0] - type(segment_relative_position_t), dimension(5) :: segment_positions + type(external_field_segment_t), dimension(5) :: external_field_segments error_cnt = 0 block integer :: i do i = 1, 5 - segment_positions(i)%position = (/i, 1, 1/) + external_field_segments(i)%position = (/i, 1, 1/) + external_field_segments(i)%direction = 1 + external_field_segments(i)%Efield_main2wire => null() + external_field_segments(i)%Efield_wire2main => null() end do end block - mtl_out = mtl_t(l1, c1, r1, g1, step_size, "line_out", segment_relative_positions = segment_positions) + mtl_out = mtl_t(l1, c1, r1, g1, step_size, "line_out", external_field_segments = external_field_segments) mtl_in = mtl_t(l1, c1, r1, g1, step_size, "line_in", & parent_name = "line_out", & conductor_in_parent = 1, & - segment_relative_positions = segment_positions) + external_field_segments = external_field_segments) allocate(levels(1)%lines(1)) allocate(levels(2)%lines(1)) diff --git a/test/mtln/testingTools.F90 b/test/mtln/testingTools.F90 index ffebd140..6e418af1 100644 --- a/test/mtln/testingTools.F90 +++ b/test/mtln/testingTools.F90 @@ -60,7 +60,7 @@ type(mtl_t) function buildLineWithNConductors(n,name, parent_name, conductor_in_ integer, intent(in), optional :: conductor_in_parent real, allocatable, dimension(:,:) :: lpul, cpul, rpul, gpul real, dimension(5) :: step_size = [20.0, 20.0, 20.0, 20.0, 20.0] - type(segment_relative_position_t), dimension(5) :: segment_positions + type(external_field_segment_t), dimension(5) :: external_field_segments integer :: i,j allocate(lpul(n,n), source = 0.0) @@ -69,7 +69,10 @@ type(mtl_t) function buildLineWithNConductors(n,name, parent_name, conductor_in_ allocate(rpul(n,n), source = 0.0) do i = 1, 5 - segment_positions(i)%position =(/i,1,1/) + external_field_segments(i)%position =(/i,1,1/) + external_field_segments(i)%direction = DIRECTION_X_POS + external_field_segments(i)%Efield_main2wire => null() + external_field_segments(i)%Efield_wire2main => null() end do do i = 1, n @@ -87,20 +90,20 @@ type(mtl_t) function buildLineWithNConductors(n,name, parent_name, conductor_in_ end do if (present(dt) .and. .not.present(parent_name)) then res = mtl_t(lpul, cpul, rpul, gpul, step_size, name, dt = dt, & - segment_relative_positions = segment_positions) + external_field_segments = external_field_segments) else if (.not.present(dt) .and. present(parent_name) ) then res = mtl_t(lpul, cpul, rpul, gpul, step_size, name, & parent_name= parent_name, & conductor_in_parent=conductor_in_parent, & - segment_relative_positions = segment_positions) + external_field_segments = external_field_segments) else if (present(dt) .and. present(parent_name) ) then res = mtl_t(lpul, cpul, rpul, gpul, step_size, name, & parent_name= parent_name, & conductor_in_parent=conductor_in_parent, & - dt = dt, segment_relative_positions = segment_positions) + dt = dt, external_field_segments = external_field_segments) else res = mtl_t(lpul, cpul, rpul, gpul, step_size, name, & - segment_relative_positions = segment_positions) + external_field_segments = external_field_segments) end if end function diff --git a/test/pyWrapper/test_full_system.py b/test/pyWrapper/test_full_system.py new file mode 100644 index 00000000..4efcafc4 --- /dev/null +++ b/test/pyWrapper/test_full_system.py @@ -0,0 +1,40 @@ +from utils import * + +def test_holland(tmp_path): + case = 'holland1981' + makeCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') + makeCopy(tmp_path, CASE_FOLDER + case + '.fdtd.json') + fn = tmp_path._str + '/' + case + '.fdtd.json' + + solver = FDTD(input_filename = fn, path_to_exe=SEMBA_EXE) + solver.run() + probe_files = solver.getSolvedProbeFilenames("mid_point") + + assert solver.hasFinishedSuccessfully() == True + p_expected = Probe(OUTPUT_FOLDER+'holland1981.fdtd_mid_point_Wz_11_11_12_s2.dat') + p_solved = Probe(probe_files[0]) + + assert np.allclose(p_expected.df.to_numpy(), p_solved.df.to_numpy()) + + +def test_sphere(tmp_path): + case = 'sphere' + input_json = getCase(case) + input_json['general']['numberOfSteps'] = 200 + input_json['probes'][0]['domain']['numberOfFrequencies'] = 100 + + fn = tmp_path._str + '/' + case + '.fdtd.json' + with open(fn, 'w') as modified_json: + json.dump(input_json, modified_json) + + makeCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') + + solver = FDTD(input_filename = fn, path_to_exe=SEMBA_EXE) + solver.run() + probe_files = solver.getSolvedProbeFilenames("Far") # semba-fdtd seems to always use the name Far for "far field" probes. + + assert solver.hasFinishedSuccessfully() == True + assert len(probe_files) == 1 + + p = Probe(probe_files[0]) + assert p.type == 'farField' \ No newline at end of file diff --git a/test/pyWrapper/test_pyWrapper.py b/test/pyWrapper/test_integration.py similarity index 70% rename from test/pyWrapper/test_pyWrapper.py rename to test/pyWrapper/test_integration.py index 2d4dccdd..91753646 100644 --- a/test/pyWrapper/test_pyWrapper.py +++ b/test/pyWrapper/test_integration.py @@ -1,5 +1,23 @@ from utils import * +def test_read_wire_probe(): + p = Probe(OUTPUT_FOLDER + 'holland1981.fdtd_mid_point_Wz_11_11_12_s2.dat') + + assert p.case_name == 'holland1981' + assert p.name == 'mid_point' + assert p.type == 'wire' + assert np.all(p.cell == np.array([11, 11, 12])) + assert p.segment_tag == 2 + + assert len(p['time']) == 1001 + assert p['time'][0] == 0.0 + assert p['time'].iat[-1] == 0.99999999600419720E-009 + + assert len(p['current']) == 1001 + assert p['current'][0] == 0.0 + assert p['current'].iat[-1] == -0.228888888E-021 + + def test_probes_output_exists(tmp_path): case = 'holland1981' input_json = getCase(case) @@ -8,7 +26,7 @@ def test_probes_output_exists(tmp_path): with open(fn, 'w') as modified_json: json.dump(input_json, modified_json) - makeTemporaryCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') + makeCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') solver = FDTD(input_filename = fn, path_to_exe=SEMBA_EXE) solver.run() @@ -29,7 +47,7 @@ def test_probes_output_number_of_steps(tmp_path): with open(fn, 'w') as modified_json: json.dump(input_json, modified_json) - makeTemporaryCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') + makeCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') solver = FDTD(input_filename = fn, path_to_exe=SEMBA_EXE) solver.run() @@ -42,8 +60,8 @@ def test_probes_output_number_of_steps(tmp_path): def test_holland(tmp_path): case = 'holland1981' - makeTemporaryCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') - makeTemporaryCopy(tmp_path, CASE_FOLDER + case + '.fdtd.json') + makeCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') + makeCopy(tmp_path, CASE_FOLDER + case + '.fdtd.json') fn = tmp_path._str + '/' + case + '.fdtd.json' solver = FDTD(input_filename = fn, path_to_exe=SEMBA_EXE) @@ -69,7 +87,7 @@ def test_towel_hanger(tmp_path): with open(fn, 'w') as modified_json: json.dump(input_json, modified_json) - makeTemporaryCopy(tmp_path, EXCITATIONS_FOLDER+'ramp.exc') + makeCopy(tmp_path, EXCITATIONS_FOLDER+'ramp.exc') solver = FDTD(input_filename = fn, path_to_exe=SEMBA_EXE) solver.run() @@ -88,22 +106,30 @@ def test_towel_hanger(tmp_path): # probe_files[0]) -def test_sphere(tmp_path): +def test_read_far_field_probe(tmp_path): case = 'sphere' input_json = getCase(case) input_json['general']['numberOfSteps'] = 1 - input_json['general']['timeStep'] = 3.0E-013 + input_json['probes'][0]['domain']['numberOfFrequencies'] = 100 + fn = tmp_path._str + '/' + case + '.fdtd.json' with open(fn, 'w') as modified_json: json.dump(input_json, modified_json) - makeTemporaryCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') + makeCopy(tmp_path, EXCITATIONS_FOLDER+'gauss.exc') solver = FDTD(input_filename = fn, path_to_exe=SEMBA_EXE) - solver.run() - probe_files = solver.getSolvedProbeFilenames("Far") # semba-fdtd seems to always use the name Far for "far field" probes. + solver.run() + + p = Probe(solver.getSolvedProbeFilenames("Far")[0]) + assert p.case_name == 'sphere' + assert p.type == 'farField' + assert np.all(p.cell_init == np.array([2,2,2])) - assert solver.hasFinishedSuccessfully() == True - assert len(probe_files) == 1 - assert 'sphere.fdtd_Far_FF_2_2_2__77_77_77.dat' == probe_files[0] \ No newline at end of file + p = Probe(solver.getSolvedProbeFilenames("electric_field_movie")[0]) + assert p.case_name == 'sphere' + assert p.type == 'movie' + assert np.all(p.cell_init == np.array([2,2,2])) + + \ No newline at end of file diff --git a/test/pyWrapper/utils.py b/test/pyWrapper/utils.py index 00af622b..40a5c22a 100644 --- a/test/pyWrapper/utils.py +++ b/test/pyWrapper/utils.py @@ -1,4 +1,3 @@ -import pytest from src_pyWrapper.pyWrapper import * import shutil, glob, re import json @@ -17,15 +16,15 @@ def getCase(case): return json.load(open(CASE_FOLDER + case + '.fdtd.json')) -def makeTemporaryCopy(temp_dir, file_path): +def makeCopy(temp_dir, file_path): file_name = file_path.split('/')[-1] temp_path = os.path.join(temp_dir, file_name) shutil.copy2(file_path, temp_path) -def copyTemporaryInputFiles(temp_dir, input, excitation, executable): - makeTemporaryCopy(temp_dir, input) - makeTemporaryCopy(temp_dir, excitation) - makeTemporaryCopy(temp_dir, executable) +def copyInputFiles(temp_dir, input, excitation, executable): + makeCopy(temp_dir, input) + makeCopy(temp_dir, excitation) + makeCopy(temp_dir, executable) def getProbeFile(prefix, probe_name): return ([x for x in glob.glob('*dat') if re.match(prefix + '_' + probe_name + '.*dat',x)])[0] @@ -43,28 +42,7 @@ def readWithoutHeader(file_name): def compareFiles(expected_name, result_name): f_expected = readWithoutHeader(expected_name) f_result = readWithoutHeader(result_name) - return f_expected == f_result - - -def readTimeAndExcitation(exc_file): - t = np.array([]) - e = np.array([]) - with open(exc_file, 'r') as f: - for line in f: - t = np.append(t, float(line.split()[0])) - e = np.append(e, float(line.split()[1])) - return t, e - -def readProbeTimeAndValue(probe_file): - t = np.array([]) - value = np.array([]) - with open(probe_file, 'r') as f: - f.readline() - for line in f: - t = np.append(t, float(line.split()[0])) - value = np.append(value, float(line.split()[1])) - return t, value - + return f_expected == f_result diff --git a/test/smbjson/test_read_mtln.F90 b/test/smbjson/test_read_mtln.F90 index 6edf0e65..e81f046f 100644 --- a/test/smbjson/test_read_mtln.F90 +++ b/test/smbjson/test_read_mtln.F90 @@ -333,9 +333,12 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(1)%step_size(9)) expected%mtln%cables(1)%step_size = [(0.1, i = 1, 9)] - allocate(expected%mtln%cables(1)%segment_relative_positions(9)) + allocate(expected%mtln%cables(1)%external_field_segments(9)) do i = 1, 9 - expected%mtln%cables(1)%segment_relative_positions(i)%position = (/i,9,1/) + expected%mtln%cables(1)%external_field_segments(i)%position = (/i,9,1/) + expected%mtln%cables(1)%external_field_segments(i)%direction = DIRECTION_X_POS + expected%mtln%cables(1)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(1)%external_field_segments(i)%Efield_wire2main => null() end do allocate(expected%mtln%cables(1)%transfer_impedance%poles(0)) @@ -359,9 +362,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(2)%step_size(9)) expected%mtln%cables(2)%step_size = [(0.1, i = 1, 9)] - allocate(expected%mtln%cables(2)%segment_relative_positions(9)) + allocate(expected%mtln%cables(2)%external_field_segments(9)) do i = 1, 9 - expected%mtln%cables(2)%segment_relative_positions(i)%position = (/i,9,1/) + expected%mtln%cables(2)%external_field_segments(i)%position = (/i,9,1/) + expected%mtln%cables(2)%external_field_segments(i)%direction = DIRECTION_X_POS + expected%mtln%cables(2)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(2)%external_field_segments(i)%Efield_wire2main => null() + end do expected%mtln%cables(2)%transfer_impedance%direction = TRANSFER_IMPEDANCE_DIRECTION_INWARDS @@ -406,9 +413,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(3)%step_size(9)) expected%mtln%cables(3)%step_size = [(0.1, i = 1, 9)] - allocate(expected%mtln%cables(3)%segment_relative_positions(9)) + allocate(expected%mtln%cables(3)%external_field_segments(9)) do i = 1, 9 - expected%mtln%cables(3)%segment_relative_positions(i)%position = (/i,9,1/) + expected%mtln%cables(3)%external_field_segments(i)%position = (/i,9,1/) + expected%mtln%cables(3)%external_field_segments(i)%direction = DIRECTION_X_POS + expected%mtln%cables(3)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(3)%external_field_segments(i)%Efield_wire2main => null() + end do @@ -439,9 +450,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(4)%step_size(8)) expected%mtln%cables(4)%step_size = [(0.1, i = 1, 8)] - allocate(expected%mtln%cables(4)%segment_relative_positions(8)) + allocate(expected%mtln%cables(4)%external_field_segments(8)) do i = 1, 8 - expected%mtln%cables(4)%segment_relative_positions(i)%position = (/9+i,9,1/) + expected%mtln%cables(4)%external_field_segments(i)%position = (/9+i,9,1/) + expected%mtln%cables(4)%external_field_segments(i)%direction = DIRECTION_X_POS + expected%mtln%cables(4)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(4)%external_field_segments(i)%Efield_wire2main => null() + end do @@ -466,9 +481,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(5)%step_size(8)) expected%mtln%cables(5)%step_size = [(0.1, i = 1, 8)] - allocate(expected%mtln%cables(5)%segment_relative_positions(8)) + allocate(expected%mtln%cables(5)%external_field_segments(8)) do i = 1, 8 - expected%mtln%cables(5)%segment_relative_positions(i)%position = (/9+i,9,1/) + expected%mtln%cables(5)%external_field_segments(i)%position = (/9+i,9,1/) + expected%mtln%cables(5)%external_field_segments(i)%direction = DIRECTION_X_POS + expected%mtln%cables(5)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(5)%external_field_segments(i)%Efield_wire2main => null() + end do @@ -502,9 +521,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(6)%step_size(8)) expected%mtln%cables(6)%step_size = [(0.1, i = 1, 8)] - allocate(expected%mtln%cables(6)%segment_relative_positions(8)) + allocate(expected%mtln%cables(6)%external_field_segments(8)) do i = 1, 8 - expected%mtln%cables(6)%segment_relative_positions(i)%position = (/9+i,9,1/) + expected%mtln%cables(6)%external_field_segments(i)%position = (/9+i,9,1/) + expected%mtln%cables(6)%external_field_segments(i)%direction = DIRECTION_X_POS + expected%mtln%cables(6)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(6)%external_field_segments(i)%Efield_wire2main => null() + end do @@ -535,9 +558,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(7)%step_size(7)) expected%mtln%cables(7)%step_size = [(0.1, i = 1, 7)] - allocate(expected%mtln%cables(7)%segment_relative_positions(7)) + allocate(expected%mtln%cables(7)%external_field_segments(7)) do i = 1,7 - expected%mtln%cables(7)%segment_relative_positions(i)%position = (/10,9-i,1/) + expected%mtln%cables(7)%external_field_segments(i)%position = (/10,9-i,1/) + expected%mtln%cables(7)%external_field_segments(i)%direction = DIRECTION_Y_NEG + expected%mtln%cables(7)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(7)%external_field_segments(i)%Efield_wire2main => null() + end do @@ -562,9 +589,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(8)%step_size(7)) expected%mtln%cables(8)%step_size = [(0.1, i = 1, 7)] - allocate(expected%mtln%cables(8)%segment_relative_positions(7)) + allocate(expected%mtln%cables(8)%external_field_segments(7)) do i = 1,7 - expected%mtln%cables(8)%segment_relative_positions(i)%position = (/10,9-i,1/) + expected%mtln%cables(8)%external_field_segments(i)%position = (/10,9-i,1/) + expected%mtln%cables(8)%external_field_segments(i)%direction = DIRECTION_Y_NEG + expected%mtln%cables(8)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(8)%external_field_segments(i)%Efield_wire2main => null() + end do expected%mtln%cables(8)%transfer_impedance%direction = TRANSFER_IMPEDANCE_DIRECTION_INWARDS @@ -604,9 +635,13 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(9)%step_size(7)) expected%mtln%cables(9)%step_size = [(0.1, i = 1, 7)] - allocate(expected%mtln%cables(9)%segment_relative_positions(7)) + allocate(expected%mtln%cables(9)%external_field_segments(7)) do i = 1,7 - expected%mtln%cables(9)%segment_relative_positions(i)%position = (/10,9-i,1/) + expected%mtln%cables(9)%external_field_segments(i)%position = (/10,9-i,1/) + expected%mtln%cables(9)%external_field_segments(i)%direction = DIRECTION_Y_NEG + expected%mtln%cables(9)%external_field_segments(i)%Efield_main2wire => null() + expected%mtln%cables(9)%external_field_segments(i)%Efield_wire2main => null() + end do expected%mtln%cables(9)%transfer_impedance%direction = TRANSFER_IMPEDANCE_DIRECTION_INWARDS diff --git a/test/smbjson/test_read_planewave.F90 b/test/smbjson/test_read_planewave.F90 index 5d544237..7b80a792 100644 --- a/test/smbjson/test_read_planewave.F90 +++ b/test/smbjson/test_read_planewave.F90 @@ -4,17 +4,17 @@ integer function test_read_planewave() bind (C) result(err) implicit none - character(len=*),parameter :: filename = PATH_TO_TEST_DATA//'cases/planewave.fdtd.json' - type(Parseador) :: problem, expected + character(len=*), parameter :: filename = PATH_TO_TEST_DATA//'cases/planewave.fdtd.json' + type(Parseador) :: pr, ex type(parser_t) :: parser logical :: areSame err = 0 - expected = expectedProblemDescription() + ex = expectedProblemDescription() parser = parser_t(filename) - problem = parser%readProblemDescription() - call expect_eq(err, expected, problem) - + pr = parser%readProblemDescription() + call expect_eq(err, ex, pr) + contains function expectedProblemDescription() result (expected) type(Parseador) :: expected diff --git a/test/smbjson/test_read_shieldedPair.F90 b/test/smbjson/test_read_shieldedPair.F90 index fbb25a51..d2b5c893 100644 --- a/test/smbjson/test_read_shieldedPair.F90 +++ b/test/smbjson/test_read_shieldedPair.F90 @@ -148,9 +148,13 @@ function expectedProblemDescription() result (expected) expected%mtln%cables(2)%conductance_per_meter = reshape(source=[0.0], shape=[1,1]) allocate(expected%mtln%cables(2)%step_size(18)) expected%mtln%cables(2)%step_size = [(0.03, i = 1, 18)] - allocate(expected%mtln%cables(2)%segment_relative_positions(18)) + allocate(expected%mtln%cables(2)%external_field_segments(18)) do i = 1, 18 - expected%mtln%cables(2)%segment_relative_positions(i)%position = (/1,1,i/) + expected%mtln%cables(2)%external_field_segments(i)%position = (/1,1,i/) + expected%mtln%cables(2)%external_field_segments(i)%direction = DIRECTION_Z_POS + expected%mtln%cables(2)%external_field_segments(i)%Efield_wire2main => null() + expected%mtln%cables(2)%external_field_segments(i)%Efield_main2wire => null() + end do allocate(expected%mtln%cables(2)%transfer_impedance%poles(0)) @@ -177,9 +181,12 @@ function expectedProblemDescription() result (expected) allocate(expected%mtln%cables(1)%step_size(18)) expected%mtln%cables(1)%step_size = [(0.03, i = 1, 18)] - allocate(expected%mtln%cables(1)%segment_relative_positions(18)) + allocate(expected%mtln%cables(1)%external_field_segments(18)) do i = 1, 18 - expected%mtln%cables(1)%segment_relative_positions(i)%position = (/1,1,i/) + expected%mtln%cables(1)%external_field_segments(i)%position = (/1,1,i/) + expected%mtln%cables(1)%external_field_segments(i)%direction = DIRECTION_Z_POS + expected%mtln%cables(1)%external_field_segments(i)%Efield_wire2main => null() + expected%mtln%cables(1)%external_field_segments(i)%Efield_main2wire => null() end do expected%mtln%cables(1)%transfer_impedance%direction = TRANSFER_IMPEDANCE_DIRECTION_INWARDS expected%mtln%cables(1)%transfer_impedance%resistive_term = 0.0 diff --git a/test/smbjson/test_read_sphere.F90 b/test/smbjson/test_read_sphere.F90 index 53e934db..6505aa1d 100644 --- a/test/smbjson/test_read_sphere.F90 +++ b/test/smbjson/test_read_sphere.F90 @@ -14,16 +14,7 @@ integer function test_read_sphere() bind (C) result(err) parser = parser_t(filename) pr = parser%readProblemDescription() - if (.not. ex%general == pr%general) call testFails(err, 'Expected and read "general" do not match') - if (.not. ex%matriz == pr%matriz) call testFails(err, 'Expected and read "media matrix" do not match') - if (.not. ex%despl == pr%despl) call testFails(err, 'Expected and read "grid" do not match') - if (.not. ex%front == pr%front) call testFails(err, 'Expected and read "boundary" do not match') - if (.not. ex%Mats == pr%Mats) call testFails(err, 'Expected and read "materials" do not match') - ! -- specific surfs not included DO NOT use comparison -- - if (.not. ex%plnSrc == pr%plnSrc) call testFails(err, 'Expected and read "planewave sources" do not match') - if (.not. ex%nodSrc == pr%nodSrc) call testFails(err, 'Expected and read "nodal sources" do not match') - if (.not. ex%sonda == pr%sonda) call testFails(err, 'Expected and read "new probes" do not match') - if (.not. ex%BloquePrb == pr%BloquePrb) call testFails(err, 'Expected and read "block probes" do not match') + call expect_eq(err, ex, pr, ignoreRegions=.true.) if (err == 0) write(*,*) "Read and expected inputs are equal." @@ -96,7 +87,11 @@ function expectedProblemDescription() result (ex) ex%oldSonda%probes(1)%n_FarField = 1 ex%oldSonda%probes(1)%n_FarField_max = 1 allocate(ex%oldSonda%probes(1)%FarField(1)) - ex%oldSonda%probes(1)%FarField(1)%probe%outputrequest = "Far field_log_" + ex%oldSonda%probes(1)%FarField(1)%probe%grname = " " + ex%oldSonda%probes(1)%FarField(1)%probe%outputrequest = "FarField_log_" + ex%oldSonda%probes(1)%FarField(1)%probe%tstart = 0.0 + ex%oldSonda%probes(1)%FarField(1)%probe%tstop = 0.0 + ex%oldSonda%probes(1)%FarField(1)%probe%tstep = 0.0 ex%oldSonda%probes(1)%FarField(1)%probe%fstart = 1e6 ex%oldSonda%probes(1)%FarField(1)%probe%fstop = 1e9 ex%oldSonda%probes(1)%FarField(1)%probe%fstep = 1e6*5 @@ -104,9 +99,10 @@ function expectedProblemDescription() result (ex) allocate(ex%oldSonda%probes(1)%FarField(1)%probe%i(2)) allocate(ex%oldSonda%probes(1)%FarField(1)%probe%j(2)) allocate(ex%oldSonda%probes(1)%FarField(1)%probe%k(2)) - ex%oldSonda%probes(1)%FarField(1)%probe%i = [2, 78] - ex%oldSonda%probes(1)%FarField(1)%probe%j = [2, 78] - ex%oldSonda%probes(1)%FarField(1)%probe%k = [2, 78] + allocate(ex%oldSonda%probes(1)%FarField(1)%probe%node(0)) + ex%oldSonda%probes(1)%FarField(1)%probe%i = [2, 77] + ex%oldSonda%probes(1)%FarField(1)%probe%j = [2, 77] + ex%oldSonda%probes(1)%FarField(1)%probe%k = [2, 77] ex%oldSonda%probes(1)%FarField(1)%probe%n_cord = 2 ex%oldSonda%probes(1)%FarField(1)%probe%n_cord_max = 2 ex%oldSonda%probes(1)%FarField(1)%probe%thetastart = 0.0 @@ -116,6 +112,33 @@ function expectedProblemDescription() result (ex) ex%oldSonda%probes(1)%FarField(1)%probe%phistop = 360.0 ex%oldSonda%probes(1)%FarField(1)%probe%phistep = 90.0 + allocate(ex%VolPrb) + ex%VolPrb%length = 1 + ex%VolPrb%length_max = 1 + ex%VolPrb%len_cor_max = 2 + allocate(ex%VolPrb%collection(1)) + allocate(ex%VolPrb%collection(1)%cordinates(1)) + ex%VolPrb%collection(1)%len_cor = 1 + ex%VolPrb%collection(1)%cordinates(1)%Xi = 2 + ex%VolPrb%collection(1)%cordinates(1)%Xe = 77 + ex%VolPrb%collection(1)%cordinates(1)%Yi = 2 + ex%VolPrb%collection(1)%cordinates(1)%Ye = 77 + ex%VolPrb%collection(1)%cordinates(1)%Zi = 2 + ex%VolPrb%collection(1)%cordinates(1)%Ze = 77 + ex%VolPrb%collection(1)%cordinates(1)%or = iExC + ex%VolPrb%collection(1)%cordinates(1)%xtrancos = 1 + ex%VolPrb%collection(1)%cordinates(1)%ytrancos = 1 + ex%VolPrb%collection(1)%cordinates(1)%ztrancos = 1 + ex%VolPrb%collection(1)%cordinates(1)%tag = "" + ex%VolPrb%collection(1)%tstart = 0.0 + ex%VolPrb%collection(1)%tstop = 0.0 + ex%VolPrb%collection(1)%tstep = 1e-9 + ex%VolPrb%collection(1)%fstart = 0.0 + ex%VolPrb%collection(1)%fstop = 0.0 + ex%VolPrb%collection(1)%fstep = 0.0 + ex%VolPrb%collection(1)%outputrequest = "electric_field_movie" + ex%VolPrb%collection(1)%filename = " " + ex%VolPrb%collection(1)%type2 = NP_T2_TIME end function end function diff --git a/test/smbjson/testingTools.F90 b/test/smbjson/testingTools.F90 index 7e921e70..dcb6334c 100644 --- a/test/smbjson/testingTools.F90 +++ b/test/smbjson/testingTools.F90 @@ -11,41 +11,47 @@ subroutine expect_eq_int(err, ex, pr, msg) if (ex /= pr) call testFails(err, msg) end subroutine - subroutine expect_eq(err, ex, pr) + subroutine expect_eq(err, ex, pr, ignoreRegions) integer, intent(inout) :: err type(Parseador), intent(in) :: ex, pr + logical, optional, intent(in) :: ignoreRegions + logical :: checkRegions + if (.not. present(ignoreRegions)) then + checkRegions = .true. + else + if (ignoreRegions) then + checkRegions = .false. + else + checkRegions = .true. + end if + end if ! Basics if (.not. ex%general == pr%general) call testFails(err, 'Expected and read "general" do not match') if (.not. ex%matriz == pr%matriz) call testFails(err, 'Expected and read "media matrix" do not match') if (.not. ex%despl == pr%despl) call testFails(err, 'Expected and read "grid" do not match') if (.not. ex%front == pr%front) call testFails(err, 'Expected and read "boundary" do not match') + ! Materials if (.not. ex%Mats == pr%Mats) call testFails(err, 'Expected and read "materials" do not match') - if (.not. ex%pecRegs == pr%pecRegs) call testFails(err, 'Expected and read "pec regions" do not match') - if (.not. ex%pmcRegs == pr%pmcRegs) call testFails(err, 'Expected and read "pmc regions" do not match') - ! if (.not. ex%DielRegs == pr%DielRegs) & - ! call testFails(error_cnt, 'Expected and read "dielectric regions" do not match') - ! if (.not. ex%LossyThinSurfs == pr%LossyThinSurfs) & - ! call testFails(error_cnt, 'Expected and read "lossy thin surfs" do not match') - ! if (.not. ex%frqDepMats == pr%frqDepMats) & - ! call testFails(error_cnt, 'Expected and read "frq. dep. materials" do not match') - ! if (.not. ex%aniMats == pr%aniMats) & - ! call testFails(error_cnt, 'Expected and read "anisotropic materials" do not match') + if (checkRegions) then + if (.not. ex%pecRegs == pr%pecRegs) call testFails(err, 'Expected and read "pec regions" do not match') + if (.not. ex%pmcRegs == pr%pmcRegs) call testFails(err, 'Expected and read "pmc regions" do not match') + end if + ! Sources - ! if (.not. ex%boxSrc == pr%boxSrc) call testFails(error_cnt, 'Expected and read "box sources" do not match') + ! if (.not. ex%boxSrc == pr%boxSrc) call testFails(err, 'Expected and read "box sources" do not match') if (.not. ex%plnSrc == pr%plnSrc) call testFails(err, 'Expected and read "planewave sources" do not match') if (.not. ex%nodSrc == pr%nodSrc) call testFails(err, 'Expected and read "nodal sources" do not match') + ! Probes if (.not. ex%oldSONDA == pr%oldSonda) call testFails(err, 'Expected and read "old probes" do not match') if (.not. ex%sonda == pr%sonda) call testFails(err, 'Expected and read "new probes" do not match') if (.not. ex%BloquePrb == pr%BloquePrb) call testFails(err, 'Expected and read "block probes" do not match') - ! if (.not. ex%VolPrb == pr%VolPrb) call testFails(error_cnt, 'Expected and read "vol probes" do not match') + if (.not. ex%VolPrb == pr%VolPrb) call testFails(err, 'Expected and read "vol probes" do not match') ! Thin elements if (.not. ex%tWires == pr%tWires) call testFails(err, 'Expected and read "thin wires" do not match') - ! if (.not. ex%sWires == pr%sWires) call testFails(error_cnt, 'Expected and read "slanted wires" do not match') - ! if (.not. ex%tSlots == pr%tSlots) call testFails(error_cnt, 'Expected and read "thin slots" do not match') - + if (err == 0) write(*,*) "Read and expected inputs are equal." end subroutine diff --git a/testData/cases/sphere.fdtd.json b/testData/cases/sphere.fdtd.json index ed9ed5f6..72172030 100644 --- a/testData/cases/sphere.fdtd.json +++ b/testData/cases/sphere.fdtd.json @@ -40951,6 +40951,17 @@ "numberOfFrequencies": 5, "frequencySpacing": "logarithmic" } + }, + { + "name": "electric_field_movie", + "type": "movie", + "field": "electric", + "components": ["x"], + "elementIds": [2], + "domain": { + "type": "time", + "samplingPeriod": 1e-9 + } } ] } \ No newline at end of file