349 mtln probe flush

src_main_pub/observation.F90

@@ -114,7 +114,7 @@ module Observa_m
   !Required at least in tests
   public fieldo
 #ifdef CompileWithMTLN
-  public FlushMTLNObservationFiles
+  public InitObservationMTLN, UpdateObservationMTLN, CloseObservationFilesMTLN
 #endif
 contains
 
@@ -2182,6 +2182,10 @@ subroutine InitObservation(sgg, media, tag_numbers, &
             !!!!        end if !del time domain !NO ES PRECISO 25/02/14
         end do !del ii=1,numberrequest
 
+#ifdef CompileWithMTLN
+        call InitObservationMTLN(control%nEntradaRoot)
+#endif
+
         write (19, '(a)') '!END '
         close (19)
 
@@ -2636,7 +2640,9 @@ subroutine CloseObservationFiles(sgg, layoutnumber, num_procs, singlefilewrite,
           end if !no hay singlewrite para sondas freqdomain
         end do
       end do
-
+#ifdef CompileWithMTLN
+      call CloseObservationFilesMTLN()
+#endif
       return
     end subroutine CloseObservationFiles
 
@@ -3304,6 +3310,9 @@ subroutine UpdateObservation(sgg, media, tag_numbers, &
           end if !del nothing
         end do loop_obser
       end do
+#ifdef CompileWithMTLN
+      call UpdateObservationMTLN(nTime)
+#endif
       !---------------------------> acaba UpdateObservation <-----------------------------------------
       return
 
@@ -4074,18 +4083,13 @@ subroutine FlushObservationFiles(sgg,nInit,FinalInstant,layoutnumber,num_procs,
     end subroutine FlushObservationFiles
 
 #ifdef CompileWithMTLN
-    subroutine FlushMTLNObservationFiles(nEntradaRoot, mtlnProblem)
+    subroutine InitObservationMTLN(nEntradaRoot)
       character(len=*), intent(in) :: nEntradaRoot
-      logical, intent(in) :: mtlnProblem
       type(mtln_solver_t), pointer :: mtln_solver
-      integer :: i, j, k, n
-      integer :: unit
-      character(len=bufsize) :: temp
+      integer :: i, j, k, unit
       character(len=bufsize) :: path
+      character(len=bufsize) :: temp
       character(len=:), allocatable :: buffer
-#ifdef CompileWithMPI
-      integer(kind=4) :: ierr
-#endif
 
       mtln_solver => GetSolverPtr()
       if (.not. allocated(mtln_solver%bundles)) return
@@ -4095,31 +4099,65 @@ subroutine FlushMTLNObservationFiles(nEntradaRoot, mtlnProblem)
 #ifdef CompileWithMPI
           if (.not. mtln_solver%bundles(i)%probes(j)%in_layer) cycle
 #endif          
+          mtln_solver%bundles(i)%probes(j)%unit = unit
           path = trim(trim(nEntradaRoot)//"_"//trim(mtln_solver%bundles(i)%probes(j)%name)//".dat")
           open (unit=unit, file=trim(path))
           write (*, *) 'name: ', trim(mtln_solver%bundles(i)%probes(j)%name)
           buffer = "time"
-
           do k = 1, size(mtln_solver%bundles(i)%probes(j)%val, 2)
             write (temp, *) k
             buffer = buffer//" "//"conductor_"//trim(adjustl(temp))
           end do
-          write (unit, *) trim(buffer)
-          do k = 1, size(mtln_solver%bundles(i)%probes(j)%val, 1)
-            buffer = ""
-            write (temp, *) mtln_solver%bundles(i)%probes(j)%t(k)
-            buffer = buffer//trim(temp)
-            do n = 1, size(mtln_solver%bundles(i)%probes(j)%val, 2)
-              write (temp, *) mtln_solver%bundles(i)%probes(j)%val(k, n)
-              buffer = buffer//" "//trim(temp)
-            end do
-            write (unit, '(a)') trim(buffer)
-          end do
-          close (unit)
+          write (unit, '(a)') trim(buffer)
           unit = unit + 1
         end do
       end do
+    end subroutine
 
+    subroutine CloseObservationFilesMTLN()
+      type(mtln_solver_t), pointer :: mtln_solver
+      integer :: i, j, k
+      mtln_solver => GetSolverPtr()
+      if (.not. allocated(mtln_solver%bundles)) return
+      do i = 1, size(mtln_solver%bundles)
+        do j = 1, size(mtln_solver%bundles(i)%probes)
+#ifdef CompileWithMPI
+          if (.not. mtln_solver%bundles(i)%probes(j)%in_layer) cycle
+#endif          
+          close(mtln_solver%bundles(i)%probes(k)%unit)
+        end do
+      end do
+    end subroutine
+
+    subroutine UpdateObservationMTLN(step)
+      integer, intent(in) :: step
+      type(mtln_solver_t), pointer :: mtln_solver
+      integer :: i, j, k, n
+      integer :: unit
+      character(len=bufsize) :: temp
+      character(len=bufsize) :: path
+      character(len=:), allocatable :: buffer
+#ifdef CompileWithMPI
+      integer(kind=4) :: ierr
+#endif
+
+      mtln_solver => GetSolverPtr()
+      if (.not. allocated(mtln_solver%bundles)) return
+      do i = 1, size(mtln_solver%bundles)
+        do j = 1, size(mtln_solver%bundles(i)%probes)
+#ifdef CompileWithMPI
+          if (.not. mtln_solver%bundles(i)%probes(j)%in_layer) cycle
+#endif          
+          buffer = ""
+          write(temp, *) mtln_solver%bundles(i)%probes(j)%t(step+1)
+          buffer = buffer//trim(temp)
+          do n = 1, size(mtln_solver%bundles(i)%probes(j)%val, 2)
+            write (temp, *) mtln_solver%bundles(i)%probes(j)%val(step+1, n)
+            buffer = buffer//" "//trim(temp)
+          end do
+          write (mtln_solver%bundles(i)%probes(j)%unit, '(a)') trim(buffer)
+        end do
+      end do
     end subroutine
 #endif
 

src_main_pub/timestepping.F90

@@ -73,7 +73,6 @@ module Solver_m
    use P_rescale
 #endif              
 #ifdef CompileWithMTLN
-   ! use mtln_solver_mod, mtln_solver_t => mtln_t
    use mtln_types_m, only: mtln_t
    use Wire_bundles_mtln_m
 #endif
@@ -280,10 +279,9 @@ subroutine launch_mtln_simulation(this, mtln_parsed, nEntradaRoot, layoutnumber)
       type(mtln_t) :: mtln_parsed
       character(len=*), intent(in) :: nEntradaRoot
       integer(kind=4), intent(in) :: layoutnumber
-
-      call solveMTLNProblem(mtln_parsed)
+      
+      call solveMTLNProblem(mtln_parsed, nEntradaRoot)
       call reportSimulationEnd(layoutnumber)
-      call FlushMTLNObservationFiles(nEntradaRoot, mtlnProblem = .true.)
    end subroutine
 #endif
 
@@ -1466,7 +1464,6 @@ subroutine initializeObservation()
          call InitObservation (this%sgg,this%media,this%tag_numbers, &
                                  this%thereAre%Observation,this%thereAre%wires,this%thereAre%FarFields,this%initialtimestep,this%lastexecutedtime, &
                                  this%sinPML_fullsize,this%eps0,this%mu0,this%bounds, this%control)
-
          l_auxinput=this%thereAre%Observation.or.this%thereAre%FarFields
          l_auxoutput=l_auxinput
 
@@ -2715,9 +2712,6 @@ subroutine solver_end(this)
       if (this%thereAre%Observation) then
          call FlushObservationFiles(this%sgg,this%ini_save, this%n,this%control%layoutnumber, this%control%num_procs, dxe, dye, dze, dxh, dyh, dzh,this%bounds,this%control%singlefilewrite,this%control%facesNF2FF,.TRUE.)
          call CloseObservationFiles(this%sgg,this%control%layoutnumber,this%control%num_procs,this%control%singlefilewrite,this%initialtimestep,this%lastexecutedtime,this%control%resume) !dump the remaining to disk
-#ifdef CompileWithMTLN      
-         call FlushMTLNObservationFiles(this%control%nentradaroot, mtlnProblem = .false.)
-#endif
       end if
 
       if (this%thereAre%FarFields) then

src_mtln/mtln_solver.F90

@@ -36,7 +36,9 @@ module mtln_solver_m
 
         procedure :: runUntil
         procedure :: run => mtln_run
-
+        procedure :: initObservation => mtln_initObservation
+        procedure :: updateObservation => mtln_updateObservation
+        procedure :: closeObservation => mtln_closeObservation
     end type mtln_t
 
 
@@ -270,8 +272,9 @@ subroutine runUntil(this, final_time)
             call this%advanceBundlesVoltage()
             call this%advanceNWVoltage()
             call this%advanceBundlesCurrent()
-            call this%advanceTime()
             call this%updateProbes()
+            call this%advanceTime()
+            call this%updateObservation(i)
         end do
 
     end subroutine
@@ -281,15 +284,89 @@ subroutine mtln_run(this)
         real(kind=RKIND_TIEMPO) :: time
         integer :: i
 
-        call this%updatePULTerms()
         do i = 0, this%getTimeRange(this%final_time)
             call this%advanceBundlesVoltage()
             call this%advanceNWVoltage()
             call this%advanceBundlesCurrent()
-            call this%advanceTime()
             call this%updateProbes()
+            call this%advanceTime()
+            call this%updateObservation(i)
         end do
 
     end subroutine
 
+    subroutine mtln_initObservation(this, nEntradaRoot)
+        class(mtln_t) :: this
+        character(len=*), intent(in) :: nEntradaRoot
+        integer :: i, j, k, unit
+        character(len=bufsize) :: path
+        character(len=bufsize) :: temp
+        character(len=:), allocatable :: buffer
+
+        if (.not. allocated(this%bundles)) return
+        unit = 2000
+        do i = 1, size(this%bundles)
+            do j = 1, size(this%bundles(i)%probes)
+#ifdef CompileWithMPI
+                if (.not. this%bundles(i)%probes(j)%in_layer) cycle
+#endif          
+                this%bundles(i)%probes(j)%unit = unit
+                path = trim(trim(nEntradaRoot)//"_"//trim(this%bundles(i)%probes(j)%name)//".dat")
+                open (unit=unit, file=trim(path))
+                write (*, *) 'name: ', trim(this%bundles(i)%probes(j)%name)
+                buffer = "time"
+                do k = 1, size(this%bundles(i)%probes(j)%val, 2)
+                    write (temp, *) k
+                    buffer = buffer//" "//"conductor_"//trim(adjustl(temp))
+                end do
+                write (unit, '(a)') trim(buffer)
+                unit = unit + 1
+            end do
+        end do
+    end subroutine
+
+    subroutine mtln_updateObservation(this, step)
+        class(mtln_t) :: this
+        integer, intent(in) :: step
+        integer :: i, j, k, n
+        integer :: unit
+        character(len=bufsize) :: temp
+        character(len=bufsize) :: path
+        character(len=:), allocatable :: buffer
+#ifdef CompileWithMPI
+        integer(kind=4) :: ierr
+#endif
+        if (.not. allocated(this%bundles)) return
+        do i = 1, size(this%bundles)
+            do j = 1, size(this%bundles(i)%probes)
+#ifdef CompileWithMPI
+                if (.not. this%bundles(i)%probes(j)%in_layer) cycle
+#endif          
+                buffer = ""
+                write(temp, *) this%bundles(i)%probes(j)%t(step+1)
+                buffer = buffer//trim(temp)
+                do n = 1, size(this%bundles(i)%probes(j)%val, 2)
+                    write (temp, *) this%bundles(i)%probes(j)%val(step+1, n)
+                    buffer = buffer//" "//trim(temp)
+                end do
+                write (this%bundles(i)%probes(j)%unit, '(a)') trim(buffer)
+            end do
+        end do
+    end subroutine
+
+    subroutine mtln_closeObservation(this)
+        class(mtln_t) :: this
+        integer :: i, j, k
+        if (.not. allocated(this%bundles)) return
+        do i = 1, size(this%bundles)
+            do j = 1, size(this%bundles(i)%probes)
+#ifdef CompileWithMPI
+                if (.not. this%bundles(i)%probes(j)%in_layer) cycle
+#endif          
+                close(this%bundles(i)%probes(k)%unit)
+            end do
+      end do
+    end subroutine
+
+
 end module mtln_solver_m

src_mtln/probes.F90

@@ -15,7 +15,7 @@ module probes_m
         real(kind=RKIND), allocatable, dimension(:) :: t
         real(kind=RKIND), allocatable, dimension(:,:) :: val
         real(kind=RKIND_TIEMPO) :: dt
-        integer :: index, current_frame
+        integer :: index, current_frame, unit
         character(len=:), allocatable :: name
         logical :: in_layer = .true.
 

src_wires_pub/wires_mtln.F90

@@ -197,16 +197,20 @@ function GetSolverPtr() result(res)
       return
    end function
 
-   subroutine solveMTLNProblem(mtln_parsed)
+   subroutine solveMTLNProblem(mtln_parsed, nEntradaRoot)
       type(mtln_t) :: mtln_parsed
+      character(len=*), intent(in) :: nEntradaRoot
       mtln_solver = mtlnCtor(mtln_parsed)
+      call mtln_solver%updatePULTerms()
+      call mtln_solver%initObservation(nEntradaRoot)
       call mtln_solver%run()
+      call mtln_solver%closeObservation()
    end subroutine
 
    subroutine reportSimulationEnd(layoutnumber)
       character(len=bufsize) :: dubuf
       integer(kind=4), intent(in) :: layoutnumber
-      write(dubuf, *) 'MTLN simulation finished. Init flusing probe data to output files'
+      write(dubuf, *) 'MTLN simulation finished.'
       call print11(layoutnumber,dubuf)
 
    end subroutine

test/pyWrapper/test_full_system.py

@@ -1528,4 +1528,3 @@ def test_bulk_current_outputs(tmp_path):
     assert probeBulkZPlane.direction == 'z'
     assert probeBulkYPoint.direction == 'y'
     assert probeBulkZVolume.direction == 'z'
-

testData/cases/paul/paul_prepost.py

@@ -0,0 +1,40 @@
+# %%
+import numpy as np
+from numpy.fft import *
+import matplotlib.pyplot as plt
+import scipy.constants
+
+import sys, os
+sys.path.append(os.path.join(os.path.dirname(__file__), '../../../', 'src_pyWrapper'))
+SEMBA_EXE = '../../../build/bin/semba-fdtd'
+OUTPUTS_FOLDER = '../../outputs/'
+CASES_FOLDER = '../../cases/'
+from pyWrapper import *
+#%%
+fn = CASES_FOLDER + 'paul/paul_8_6_triangle.fdtd.json'
+
+solver = FDTD(input_filename=fn, path_to_exe=SEMBA_EXE)
+solver.run()
+
+#%%
+p_expected = Probe(
+    OUTPUTS_FOLDER+'paul_8_6_triangle.fdtd_start_voltage_wire_V_5_5_1.dat')
+
+probe_voltage = solver.getSolvedProbeFilenames("start_voltage")[0]
+probe_current = solver.getSolvedProbeFilenames("end_current")[0]
+probe_files = [probe_voltage, probe_current]
+p_solved = Probe(probe_files[0])
+
+#%%
+plt.plot(p_expected['time'], p_expected['voltage_0'], label = 'reference')
+plt.plot(p_solved['time'], p_solved['voltage_0'], '--',label = 'solved')
+plt.xlabel('time [s]')
+plt.ylabel('V [V]')
+plt.legend()
+#%%
+solved = np.interp(p_expected['time'].to_numpy(), 
+                    p_solved['time'].to_numpy(), 
+                    p_solved['voltage_0'].to_numpy())
+assert np.corrcoef(solved, p_expected['voltage_0'])[0,1] > 0.999
+
+# %%