Merge pull request #176 from OpenSEMBA/dev

Dev

Author: lmdiazangulo

doc/smbjson.md

@@ -269,14 +269,55 @@ A `material` with `type` `isotropic` represents an isotropic material with const
 
 ```json
 {
-    "name": "teflon"
+    "name": "teflon",
     "id": 1, 
     "type": "isotropic",
     "relativePermittivity": 2.5,
     "electricConducitivity": 1e-6
 } 
 ```
 
+### `lumped`
+
+A `material` with `type` `lumped` represents a lumped circuit `model`, e.g a resistor. Lumped materials can only be assigned to `cell` `elements` with `intervals` describing oriented lines. If multiple cells are assigned to a lumped element, only the first one of them will be treated as a lumped by the solver, the other cells will be treated as a PEC material.
+The specific behavior is described using the `<model>` keyword, described below. `resistor`, `inductor` and `capacitor` are based, with some additions, on the following reference
+
+```
+    Liu, Y., Mittra, R., Su, T., Yang, X., & Yu, W. (2006). Parallel finite-difference time-domain method. Artech.
+    Chapter 3, Section 5.
+```
+
+#### `resistor` model
+
+Defined by:
++ `<resistance>` a positive real number.
++ `[startingTime]` and `[endTime]` are the times in which the resistor will be active. Default to $0.0$ and $1.0$ seconds, respectively. When deactivated, the backgroung material properties will be used, i.e. the edge will be equivalent to an open circuit at low frequencies.
+
+**Example:**
+
+```json
+{
+    "name": "100_ohm_resistor",
+    "id": 1, 
+    "type": "lumped",
+    "model": "resistor",
+    "resistance": 100,
+    "startingTime": 0.0,
+    "endTime": 1.0
+} 
+```
+
+#### `inductor` model
+A series $LR$ circuit, $R$ is optional:
++ `<inductance>` a positive real number
++ `[resistance]` a positive or zero real number. Defaults to $0.0$.
+
+#### `capacitor` model
+A parallel $CR$ circuit:
++ `<capacitance>` a positive real number
++ `<resistance>` a positive real number.
+
+
 ### `multilayeredSurface`
 
 A `multilayeredSurface` must contain the entry `<layers>` which is an array indicating materials which are described in the same way as [isotropic materials](#isotropic) and a `<thickness>`.

src_json_parser/smbjson.F90

@@ -547,7 +547,8 @@ subroutine fillDielectricsOfCellType(res, cellType)
          integer :: i, j
          integer :: nCs, nDielectrics
 
-         mAs = this%getMaterialAssociations([J_MAT_TYPE_ISOTROPIC])
+         mAs = this%getMaterialAssociations( &
+            [J_MAT_TYPE_ISOTROPIC, J_MAT_TYPE_LUMPED])
          if (size(mAs) == 0) then
             allocate(res(0))
             return
@@ -567,11 +568,19 @@ subroutine fillDielectricsOfCellType(res, cellType)
          if (nDielectrics == 0) return
 
          j = 0
+         mAs = this%getMaterialAssociations([J_MAT_TYPE_ISOTROPIC])
          do i = 1, size(mAs)       
             if (.not. containsCellRegionsWithType(mAs(i), cellType)) cycle
             j = j + 1
             res(j) = readDielectric(mAs(i), cellType)
          end do
+
+         mAs = this%getMaterialAssociations([J_MAT_TYPE_LUMPED])
+         do i = 1, size(mAs)
+            if (.not. containsCellRegionsWithType(mAs(i), cellType)) cycle
+            j = j + 1
+            res(j) = readLumped(mAs(i), cellType)
+         end do
       end subroutine
 
       function readDielectric(mA, cellType) result(res)
@@ -597,6 +606,78 @@ function readDielectric(mA, cellType) result(res)
 
       end function
 
+
+      function readLumped(mA, cellType) result(res)
+         type(materialAssociation_t), intent(in) :: mA
+         integer, intent(in) :: cellType
+         type(Dielectric_t) :: res
+         type(cell_region_t) :: cR
+         type (coords), dimension(:), allocatable :: coords
+         type(json_value_ptr) :: matPtr
+         integer :: e, j
+         character(len=:), allocatable :: model
+         logical :: found
+
+         allocate(res%c1P(0))
+         res%n_c1p = 0
+         call this%matAssToCoords(res%c2p, mA, cellType)
+         res%n_c2p = size(res%c2p)
+         
+         matPtr = this%matTable%getId(mA%materialId)
+         
+         ! Get the model type
+         model = this%getStrAt(matPtr%p, J_MAT_LUMPED_MODEL, found)
+         if (.not. found) then
+            write(error_unit, *) "ERROR reading lumped material: ", mA%materialId, " model not found."
+            stop
+         end if
+         
+         ! Not really needed for resistor, inductor, or capacitor. 
+         ! But avoids error in lumped initialization.
+         res%orient = 1
+         res%DiodOri = 1
+
+         res%eps = EPSILON_VACUUM
+         res%mu = MU_VACUUM
+
+         ! Handle resistor model
+         select case (model)
+          case (J_MAT_LUMPED_MODEL_RESISTOR)
+            res%resistor = .true.
+            res%R = this%getRealAt(matPtr%p, J_MAT_LUMPED_RESISTANCE, found)
+            if (.not. found) then
+               write(error_unit, *) "ERROR reading lumped material: ", mA%materialId, " resistance not found."
+               stop
+            end if
+            res%Rtime_on = this%getRealAt(matPtr%p, J_MAT_LUMPED_STARTING_TIME, default=0.0)
+            res%Rtime_off = this%getRealAt(matPtr%p, J_MAT_LUMPED_END_TIME, default=1.0)
+          case (J_MAT_LUMPED_MODEL_INDUCTOR)
+            res%inductor = .true.
+            res%L = this%getRealAt(matPtr%p, J_MAT_LUMPED_INDUCTANCE, found)
+            if (.not. found) then
+               write(error_unit, *) "ERROR reading lumped material: ", mA%materialId, " inductance not found."
+               stop
+            end if
+            res%R = this%getRealAt(matPtr%p, J_MAT_LUMPED_RESISTANCE, default=0.0)
+          case (J_MAT_LUMPED_MODEL_CAPACITOR)
+            res%capacitor = .true.
+            res%C = this%getRealAt(matPtr%p, J_MAT_LUMPED_CAPACITANCE, found)
+            if (.not. found) then
+               write(error_unit, *) "ERROR reading lumped material: ", mA%materialId, " capacitance not found."
+               stop
+            end if
+            res%R = this%getRealAt(matPtr%p, J_MAT_LUMPED_RESISTANCE, found)
+            if (.not. found) then
+               write(error_unit, *) "ERROR reading lumped material: ", mA%materialId, " resistance not found."
+               stop
+            end if
+          case default
+            write(error_unit, *) "ERROR reading lumped material: ", mA%materialId, " invalid model."
+            stop
+          end select
+
+      end function
+
       logical function containsCellRegionsWithType(mA, cellType)
          integer, intent(in) :: cellType
          type(materialAssociation_t), intent(in) :: mA

src_json_parser/smbjson_labels.F90

@@ -30,6 +30,7 @@ module smbjson_labels_mod
    character (len=*), parameter :: J_MAT_TYPE_PEC = "pec"
    character (len=*), parameter :: J_MAT_TYPE_PMC = "pmc"
    character (len=*), parameter :: J_MAT_TYPE_ISOTROPIC = "isotropic"
+   character (len=*), parameter :: J_MAT_TYPE_LUMPED = "lumped"
    character (len=*), parameter :: J_MAT_TYPE_MULTILAYERED_SURFACE = "multilayeredSurface"
    character (len=*), parameter :: J_MAT_TYPE_SLOT = "thinSlot"
    character (len=*), parameter :: J_MAT_TYPE_WIRE = "wire"
@@ -46,7 +47,17 @@ module smbjson_labels_mod
    character (len=*), parameter :: J_MAT_WIRE_DIELECTRIC_RADIUS = "radius"
    character (len=*), parameter :: J_MAT_WIRE_DIELECTRIC_PERMITTIVITY = "relativePermittivity"
    character (len=*), parameter :: J_MAT_WIRE_PASS = "isPassthrough"
-   
+
+   character (len=*), parameter :: J_MAT_LUMPED_MODEL = "model"
+   character (len=*), parameter :: J_MAT_LUMPED_MODEL_RESISTOR = "resistor"
+   character (len=*), parameter :: J_MAT_LUMPED_MODEL_INDUCTOR = "inductor"
+   character (len=*), parameter :: J_MAT_LUMPED_MODEL_CAPACITOR = "capacitor"
+   character (len=*), parameter :: J_MAT_LUMPED_RESISTANCE = "resistance"
+   character (len=*), parameter :: J_MAT_LUMPED_STARTING_TIME = "startingTime"
+   character (len=*), parameter :: J_MAT_LUMPED_END_TIME = "endTime"
+   character (len=*), parameter :: J_MAT_LUMPED_INDUCTANCE = "inductance"
+   character (len=*), parameter :: J_MAT_LUMPED_CAPACITANCE = "capacitance" 
+      
    character (len=*), parameter :: J_MAT_TERM_TERMINATIONS = "terminations"
    character (len=*), parameter :: J_MAT_TERM_TYPE_OPEN = "open"
    character (len=*), parameter :: J_MAT_TERM_TYPE_SHORT = "short"

src_main_pub/lumped.F90

@@ -1,12 +1,8 @@
-
-    
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Module Lumped
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!!!!!17/08/15 update!!!!!!!!!!!
 !!!Elimino el tratamiento de los campos magneticos de Lumped para programar un multiLumped 
 !!!solo teniendo en cuenta los parametros efectivos y sin actualizar los magneticos.
-!!!Mangento en el fichero Lumped_pre170815_noupdateababienH.F90 la version antigua
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 module Lumped
@@ -20,17 +16,14 @@ module Lumped
 
    implicit none
 
-   !!!!!valiables locales
    type (LumpedElem_t), save, target   ::  LumpElem
 
-!!!variables globales del modulo
    REAL (KIND=RKIND), save           ::  eps0,mu0,zvac,cluz
 
    private
-
 !!!
 !!!
-   public LumpedElem_t,Nodes_t !el tipo es publico
+   public LumpedElem_t,Nodes_t 
    public AdvanceLumpedE,InitLumped,DestroyLumped,StoreFieldsLumpeds,calc_lumpedconstants,Getlumped
 
 contains
@@ -192,7 +185,7 @@ subroutine InitLumped(sgg,sggMiEx,sggMiEy,sggMiEz,Ex,Ey,Ez,Hx,Hy,Hz,&
             lumped_ => LumpElem%Nodes(conta)
             lumped_%EfieldPrevPrev=0.0_RKIND
             lumped_%EfieldPrev    =0.0_RKIND
-            lumped_%Jcur          =0.0_RKIND !olvide inicializar jcur 071118
+            lumped_%Jcur          =0.0_RKIND
          end do
 #ifdef CompileWithStochastic
          if (stochastic) then
@@ -207,13 +200,13 @@ subroutine InitLumped(sgg,sggMiEx,sggMiEy,sggMiEz,Ex,Ey,Ez,Hx,Hy,Hz,&
       else  
         do conta=1,LumpElem%numnodes
             lumped_ => LumpElem%Nodes(conta)
-            read(14) lumped_%EfieldPrevPrev,lumped_%EfieldPrev,lumped_%Jcur !olvide almacenar jcur 071118
+            read(14) lumped_%EfieldPrevPrev,lumped_%EfieldPrev,lumped_%Jcur 
         end do  
 #ifdef CompileWithStochastic
          if (stochastic) then
              do conta=1,LumpElem%numnodes
                 lumped_ => LumpElem%Nodes(conta)
-                read(14) lumped_%EfieldPrevPrev_for_devia,lumped_%EfieldPrev_for_devia,lumped_%Jcur_for_devia !olvide almacenar jcur 071118
+                read(14) lumped_%EfieldPrevPrev_for_devia,lumped_%EfieldPrev_for_devia,lumped_%Jcur_for_devia
              end do
          endif
 #endif
@@ -238,7 +231,8 @@ subroutine AdvanceLumpedE(sgg,timestep,simu_devia,stochastic)
          if (sgg%med(lumped_%jmed)%lumped(1)%inductor) then
              lumped_%EfieldPrevPrev = lumped_%EfieldPrev
              lumped_%EfieldPrev     = lumped_%Efield
-             lumped_%Jcur = lumped_%Jcur + lumped_%sigmaEffResistInduct * (lumped_%EfieldPrev + lumped_%EfieldPrevPrev)
+             !!! The evolution of the current must be modified by a coefficient different from that shown in Mittra pag65.
+             lumped_%Jcur = lumped_%currentCoeff * lumped_%Jcur + lumped_%sigmaEffResistInduct * (lumped_%EfieldPrev + lumped_%EfieldPrevPrev)
          else
              lumped_%Jcur=0.0_RKIND
          endif
@@ -250,14 +244,14 @@ subroutine AdvanceLumpedE(sgg,timestep,simu_devia,stochastic)
          else !debe entrar aqui si es un resistor, inductor o capacitor
             if (sgg%med(lumped_%jmed)%lumped(1)%resistor) then
                 if ((timestep*sgg%dt >= sgg%Med(lumped_%jmed)%Lumped(1)%Rtime_on).and.(timestep*sgg%dt <= sgg%Med(lumped_%jmed)%Lumped(1)%Rtime_off)) then
-                   lumped_%Efield = lumped_%G1 * lumped_%Efield +  (lumped_%G2a *(lumped_%Ha_Plus   - lumped_%Ha_Minu    ) - lumped_%G2b *(lumped_%Hb_Plus     - lumped_%Hb_Minu  ) ) + &
+                   lumped_%Efield = lumped_%G1 * lumped_%Efield +  (lumped_%G2a *(lumped_%Ha_Plus   - lumped_%Ha_Minu    ) - lumped_%G2b *(lumped_%Hb_Plus     - lumped_%Hb_Minu  ) ) - &
                                     lumped_%GJ * lumped_%Jcur
                 else
                    lumped_%Efield = lumped_%G1_usual * lumped_%Efield +(lumped_%G2a_usual *(lumped_%Ha_Plus - lumped_%Ha_Minu) - &
                                                                         lumped_%G2b_usual *(lumped_%Hb_Plus - lumped_%Hb_Minu))
                 endif
             else !inductor o capacitor
-                lumped_%Efield = lumped_%G1 * lumped_%Efield +  (lumped_%G2a *(lumped_%Ha_Plus   - lumped_%Ha_Minu    ) - lumped_%G2b *(lumped_%Hb_Plus     - lumped_%Hb_Minu  ) ) + &
+                lumped_%Efield = lumped_%G1 * lumped_%Efield +  (lumped_%G2a *(lumped_%Ha_Plus   - lumped_%Ha_Minu    ) - lumped_%G2b *(lumped_%Hb_Plus     - lumped_%Hb_Minu  ) ) - &
                                  lumped_%GJ * lumped_%Jcur
             endif
          endif                     
@@ -282,7 +276,7 @@ subroutine calc_lumpedconstants(sgg,eps00,mu00)
       real (kind=RKIND) :: epsilon,sigma,g1,g2,Resist,Induct,Capaci,sigmaeff,epsiloneff,DiodB,DiodIsat
       real (kind=RKIND) :: g1_usual,g2_usual
       real (kind=RKIND) :: epsilonEffCapac    ,sigmaEffResistInduct,sigmaEffResist   ,sigmaEffResistCapac ,sigmaEffResistDiode, &
-                           alignedDeltaE,transversalDeltaHa,transversalDeltaHb 
+                           alignedDeltaE,transversalDeltaHa,transversalDeltaHb,  currentCoeff
 
       character(len=BUFSIZE) :: buff
 !
@@ -307,13 +301,13 @@ subroutine calc_lumpedconstants(sgg,eps00,mu00)
             alignedDeltaE      =lumped_%alignedDeltaE      
             transversalDeltaHa =lumped_%transversalDeltaHa 
             transversalDeltaHb =lumped_%transversalDeltaHb     
-               
 !
             epsilonEffCapac      = alignedDeltaE * Capaci / (            transversalDeltaHa * transversalDeltaHb)
             sigmaEffResistInduct = alignedDeltaE * sgg%dt / (2.0_RKIND * transversalDeltaHa * transversalDeltaHb * (Induct + Resist * sgg%dt /2.0_RKIND))
             sigmaEffResist       = alignedDeltaE          / (   Resist * transversalDeltaHa * transversalDeltaHb)
             sigmaEffResistCapac  = sigmaEffResist
             sigmaEffResistDiode  = sigmaEffResist
+            currentCoeff         = (Induct - Resist * sgg%dt /2.0_RKIND) / (Induct + Resist * sgg%dt /2.0_RKIND)
 
 !!!! Mittra pag65   Parallel Finite-Difference Time-Domain Method
             if (sgg%med(jmed)%lumped(1)%resistor) then
@@ -348,7 +342,6 @@ subroutine calc_lumpedconstants(sgg,eps00,mu00)
                (epsilonEff/sgg%dt   + SigmaEff/2.0_RKIND  ) 
 
 
-    !!!!repensar si es necesario 27/08/15
             !!!lo he comentado a 050122 por consistencia con stochastic
             !!if (g1 < 0.0_RKIND) then !exponential time stepping
             !!    g1=exp(- SigmaEff * sgg%dt / (epsilonEff ))
@@ -357,13 +350,13 @@ subroutine calc_lumpedconstants(sgg,eps00,mu00)
             lumped_%g1=g1 
             lumped_%G2a= G2 / transversalDeltaHa 
             lumped_%G2b= G2 / transversalDeltaHb
-            lumped_%GJ = G2 !!!!Mittra pag65
+            !!! The current derivation on page 65 of Mittra is incorrect and leads to saturation. 
+            !!! The correct coefficient is the one shown here.
+            lumped_%GJ = G2 * (1 + currentCoeff) / 2 !!!!Mittra pag65
             lumped_%sigmaEffResistInduct = sigmaEffResistInduct
+            lumped_%currentCoeff = currentCoeff
 
-            
-            
             !!!!usual para resistencia que se encienden/apagan 200319
-
             G1_usual=(1.0_RKIND  - Sigma * sgg%dt / (2.0_RKIND * epsilon) ) / &
                 (1.0_RKIND  + Sigma * sgg%dt / (2.0_RKIND * epsilon) ) 
             G2_usual=  sgg%dt / epsilon                        / &
@@ -377,7 +370,6 @@ subroutine calc_lumpedconstants(sgg,eps00,mu00)
             lumped_%G2a_usual= G2_usual / transversalDeltaHa 
             lumped_%G2b_usual= G2_usual / transversalDeltaHb
 
-                
            !!!only for diodes 
             if (orient>0.0) then
                 lumped_%diodeB    = lumped_%diodeB * alignedDeltaE / 2.0_RKIND

src_main_pub/lumped_types.F90

@@ -5,7 +5,7 @@ module lumped_vars
    use fdetypes   
 
    TYPE, public  ::  Nodes_t
-      REAL (KIND=RKIND) :: EfieldPrev,EfieldPrevPrev,Jcur,sigmaEffResistInduct,alignedDeltaE ,transversalDeltaHa ,transversalDeltaHb 
+      REAL (KIND=RKIND) :: EfieldPrev,EfieldPrevPrev,Jcur,sigmaEffResistInduct,alignedDeltaE ,transversalDeltaHa ,transversalDeltaHb, currentCoeff 
       REAL (KIND=RKIND) :: diodepreA,diodeB
       REAL (KIND=RKIND), pointer ::  Efield,Ha_Plus,Ha_Minu,Hb_Plus,Hb_Minu
       REAL (KIND=RKIND)          ::  g1

src_main_pub/semba_fdtd.F90

@@ -722,8 +722,8 @@ PROGRAM SEMBA_FDTD_launcher
             write(thefileno,'(a)') '# (  8.0 ,  8.0 ) '//trim(adjustl('Thin wire segments colliding with structure                             (Line)'))
             write(thefileno,'(a)') '# (  8.5 ,  8.5 ) '//trim(adjustl('Soft/Hard Nodal CURRENT/FIELD ELECTRIC DENSITY SOURCE                   (Line)'))
             write(thefileno,'(a)') '# (  9.0 ,  9.0 ) '//trim(adjustl('Soft/Hard Nodal CURRENT/FIELD MAGNETIC DENSITY SOURCE                   (Line)'))
-            write(thefileno,'(a)') '# (   10 ,   11 ) '//trim(adjustl('LeftEnd/RightEnd/Ending  segment                                                 (Wire)'))
-            write(thefileno,'(a)') '# (   20 ,   20 ) '//trim(adjustl('Intermediate segment +number_holland_parallel or +number_berenger       (Wire) '))
+            write(thefileno,'(a)') '# (   10 ,   11 ) '//trim(adjustl('LeftEnd/RightEnd/Ending wire segment                                                 (Wire)'))
+            write(thefileno,'(a)') '# (   20 ,   20 ) '//trim(adjustl('Intermediate wire segment +number_holland_parallel or +number_berenger       (Wire) '))
             write(thefileno,'(a)') '# (  400 ,  499 ) '//trim(adjustl('Thin slot (+indexmedium)                                                (Surface)'))
             write(thefileno,'(a)') '# ( -0.5 , -0.5 ) '//trim(adjustl('Other types of media                                                    (Line)'))
             write(thefileno,'(a)') '# ( -1.0 , -1.0 ) '//trim(adjustl('Other types of media                                                    (Surface)'))

src_pyWrapper/pyWrapper.py

@@ -390,7 +390,7 @@ def getCurrentVTKMap(self):
     def getMaterialProperties(self, materialName):
         if 'materials' in self._input:
             for idx, element in enumerate(self._input['materials']):
-                if element["name"] == materialName:
+                if element.get("name") == materialName:
                     return self._input['materials'][idx]