diff --git a/src_output/bulkProbeOutput.F90 b/src_output/bulkProbeOutput.F90 new file mode 100644 index 00000000..d97626c8 --- /dev/null +++ b/src_output/bulkProbeOutput.F90 @@ -0,0 +1,183 @@ +module bulkProbeOutput_m + use FDETYPES_m + use utils_m + use outputTypes_m + use outputUtils_m + implicit none + +contains + + subroutine init_bulk_probe_output(this, lowerBound, upperBound, field, domain, outputTypeExtension, mpidir) + type(bulk_current_probe_output_t), intent(out) :: this + type(cell_coordinate_t), intent(in) :: lowerBound, upperBound + integer(kind=SINGLE), intent(in) :: mpidir, field + character(len=BUFSIZE), intent(in) :: outputTypeExtension + type(domain_t), intent(in) :: domain + + integer(kind=SINGLE) :: i + + this%mainCoords = lowerBound + this%auxCoords = upperBound + this%component = field + + this%domain = domain + this%path = get_output_path() + + call alloc_and_init(this%timeStep, BuffObse, 0.0_RKIND_tiempo) + call alloc_and_init(this%valueForTime, BuffObse, 0.0_RKIND) + call create_data_file(this%filePathTime, this%path, timeExtension, datFileExtension) + + contains + + function get_output_path() result(outputPath) + character(len=BUFSIZE) :: probeBoundsExtension, prefixFieldExtension + character(len=BUFSIZE) :: outputPath + probeBoundsExtension = get_coordinates_extension(this%mainCoords, this%auxCoords, mpidir) + prefixFieldExtension = get_prefix_extension(field, mpidir) + outputPath = & + trim(adjustl(outputTypeExtension))//'_'//trim(adjustl(prefixFieldExtension))//'_'//trim(adjustl(probeBoundsExtension)) + return + end function get_output_path + + end subroutine init_bulk_probe_output + + subroutine update_bulk_probe_output(this, step, field) + type(bulk_current_probe_output_t), intent(inout) :: this + real(kind=RKIND_tiempo), intent(in) :: step + type(field_data_t), intent(in) :: field + + integer(kind=SINGLE) :: i1_m, i2_m, j1_m, j2_m, k1_m, k2_m + integer(kind=SINGLE) :: i1, i2, j1, j2, k1, k2 + integer(kind=SINGLE) :: iii, jjj, kkk + + real(kind=RKIND), pointer, dimension(:, :, :) :: xF, yF, zF + real(kind=RKIND), pointer, dimension(:) :: dx, dy, dz + + i1_m = this%mainCoords%x + j1_m = this%mainCoords%y + k1_m = this%mainCoords%z + i2_m = this%auxCoords%x + j2_m = this%auxCoords%y + k2_m = this%auxCoords%z + + i1 = i1_m + j1 = j1_m + k1 = k1_m + i2 = i2_m + j2 = j2_m + k2 = k2_m + + xF => field%x + yF => field%y + zF => field%z + dx => field%deltaX + dy => field%deltaY + dz => field%deltaZ + + this%nTime = this%nTime + 1 + this%timeStep(this%nTime) = step + this%valueForTime(this%nTime) = 0.0_RKIND !Clear uninitialized value + selectcase (this%component) + case (iBloqueJx) + do JJJ = j1, j2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (yF(i1_m, JJJ, k1_m - 1) - yF(i1_m, JJJ, k2_m))*dy(JJJ) + end do + do KKK = k1, k2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (-zF(i1_m, j1_m - 1, KKK) + zF(i1_m, j2_m, KKK))*dz(KKK) + end do + + case (iBloqueJy) + do KKK = k1, k2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (-zF(i2_m, j1_m, KKK) + zF(i1_m - 1, j1_m, KKK))*dz(KKK) + end do + do III = i1, i2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (xF(III, j1_m, k2_m) - xF(III, j1_m, k1_m - 1))*dx(III) + end do + + case (iBloqueJz) + do III = i1, i2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (xF(III, j1_m - 1, k1_m) - xF(III, j2_m, k1_m))*dx(III) + end do + do JJJ = j1, j2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (-yF(i1_m - 1, JJJ, k1_m) + yF(i2_m, JJJ, k1_m))*dy(JJJ) + end do + + case (iBloqueMx) + do JJJ = j1, j2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (-yF(i1_m, JJJ, k1_m) + yF(i1_m, JJJ, k2_m + 1))*dy(JJJ) + end do + do KKK = k1, k2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (zF(i1_m, j1_m, KKK) - zF(i1_m, j2_m + 1, KKK))*dz(KKK) + end do + + case (iBloqueMy) + do KKK = k1, k2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (zF(i2_m + 1, j1_m, KKK) - zF(i1_m, j1_m, KKK))*dz(KKK) + end do + do III = i1, i2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (-xF(III, j1_m, k2_m + 1) + xF(III, j1_m, k1_m))*dx(III) + end do + + case (iBloqueMz) + do III = i1, i2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (-xF(III, j1_m, k1_m) + xF(III, j2_m + 1, k1_m))*dx(III) + end do + do JJJ = j1, j2 + this%valueForTime(this%nTime) = & + this%valueForTime(this%nTime) + & + (yF(i1_m, JJJ, k1_m) - yF(i2_m + 1, JJJ, k1_m))*dy(JJJ) + end do + + end select + + end subroutine update_bulk_probe_output + + subroutine flush_bulk_probe_output(this) + type(bulk_current_probe_output_t), intent(inout) :: this + integer :: i + integer :: unit + if (this%nTime <= 0) then + print *, "No data to write." + return + end if + + open (unit=unit, file=this%filePathTime, status="old", action="write", position="append") + + do i = 1, this%nTime + write (unit, fmt) this%timeStep(i), this%valueForTime(i) + end do + + close (unit) + call clear_time_data() + contains + subroutine clear_time_data() + this%timeStep = 0.0_RKIND_tiempo + this%valueForTime = 0.0_RKIND + + this%nTime = 0 + end subroutine clear_time_data + end subroutine flush_bulk_probe_output + +end module bulkProbeOutput_m diff --git a/src_output/farFieldProbeOutput.F90 b/src_output/farFieldProbeOutput.F90 new file mode 100644 index 00000000..59b21a72 --- /dev/null +++ b/src_output/farFieldProbeOutput.F90 @@ -0,0 +1,94 @@ +module farFieldOutput_m + use outputTypes_m + use report_m + use outputUtils_m + use farfield_m + implicit none + private + !=========================== + ! Public interface summary + !=========================== + public :: init_farField_probe_output + !public :: create_farField_probe_output + public :: update_farField_probe_output + public :: flush_farField_probe_output + !=========================== +contains + + subroutine init_farField_probe_output(this, sgg, lowerBound, upperBound, field, domain, sphericRange, outputTypeExtension, fileNormalize,control, problemInfo, eps0, mu0) + type(far_field_probe_output_t), intent(out) :: this + type(domain_t), intent(in) :: domain + type(SGGFDTDINFO_t), intent(in) :: sgg + type(cell_coordinate_t), intent(in) :: lowerBound, upperBound + integer(kind=SINGLE), intent(in) :: field + type(spheric_domain_t), intent(in) :: sphericRange + type(sim_control_t), intent(in) :: control + character(len=*), intent(in) :: fileNormalize, outputTypeExtension + type(problem_info_t), intent(in) :: problemInfo + real(kind=RKIND), intent(in) :: mu0, eps0 + + if (domain%domainType /= TIME_DOMAIN) call StopOnError(0, 0, "Unexpected domain type for farField probe") + + this%domain = domain + this%sphericRange = sphericRange + this%component = field + this%path = get_output_path() + call InitFarField(sgg, & + problemInfo%geometryToMaterialData%sggMiEx, problemInfo%geometryToMaterialData%sggMiEy, problemInfo%geometryToMaterialData%sggMiEz, & + problemInfo%geometryToMaterialData%sggMiHx, problemInfo%geometryToMaterialData%sggMiHy, problemInfo%geometryToMaterialData%sggMiHz, & + control%layoutnumber, control%num_procs, problemInfo%simulationBounds, control%resume, & + 2025, this%path, & + lowerBound%x, upperBound%x, & + lowerBound%y, upperBound%y, & + lowerBound%z, upperBound%z, & + domain%fstart, domain%fstop, domain%fstep, & + sphericRange%phiStart, sphericRange%phiStop, sphericRange%phiStep, & + sphericRange%thetaStart, sphericRange%thetaStop, sphericRange%thetaStep, & + fileNormalize, problemInfo%problemDimension, & + control%facesNF2FF, control%NF2FFDecim, & +#ifdef CompileWithMPI + 0, 0, & +#endif + eps0, mu0) + + contains + function get_output_path() result(outputPath) + character(len=BUFSIZE) :: probeBoundsExtension, prefixFieldExtension + character(len=BUFSIZE) :: outputPath + probeBoundsExtension = get_coordinates_extension(this%mainCoords, this%auxCoords, control%mpidir) + prefixFieldExtension = get_prefix_extension(field, control%mpidir) + outputPath = & + trim(adjustl(outputTypeExtension))//'_'//trim(adjustl(prefixFieldExtension))//'_'//trim(adjustl(probeBoundsExtension)) + return + end function get_output_path + + end subroutine init_farField_probe_output + + subroutine update_farField_probe_output(this, ntime, bounds, fieldsReference) + type(far_field_probe_output_t), intent(inout) :: this + type(fields_reference_t), intent(in) :: fieldsReference + integer(kind=SINGLE), intent(in) :: ntime + type(bounds_t), intent(in) :: bounds + call UpdateFarField(ntime, bounds, & + fieldsReference%E%x, fieldsReference%E%y, fieldsReference%E%z, & + fieldsReference%H%x, fieldsReference%H%y, fieldsReference%H%z) + end subroutine update_farField_probe_output + + subroutine flush_farField_probe_output(this, simlulationTimeArray, timeIndex, control, fieldsReference, bounds) + type(far_field_probe_output_t), intent(out) :: this + integer, intent(in) :: timeIndex + real(KIND=RKIND_tiempo), pointer, dimension(:), intent(in) :: simlulationTimeArray + type(sim_control_t), intent(in) :: control + type(fields_reference_t), pointer, intent(in) :: fieldsReference + type(bounds_t), intent(in) :: bounds + + real(kind=RKIND_tiempo) :: flushTime + + flushTime = simlulationTimeArray(timeIndex) + call FlushFarfield(control%layoutnumber, control%num_procs, bounds, & + fieldsReference%E%deltaX, fieldsReference%E%deltaY, fieldsReference%E%deltaZ, & + fieldsReference%H%deltaX, fieldsReference%H%deltaY, fieldsReference%H%deltaZ, & + control%facesNF2FF, flushTime) + end subroutine flush_farfield_probe_output + +end module farFieldOutput_m diff --git a/src_output/pointProbeOutput.F90 b/src_output/pointProbeOutput.F90 new file mode 100644 index 00000000..5dd4e815 --- /dev/null +++ b/src_output/pointProbeOutput.F90 @@ -0,0 +1,162 @@ +module pointProbeOutput_m + use FDETYPES_m + use utils_m + use outputTypes_m + use domain_m + use outputUtils_m + + implicit none + + private + + public :: init_point_probe_output + public :: update_point_probe_output + public :: flush_point_probe_output + +contains + subroutine init_point_probe_output(this, coordinates, field, domain, outputTypeExtension, mpidir, timeInterval) + type(point_probe_output_t), intent(out) :: this + type(cell_coordinate_t) :: coordinates + integer(kind=SINGLE), intent(in) :: mpidir, field + character(len=*), intent(in) :: outputTypeExtension + type(domain_t), intent(in) :: domain + + real(kind=RKIND_tiempo), intent(in) :: timeInterval + + integer(kind=SINGLE) :: i + + this%mainCoords = coordinates + + this%component = field + + this%domain = domain + this%path = get_output_path() + + if (any(this%domain%domainType == (/TIME_DOMAIN, BOTH_DOMAIN/))) then + call alloc_and_init(this%timeStep, BUFSIZE, 0.0_RKIND_tiempo) + call alloc_and_init(this%valueForTime, BUFSIZE, 0.0_RKIND) + call create_data_file(this%filePathTime, this%path, timeExtension, datFileExtension) + end if + if (any(this%domain%domainType == (/FREQUENCY_DOMAIN, BOTH_DOMAIN/))) then + this%nFreq = this%domain%fnum + allocate (this%frequencySlice(this%domain%fnum)) + call alloc_and_init(this%valueForFreq, this%domain%fnum, (0.0_CKIND, 0.0_CKIND)) + do i = 1, this%nFreq + call init_frequency_slice(this%frequencySlice, this%domain) + end do + this%valueForFreq = (0.0_RKIND, 0.0_RKIND) + + allocate (this%auxExp_E(this%nFreq)) + allocate (this%auxExp_H(this%nFreq)) + do i = 1, this%nFreq + this%auxExp_E(i) = timeInterval*(1.0E0_RKIND, 0.0E0_RKIND)*Exp(mcpi2*this%frequencySlice(i)) !el dt deberia ser algun tipo de promedio + this%auxExp_H(i) = this%auxExp_E(i)*Exp(mcpi2*this%frequencySlice(i)*timeInterval*0.5_RKIND) + end do + call create_data_file(this%filePathFreq, this%path, frequencyExtension, datFileExtension) + end if + + contains + function get_output_path() result(outputPath) + character(len=BUFSIZE) :: probeBoundsExtension, prefixFieldExtension + character(len=BUFSIZE) :: outputPath + probeBoundsExtension = get_coordinates_extension(this%mainCoords, mpidir) + prefixFieldExtension = get_prefix_extension(field, mpidir) + outputPath = & + trim(adjustl(outputTypeExtension))//'_'//trim(adjustl(prefixFieldExtension))//'_'//trim(adjustl(probeBoundsExtension)) + return + end function get_output_path + + end subroutine init_point_probe_output + + subroutine update_point_probe_output(this, step, field) + type(point_probe_output_t), intent(inout) :: this + real(kind=RKIND), pointer, dimension(:, :, :), intent(in) :: field + real(kind=RKIND_tiempo), intent(in) :: step + + integer(kind=SINGLE) :: iter + + if (any(this%domain%domainType == (/TIME_DOMAIN, BOTH_DOMAIN/))) then + this%nTime = this%nTime + 1 + this%timeStep(this%nTime) = step + this%valueForTime(this%nTime) = field(this%mainCoords%x, this%mainCoords%y, this%mainCoords%z) + end if + + if (any(this%domain%domainType == (/FREQUENCY_DOMAIN, BOTH_DOMAIN/))) then + select case(this%component) + case (iEx, iEy, iEz) + do iter = 1, this%nFreq + this%valueForFreq(iter) = & + this%valueForFreq(iter) + field(this%mainCoords%x, this%mainCoords%y, this%mainCoords%z)*(this%auxExp_E(iter)**step) + end do + case (iHx, iHy, iHz) + do iter = 1, this%nFreq + this%valueForFreq(iter) = & + this%valueForFreq(iter) + field(this%mainCoords%x, this%mainCoords%y, this%mainCoords%z)*(this%auxExp_H(iter)**step) + end do + end select + + end if + end subroutine update_point_probe_output + + subroutine flush_point_probe_output(this) + type(point_probe_output_t), intent(inout) :: this + if (any(this%domain%domainType == (/TIME_DOMAIN, BOTH_DOMAIN/))) then + call flush_time_domain(this) + call clear_time_data() + end if + if (any(this%domain%domainType == (/FREQUENCY_DOMAIN, BOTH_DOMAIN/))) then + call flush_frequency_domain(this) + end if + contains + + subroutine flush_time_domain(this) + type(point_probe_output_t), intent(in) :: this + integer :: i + integer :: unit + + if (this%nTime <= 0) then + print *, "No data to write." + return + end if + + open (unit=unit, file=this%filePathTime, status="old", action="write", position="append") + + do i = 1, this%nTime + write (unit, '(F12.6,1X,F12.6)') this%timeStep(i), this%valueForTime(i) + end do + + close (unit) + end subroutine flush_time_domain + + subroutine flush_frequency_domain(this) + type(point_probe_output_t), intent(in) :: this + integer :: i + integer :: unit + + if (.not. allocated(this%frequencySlice) .or. .not. allocated(this%valueForFreq)) then + print *, "Error: arrays not allocated." + return + end if + + if (this%nFreq <= 0) then + print *, "No data to write." + return + end if + open (unit=unit, file=this%filePathFreq, status="replace", action="write") + + do i = 1, this%nFreq + write (unit, '(F12.6,1X,F12.6,1X,F12.6)') this%frequencySlice(i), real(this%valueForFreq(i)), aimag(this%valueForFreq(i)) + end do + + close (unit) + end subroutine flush_frequency_domain + + subroutine clear_time_data() + this%timeStep = 0.0_RKIND_tiempo + this%valueForTime = 0.0_RKIND + + this%nTime = 0 + end subroutine clear_time_data + + end subroutine flush_point_probe_output +end module diff --git a/src_output/wireProbeOutput.F90 b/src_output/wireProbeOutput.F90 new file mode 100644 index 00000000..f8ef97d9 --- /dev/null +++ b/src_output/wireProbeOutput.F90 @@ -0,0 +1,459 @@ +module wireProbeOutput_m + use FDETYPES_m + use utils_m + use report_m + use outputTypes_m + use outputUtils_m + use wiresHolland_constants_m + use HollandWires_m + + implicit none + private + + !=========================== + ! Public interface + !=========================== + public :: init_wire_current_probe_output + public :: init_wire_charge_probe_output + public :: update_wire_current_probe_output + public :: update_wire_charge_probe_output + public :: flush_wire_current_probe_output + public :: flush_wire_charge_probe_output + !=========================== + + !=========================== + ! Private interface + !=========================== + private :: find_current_segment + private :: find_charge_segment + private :: build_output_path + private :: probe_bounds_extension + private :: clear_current_time_data + private :: clear_charge_time_data + private :: update_current_holland + +#ifdef CompileWithBerengerWires + private :: update_current_berenger +#endif + +#ifdef CompileWithSlantedWires + private :: update_current_slanted +#endif + !=========================== + + contains + !====================================================================== + ! INITIALIZATION + !====================================================================== + subroutine init_wire_current_probe_output(this, coordinates, node, field, domain, media, & + outputTypeExtension, mpidir, wiresflavor) + type(wire_current_probe_output_t), intent(out) :: this + type(cell_coordinate_t), intent(in) :: coordinates + type(domain_t), intent(in) :: domain + type(MediaData_t), intent(in) :: media(:) + integer(kind=SINGLE), intent(in) :: node, field, mpidir + character(len=*), intent(in) :: outputTypeExtension, wiresflavor + + this%mainCoords = coordinates + this%component = field + this%domain = domain + this%sign = 1 + + call find_current_segment(this, node, field, media, wiresflavor) + this%path = build_output_path(outputTypeExtension, field, node, mpidir, coordinates) + + call alloc_and_init(this%timeStep, BuffObse, 0.0_RKIND_tiempo) + call create_data_file(this%filePathTime, this%path, timeExtension, datFileExtension) + + end subroutine init_wire_current_probe_output + + + subroutine init_wire_charge_probe_output(this, coordinates, node, field, domain, & + outputTypeExtension, mpidir, wiresflavor) + type(wire_charge_probe_output_t), intent(out) :: this + type(cell_coordinate_t), intent(in) :: coordinates + type(domain_t), intent(in) :: domain + integer(kind=SINGLE), intent(in) :: node, field, mpidir + character(len=*), intent(in) :: outputTypeExtension, wiresflavor + + this%mainCoords = coordinates + this%component = field + this%domain = domain + this%sign = 1 + + call find_charge_segment(this, node, field, wiresflavor) + this%path = build_output_path(outputTypeExtension, field, node, mpidir, coordinates) + + call alloc_and_init(this%timeStep, BuffObse, 0.0_RKIND_tiempo) + call alloc_and_init(this%chargeValue, BuffObse, 0.0_RKIND) + call create_data_file(this%filePathTime, this%path, timeExtension, datFileExtension) + + end subroutine init_wire_charge_probe_output + + !====================================================================== + ! UPDATE + !====================================================================== + subroutine update_wire_current_probe_output(this, step, control, InvEps, InvMu) + type(wire_current_probe_output_t), intent(inout) :: this + real(kind=RKIND_tiempo), intent(in) :: step + type(sim_control_t), intent(in) :: control + real(kind=RKIND), intent(in) :: InvEps(:), InvMu(:) + + this%nTime = this%nTime + 1 + this%timeStep(this%nTime) = step + + select case (trim(control%wiresflavor)) + case ('holland','transition') + call update_current_holland(this, control, InvEps, InvMu) +#ifdef CompileWithBerengerWires + case ('berenger') + call update_current_berenger(this, InvEps, InvMu) +#endif +#ifdef CompileWithSlantedWires + case ('slanted','semistructured') + call update_current_slanted(this) +#endif + end select + end subroutine + + + subroutine update_wire_charge_probe_output(this, step) + type(wire_charge_probe_output_t), intent(inout) :: this + real(kind=RKIND_tiempo), intent(in) :: step + + this%nTime = this%nTime + 1 + this%timeStep(this%nTime) = step + this%chargeValue(this%nTime) = this%segment%ChargeMinus%ChargePresent + end subroutine + + !====================================================================== + ! FLUSH + !====================================================================== + subroutine flush_wire_current_probe_output(this) + type(wire_current_probe_output_t), intent(inout) :: this + integer :: i + integer :: unit + + open(unit, file=this%filePathTime, & + status='old', position='append') + + do i = 1, this%nTime + write(unit, fmt) this%timeStep(i), & + this%currentValues(i)%current, & + this%currentValues(i)%deltaVoltage, & + this%currentValues(i)%plusVoltage, & + this%currentValues(i)%minusVoltage, & + this%currentValues(i)%voltageDiference + end do + close(unit) + + call clear_current_time_data(this) + end subroutine + + + subroutine flush_wire_charge_probe_output(this) + type(wire_charge_probe_output_t), intent(inout) :: this + integer :: i + integer :: unit + + open(unit, file=this%filePathTime, & + status='old', position='append') + + do i = 1, this%nTime + write(unit, fmt) this%timeStep(i), this%chargeValue(i) + end do + close(unit) + + call clear_charge_time_data(this) + end subroutine + + subroutine find_current_segment(this, node, field, media, wiresflavor) + type(wire_current_probe_output_t), intent(inout) :: this + type(MediaData_t), intent(in) :: media(:) + integer(kind=SINGLE), intent(in) :: node, field + character(len=*), intent(in) :: wiresflavor + + type(Thinwires_t), pointer :: Hwireslocal + type(CurrentSegments_t), pointer :: seg +#ifdef CompileWithBerengerWires + type(TWires), pointer :: Hwireslocal_B +#endif +#ifdef CompileWithSlantedWires + type(WiresData), pointer :: Hwireslocal_S +#endif + + integer :: n, iwi, iwj, node2 + logical :: found + character(len=BUFSIZE) :: buff + + found = .false. + this%sign = 1 + + select case (trim(adjustl(wiresflavor))) + case ('holland','transition') + Hwireslocal => GetHwires() + this%segment => Hwireslocal%NullSegment + + do n = 1, Hwireslocal%NumCurrentSegments + seg => Hwireslocal%CurrentSegment(n) + if (seg%origindex == node .and. & + seg%i == iCoord .and. seg%j == jCoord .and. seg%k == kCoord .and. & + seg%tipofield*10 == field) then + found = .true. + this%segment => seg + if (seg%orientadoalreves) this%sign = -1 + exit + end if + end do + +#ifdef CompileWithBerengerWires + case ('berenger') + Hwireslocal_B => GetHwires_Berenger() + do n = 1, Hwireslocal_B%NumSegments + if (Hwireslocal_B%Segments(n)%IndexSegment == node) then + found = .true. + this%segmentBerenger => Hwireslocal_B%Segments(n) + if (Hwireslocal_B%Segments(n)%orientadoalreves) this%sign = -1 + exit + end if + end do +#endif + +#ifdef CompileWithSlantedWires + case ('slanted','semistructured') + Hwireslocal_S => GetHwires_Slanted() + do n = 1, Hwireslocal_S%NumSegments + if (Hwireslocal_S%Segments(n)%ptr%Index == node) then + found = .true. + this%segmentSlanted => Hwireslocal_S%Segments(n)%ptr + exit + end if + end do +#endif + end select + + ! --- multirabo fallback (Holland only) + if (.not. found .and. trim(adjustl(wiresflavor)) /= 'berenger') then + buscarabono: do iwi = 1, Hwireslocal%NumDifferentWires + do iwj = 1, media(Hwireslocal%WireTipoMedio(iwi))%wire(1)%numsegmentos + if (node == media(Hwireslocal%WireTipoMedio(iwi))%wire(1)%segm(iwj)%origindex .and. & + media(Hwireslocal%WireTipoMedio(iwi))%wire(1)%segm(iwj)%multirabo) then + + node2 = media(Hwireslocal%WireTipoMedio(iwi))%wire(1)%segm(iwj)%multiraboDE + do n = 1, Hwireslocal%NumCurrentSegments + seg => Hwireslocal%CurrentSegment(n) + if (seg%origindex == node2) then + found = .true. + this%segment => seg + if (seg%orientadoalreves) this%sign = -1 + exit buscarabono + end if + end do + end if + end do + end do buscarabono + end if + + if (.not. found) then + write(buff,'(a,4i7,a)') 'ERROR: WIRE probe ',node,iCoord,jCoord,kCoord,' DOES NOT EXIST' + call WarnErrReport(buff,.true.) + end if + end subroutine find_current_segment + + subroutine find_charge_segment(this, node, field, wiresflavor) + type(wire_charge_probe_output_t), intent(inout) :: this + integer(kind=SINGLE), intent(in) :: node, field + character(len=*), intent(in) :: wiresflavor + + type(Thinwires_t), pointer :: Hwireslocal + type(CurrentSegments_t), pointer :: seg + integer :: n + logical :: found + character(len=BUFSIZE) :: buff + + found = .false. + this%sign = 1 + + if (trim(adjustl(wiresflavor)) /= 'holland' .and. & + trim(adjustl(wiresflavor)) /= 'transition') then + call WarnErrReport('Charge probes only supported for holland wires', .true.) + return + end if + + Hwireslocal => GetHwires() + + do n = 1, Hwireslocal%NumCurrentSegments + seg => Hwireslocal%CurrentSegment(n) + if (seg%origindex == node .and. & + seg%i == iCoord .and. seg%j == jCoord .and. seg%k == kCoord .and. & + seg%tipofield*10000 == field) then + found = .true. + this%segment => seg + if (seg%orientadoalreves) this%sign = -1 + exit + end if + end do + + if (.not. found) then + write(buff,'(a,4i7,a)') 'ERROR: CHARGE probe ',node,iCoord,jCoord,kCoord,' DOES NOT EXIST' + call WarnErrReport(buff,.true.) + end if + end subroutine find_charge_segment + + function probe_bounds_extension(mpidir, coords) result(ext) + integer(kind=SINGLE), intent(in) :: mpidir + type(cell_coordinate_t), intent(in) :: coords + character(len=BUFSIZE) :: ext + character(len=BUFSIZE) :: ci, cj, ck + + write(ci,'(i7)') coords%x + write(cj,'(i7)') coords%y + write(ck,'(i7)') coords%z + +#if CompileWithMPI + select case (mpidir) + case (3) + ext = trim(adjustl(ci))//'_'//trim(adjustl(cj))//'_'//trim(adjustl(ck)) + case (2) + ext = trim(adjustl(cj))//'_'//trim(adjustl(ck))//'_'//trim(adjustl(ci)) + case (1) + ext = trim(adjustl(ck))//'_'//trim(adjustl(ci))//'_'//trim(adjustl(cj)) + case default + call stoponerror(0,0,'Buggy error in mpidir.') + end select +#else + ext = trim(adjustl(ci))//'_'//trim(adjustl(cj))//'_'//trim(adjustl(ck)) +#endif + end function probe_bounds_extension + + function build_output_path(outExt, field, node, mpidir, coords) result(path) + character(len=*), intent(in) :: outExt + integer(kind=SINGLE), intent(in) :: field, node, mpidir + type(cell_coordinate_t), intent(in) :: coords + character(len=BUFSIZE) :: path + character(len=BUFSIZE) :: nodeStr, fieldExt, boundsExt + + write(nodeStr,'(i7)') node + fieldExt = get_prefix_extension(field, mpidir) + boundsExt = probe_bounds_extension(mpidir, coords) + + path = trim(outExt)//'_'//trim(fieldExt)//'_'// & + trim(boundsExt)//'_s'//trim(adjustl(nodeStr)) + end function build_output_path + + subroutine clear_current_time_data(this) + type(wire_current_probe_output_t), intent(inout) :: this + + this%timeStep = 0.0_RKIND_tiempo + this%currentValues%current = 0.0_RKIND + this%currentValues%deltaVoltage = 0.0_RKIND + this%currentValues%plusVoltage = 0.0_RKIND + this%currentValues%minusVoltage = 0.0_RKIND + this%currentValues%voltageDiference = 0.0_RKIND + this%nTime = 0 + end subroutine clear_current_time_data + + subroutine clear_charge_time_data(this) + type(wire_charge_probe_output_t), intent(inout) :: this + + this%timeStep = 0.0_RKIND_tiempo + this%chargeValue = 0.0_RKIND + this%nTime = 0 + end subroutine clear_charge_time_data + + subroutine update_current_holland(this, control, InvEps, InvMu) + type(wire_current_probe_output_t), intent(inout) :: this + type(sim_control_t), intent(in) :: control + real(kind=RKIND), intent(in) :: InvEps(:), InvMu(:) + + type(CurrentSegments_t), pointer :: seg + + seg => this%segment + + this%currentValues(this%nTime)%current = & + this%sign * seg%currentpast + + this%currentValues(this%nTime)%deltaVoltage = & + - seg%Efield_wire2main * seg%delta + + if (control%wirecrank) then + this%currentValues(this%nTime)%plusVoltage = this%sign * & + (seg%ChargePlus%ChargePresent) * seg%Lind * & + (InvMu(seg%indexmed) * InvEps(seg%indexmed)) + + this%currentValues(this%nTime)%minusVoltage = this%sign * & + (seg%ChargeMinus%ChargePresent) * seg%Lind * & + (InvMu(seg%indexmed) * InvEps(seg%indexmed)) + else + this%currentValues(this%nTime)%plusVoltage = this%sign * & + ((seg%ChargePlus%ChargePresent + seg%ChargePlus%ChargePast) / 2.0_RKIND) * & + seg%Lind * (InvMu(seg%indexmed) * InvEps(seg%indexmed)) + + this%currentValues(this%nTime)%minusVoltage = this%sign * & + ((seg%ChargeMinus%ChargePresent + seg%ChargeMinus%ChargePast) / 2.0_RKIND) * & + seg%Lind * (InvMu(seg%indexmed) * InvEps(seg%indexmed)) + end if + + this%currentValues(this%nTime)%voltageDiference = & + this%currentValues(this%nTime)%plusVoltage - & + this%currentValues(this%nTime)%minusVoltage + end subroutine update_current_holland + +#ifdef CompileWithBerengerWires + subroutine update_current_berenger(this, InvEps, InvMu) + type(wire_current_probe_output_t), intent(inout) :: this + real(kind=RKIND), intent(in) :: InvEps(:), InvMu(:) + + type(TSegment), pointer :: seg + + seg => this%segmentBerenger + + this%currentValues(this%nTime)%current = & + this%sign * seg%currentpast + + this%currentValues(this%nTime)%deltaVoltage = & + - seg%field * seg%dl + + this%currentValues(this%nTime)%plusVoltage = this%sign * & + ((seg%ChargePlus + seg%ChargePlusPast) / 2.0_RKIND) * & + seg%L * (InvMu(seg%imed) * InvEps(seg%imed)) + + this%currentValues(this%nTime)%minusVoltage = this%sign * & + ((seg%ChargeMinus + seg%ChargeMinusPast) / 2.0_RKIND) * & + seg%L * (InvMu(seg%imed) * InvEps(seg%imed)) + + this%currentValues(this%nTime)%voltageDiference = & + this%currentValues(this%nTime)%plusVoltage - & + this%currentValues(this%nTime)%minusVoltage + end subroutine update_current_berenger +#endif + +#ifdef CompileWithSlantedWires + subroutine update_current_slanted(this) + type(wire_current_probe_output_t), intent(inout) :: this + + class(Segment), pointer :: seg + + seg => this%segmentSlanted + + this%currentValues(this%nTime)%current = & + seg%Currentpast + + this%currentValues(this%nTime)%deltaVoltage = & + - seg%field * seg%dl + + this%currentValues(this%nTime)%plusVoltage = & + (seg%Voltage(iPlus)%ptr%Voltage + & + seg%Voltage(iPlus)%ptr%VoltagePast) / 2.0_RKIND + + this%currentValues(this%nTime)%minusVoltage = & + (seg%Voltage(iMinus)%ptr%Voltage + & + seg%Voltage(iMinus)%ptr%VoltagePast) / 2.0_RKIND + + this%currentValues(this%nTime)%voltageDiference = & + this%currentValues(this%nTime)%plusVoltage - & + this%currentValues(this%nTime)%minusVoltage + end subroutine update_current_slanted +#endif + +end module wireProbeOutput_m