Rename of procedures in boundaries

doc/user-manual/fpakc_UserManual.tex

@@ -581,6 +581,19 @@ make
                     
                   \textbf{crossSection}: Character.
                     Complete path to the cross-section data for the ionization process.
+                    
+              \item \textbf{quasiNeutrality}: Maintains equal density between the species defined by the \textbf{incident} parameter and the others in the case by modifying the reflection coefficient and reflecting incident particles of the species.
+              The required input is:
+                
+                \textbf{incident}: Character.
+                  Name of the incident species that will be reflected.
+                    
+              \item \textbf{outflowAdaptive}: Changes the velocity of the reflected particles in the species defined by the \textbf{incident} parameter to maintain the mean velocity equal to the average mean velocity of the others in the case.
+              The required input is:
+                
+                \textbf{incident}: Character.
+                  Name of the incident species that will be reflected.
+                  
             \end{itemize}
           \end{itemize}
         \item \textbf{EM}. Array of objects.determines the boundary conditions for the electromagnetic field.

src/modules/init/moduleInput.f90

@@ -864,7 +864,7 @@ MODULE moduleInput
               CALL config%get(object // '.specificHeat', cw, found)
               IF (.NOT. found) CALL criticalError("specificHeat not found for wallTemperature boundary type", 'readBoundary')
 
-              CALL initWallTemperature(bound, Tw, cw)
+              CALL wallTemperature_init(bound, Tw, cw)
 
             CASE('ionization')
               !Neutral parameters
@@ -882,10 +882,10 @@ MODULE moduleInput
               CALL config%get(object // '.electronSecondary', electronSecondary, found)
               IF (found) THEN
                 electronSecondaryID = speciesName2Index(electronSecondary)
-                CALL initIonization(bound, me/m_ref, m0, n0, v0, T0, &
+                CALL ionization_init(bound, me/m_ref, m0, n0, v0, T0, &
                                            speciesID, effTime, crossSection, eThreshold, electronSecondaryID)
               ELSE
-                CALL initIonization(bound, me/m_ref, m0, n0, v0, T0, &
+                CALL ionization_init(bound, me/m_ref, m0, n0, v0, T0, &
                                            speciesID, effTime, crossSection, eThreshold)
 
               END IF
@@ -898,10 +898,15 @@ MODULE moduleInput
 
               s_incident = speciesName2Index(speciesName)
 
-              call initQuasiNeutrality(bound, s_incident)
+              call quasiNeutrality_init(bound, s_incident)
 
             CASE('outflowAdaptive')
-              ALLOCATE(boundaryOutflowAdaptive:: bound)
+              call config%get(object // '.incident', speciesName, found)
+              if (.not. found) call criticalError("Incident species name not found for quasiNeutrality boundary model", 'readBoundary')
+
+              s_incident = speciesName2Index(speciesName)
+
+              call outflowAdaptive_init(bound, s_incident)
 
             CASE DEFAULT
               CALL criticalError('Boundary type ' // bType // ' undefined', 'readBoundary')
@@ -1295,6 +1300,7 @@ MODULE moduleInput
               ! If the boundary for the species is linked to the one analysing, add the edges
               if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then
                 bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
+
               end if
 
             end do
@@ -1304,6 +1310,24 @@ MODULE moduleInput
           allocate(bound%alpha(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
           bound%alpha = 0.d0
 
+        type is(boundaryOutflowAdaptive)
+          ! Loop over all physical surfaces
+          do ps = 1, nPhysicalSurfaces
+            ! Loop over all species
+            do s = 1, nSpecies
+              ! If the boundary for the species is linked to the one analysing, add the edges
+              if (associated(physicalSurfaces(ps)%particles(s)%obj, bound)) then
+                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
+
+              end if
+
+            end do
+
+          end do
+
+          allocate(bound%velocity_shift(1:mesh%numEdges)) ! TODO: Change this so only the edges associated to the boundary are here
+          bound%velocity_shift = 0.d0
+
         end select
 
       end do

src/modules/mesh/moduleMesh.f90

@@ -664,13 +664,15 @@ MODULE moduleMesh
   end type boundaryParticleGeneric
 
   interface
-    module subroutine initWallTemperature(boundary, T, c)
+    ! Init interfaces
+    ! wallTemeprature
+    module subroutine wallTemperature_init(boundary, T, c)
       CLASS(boundaryParticleGeneric), ALLOCATABLE, INTENT(inout):: boundary
       REAL(8), INTENT(in):: T, c !Wall temperature and specific heat
 
-    end subroutine initWallTemperature
+    end subroutine wallTemperature_init
 
-    module subroutine initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
+    module subroutine ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
       real(8), intent(in):: mImpact
       real(8), intent(in):: m0, n0, v0(1:3), T0 !Neutral properties
@@ -680,14 +682,23 @@ MODULE moduleMesh
       real(8), intent(in):: eThreshold
       integer, optional, intent(in):: electronSecondary
 
-    end subroutine initIonization
+    end subroutine ionization_init
 
-    module subroutine initQuasiNeutrality(boundary, s_incident)
+    ! quasiNeutrality
+    module subroutine quasiNeutrality_init(boundary, s_incident)
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
       integer, intent(in):: s_incident
 
-    end subroutine initQuasiNeutrality
+    end subroutine quasiNeutrality_init
 
+    ! outflowAdaptive
+    module subroutine outflowAdaptive_init(boundary, s_incident)
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      integer, intent(in):: s_incident
+
+    end subroutine outflowAdaptive_init
+
+    ! Get the index of the species from its name
     module function boundaryParticleName_to_Index(boundaryName) result(bp)
       character(:), allocatable:: boundaryName
       integer:: bp
@@ -697,6 +708,7 @@ MODULE moduleMesh
   end interface
 
   abstract interface
+    ! Apply the effects of the boundary to the particle
     subroutine applyParticle_interface(self, edge, part)
       use moduleSpecies
       import boundaryParticleGeneric, meshEdge
@@ -715,7 +727,7 @@ MODULE moduleMesh
 
     end subroutine updateParticle_interface
 
-    ! Update the values of the particle boundary model
+    ! Write the values of the particle boundary model
     subroutine printParticle_interface(self, fileID)
       import boundaryParticleGeneric
 
@@ -729,28 +741,28 @@ MODULE moduleMesh
   !Reflecting boundary
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryReflection
     CONTAINS
-      procedure, pass:: apply => reflection
+      procedure, pass:: apply => reflection_apply
 
   END TYPE boundaryReflection
 
   !Absorption boundary
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAbsorption
     CONTAINS
-      procedure, pass:: apply => absorption
+      procedure, pass:: apply => absorption_apply
 
   END TYPE boundaryAbsorption
 
   !Transparent boundary
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryTransparent
     CONTAINS
-      procedure, pass:: apply => transparent
+      procedure, pass:: apply => transparent_apply
 
   END TYPE boundaryTransparent
 
   !Symmetry axis
   TYPE, PUBLIC, EXTENDS(boundaryParticleGeneric):: boundaryAxis
     CONTAINS
-      procedure, pass:: apply => symmetryAxis
+      procedure, pass:: apply => symmetryAxis_apply
 
   END TYPE boundaryAxis
 
@@ -759,7 +771,7 @@ MODULE moduleMesh
     !Thermal velocity of the wall: square root(Wall temperature X specific heat)
     REAL(8):: vTh
     CONTAINS
-      procedure, pass:: apply => wallTemperature
+      procedure, pass:: apply => wallTemperature_apply
 
   END TYPE boundaryWallTemperature
 
@@ -775,7 +787,7 @@ MODULE moduleMesh
     integer:: tally
     integer(KIND=OMP_LOCK_KIND):: tally_lock
     CONTAINS
-      procedure, pass:: apply => ionization
+      procedure, pass:: apply => ionization_apply
 
   END TYPE boundaryIonization
 
@@ -785,91 +797,92 @@ MODULE moduleMesh
     integer:: s_incident ! species index of the incident species
     type(meshEdgePointer), allocatable:: edges(:) !Array with edges
     contains
-      procedure, pass:: apply => quasiNeutrality
+      procedure, pass:: apply => quasiNeutrality_apply
 
   end type boundaryQuasiNeutrality
 
   !Boundary for quasi-neutral outflow adjusting reflection coefficient
   type, public, extends(boundaryParticleGeneric):: boundaryOutflowAdaptive
-    real(8):: outflowCurrent
+    real(8), allocatable:: velocity_shift(:)
-    real(8):: reflectionFraction
+    integer:: s_incident ! species index of the incident species
+    type(meshEdgePointer), allocatable:: edges(:) !Array with edges
     contains
-      procedure, pass:: apply => outflowAdaptive
+      procedure, pass:: apply => outflowAdaptive_apply
 
   end type boundaryOutflowAdaptive
 
   interface
-    module subroutine reflection(self, edge, part)
+    module subroutine reflection_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryReflection), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine reflection
+    end subroutine reflection_apply
 
-    module subroutine absorption(self, edge, part)
+    module subroutine absorption_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryAbsorption), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine absorption
+    end subroutine absorption_apply
 
-    module subroutine transparent(self, edge, part)
+    module subroutine transparent_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryTransparent), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine transparent
+    end subroutine transparent_apply
 
-    module subroutine symmetryAxis(self, edge, part)
+    module subroutine symmetryAxis_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryAxis), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine symmetryAxis
+    end subroutine symmetryAxis_apply
 
-    module subroutine wallTemperature(self, edge, part)
+    module subroutine wallTemperature_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryWallTemperature), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine wallTemperature
+    end subroutine wallTemperature_apply
 
-    module subroutine ionization(self, edge, part)
+    module subroutine ionization_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryIonization), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine ionization
+    end subroutine ionization_apply
 
-    module subroutine quasiNeutrality(self, edge, part)
+    module subroutine quasiNeutrality_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryQuasiNeutrality), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine quasiNeutrality
+    end subroutine quasiNeutrality_apply
 
-    module subroutine outflowAdaptive(self, edge, part)
+    module subroutine outflowAdaptive_apply(self, edge, part)
       use moduleSpecies
 
       class(boundaryOutflowAdaptive), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
 
-    end subroutine outflowAdaptive
+    end subroutine outflowAdaptive_apply
 
     ! Generic basic boundary conditions to use internally in the code
     module subroutine genericReflection(edge, part) 

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -26,7 +26,7 @@ submodule(moduleMesh) boundaryParticle
     
     end function boundaryParticleName_to_Index
 
-    module SUBROUTINE reflection(self, edge, part)
+    module SUBROUTINE reflection_apply(self, edge, part)
       USE moduleCaseParam
       USE moduleSpecies
       IMPLICIT NONE
@@ -50,10 +50,10 @@ submodule(moduleMesh) boundaryParticle
       !particle is assumed to be inside
       part%n_in = .TRUE.
 
-    END SUBROUTINE reflection
+    END SUBROUTINE reflection_apply
 
     !Absoption in a surface
-    module SUBROUTINE absorption(self, edge, part)
+    module SUBROUTINE absorption_apply(self, edge, part)
       USE moduleCaseParam
       USE moduleSpecies
       IMPLICIT NONE
@@ -87,10 +87,10 @@ submodule(moduleMesh) boundaryParticle
 
       END IF
 
-    END SUBROUTINE absorption
+    END SUBROUTINE absorption_apply
 
     !Transparent boundary condition
-    module SUBROUTINE transparent(self, edge, part)
+    module SUBROUTINE transparent_apply(self, edge, part)
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -101,12 +101,12 @@ submodule(moduleMesh) boundaryParticle
       !Removes particle from domain
       part%n_in = .FALSE.
 
-    END SUBROUTINE transparent
+    END SUBROUTINE transparent_apply
 
     !Symmetry axis. Reflects particles.
     !Although this function should never be called, it is set as a reflective boundary
     !to properly deal with possible particles reaching a corner and selecting this boundary.
-    module SUBROUTINE symmetryAxis(self, edge, part)
+    module SUBROUTINE symmetryAxis_apply(self, edge, part)
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -116,12 +116,12 @@ submodule(moduleMesh) boundaryParticle
 
       CALL genericReflection(edge, part)
 
-    END SUBROUTINE symmetryAxis
+    END SUBROUTINE symmetryAxis_apply
 
     ! wallTemperature
     ! During the collision, the boundary transfers part of the energy to the incident particle
     ! Init
-    module SUBROUTINE initWallTemperature(boundary, T, c)
+    module SUBROUTINE wallTemperature_init(boundary, T, c)
       USE moduleRefParam
       IMPLICIT NONE
 
@@ -138,10 +138,10 @@ submodule(moduleMesh) boundaryParticle
 
       end select
 
-    END SUBROUTINE initWallTemperature
+    END SUBROUTINE wallTemperature_init
 
     ! Apply
-    module SUBROUTINE wallTemperature(self, edge, part)
+    module SUBROUTINE wallTemperature_apply(self, edge, part)
       USE moduleSpecies
       USE moduleRandom
       IMPLICIT NONE
@@ -159,12 +159,12 @@ submodule(moduleMesh) boundaryParticle
           
       CALL genericReflection(edge, part)
 
-    END SUBROUTINE wallTemperature
+    END SUBROUTINE wallTemperature_apply
 
     ! ionization
     ! An electron will pass through the surface and create a series of ion-electron pairs based on a neutral background
     !init
-    module SUBROUTINE initIonization(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
+    module SUBROUTINE ionization_init(boundary, mImpact, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
       use moduleRefParam
       use moduleSpecies
       use moduleCaseParam
@@ -212,7 +212,7 @@ submodule(moduleMesh) boundaryParticle
         boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
         boundary%deltaV = DSQRT(boundary%eThreshold/mImpact)
 
-        boundary%update => ionization_resetTally
+        boundary%update => ionization_update
         boundary%print  => ionization_print
 
         boundary%tally = 0
@@ -220,10 +220,10 @@ submodule(moduleMesh) boundaryParticle
 
       END SELECT
 
-    END SUBROUTINE initIonization
+    END SUBROUTINE ionization_init
 
     ! Apply
-    module SUBROUTINE ionization(self, edge, part)
+    module SUBROUTINE ionization_apply(self, edge, part)
       USE moduleList
       USE moduleSpecies
       USE moduleMesh
@@ -324,10 +324,10 @@ submodule(moduleMesh) boundaryParticle
       !Removes ionizing electron regardless the number of pair created
       part%n_in = .FALSE.
 
-    END SUBROUTINE ionization
+    END SUBROUTINE ionization_apply
 
     ! Update
-    subroutine ionization_resetTally(self)
+    subroutine ionization_update(self)
       implicit none
 
       class(boundaryParticleGeneric), intent(inout):: self
@@ -338,7 +338,7 @@ submodule(moduleMesh) boundaryParticle
 
       end select
 
-    end subroutine ionization_resetTally
+    end subroutine ionization_update
 
     ! Print
     subroutine ionization_print(self, fileID)
@@ -360,7 +360,7 @@ submodule(moduleMesh) boundaryParticle
     ! quasiNeutrality
     ! Maintains n_incident = n_rest by modifying the reflection parameter of the edge.
     ! Init
-    module subroutine initQuasiNeutrality(boundary, s_incident)
+    module subroutine quasiNeutrality_init(boundary, s_incident)
       implicit none
 
       class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
@@ -379,10 +379,10 @@ submodule(moduleMesh) boundaryParticle
 
       end select
 
-    end subroutine initQuasiNeutrality
+    end subroutine quasiNeutrality_init
 
     ! Apply
-    module subroutine quasiNeutrality(self, edge, part)
+    module subroutine quasiNeutrality_apply(self, edge, part)
       use moduleRandom
       implicit none
 
@@ -398,7 +398,7 @@ submodule(moduleMesh) boundaryParticle
 
       end if
 
-    end subroutine quasiNeutrality
+    end subroutine quasiNeutrality_apply
 
     ! Update
     subroutine quasiNeutrality_update(self)
@@ -416,7 +416,7 @@ submodule(moduleMesh) boundaryParticle
       select type(self)
       type is(boundaryQuasiNeutrality)
         do e = 1, size(self%edges)
-          edge => mesh%edges(e)%obj
+          edge => self%edges(e)%obj
 
           density_incident = 0.d0
           density_rest     = 0.d0
@@ -425,7 +425,8 @@ submodule(moduleMesh) boundaryParticle
           allocate(density_nodes(1:edge%nNodes))
           do s = 1, nSpecies
             do n = 1, edge%nNodes
-              node => mesh%nodes(n)%obj
+              node => mesh%nodes(nodes(n))%obj
+
               density_nodes(n) = node%output(s)%den
 
             end do
@@ -479,7 +480,7 @@ submodule(moduleMesh) boundaryParticle
       type is(boundaryQuasiNeutrality)
         write(fileID, '(A,",",A)') '"Edge id"', '"alpha"'
         do e = 1, size(self%edges)
-          write(fileID, '('//fmtColInt//','//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)
+          write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)
 
         end do
 
@@ -489,7 +490,30 @@ submodule(moduleMesh) boundaryParticle
 
     ! outflowAdaptive
     ! Adjust the reflection coefficient of the boundary to maintain a quasi-neutral outflow
-    module subroutine outflowAdaptive(self, edge, part)
+    ! Init
+    subroutine outflowAdaptive_init(boundary, s_incident)
+      implicit none
+
+      class(boundaryParticleGeneric), allocatable, intent(inout):: boundary
+      integer, intent(in):: s_incident
+
+      allocate(boundaryOutflowAdaptive:: boundary)
+      
+      select type(boundary)
+      type is(boundaryOutflowAdaptive)
+        allocate(boundary%edges(0))
+
+        boundary%s_incident = s_incident
+
+        boundary%update => outflowAdaptive_update
+        boundary%print  => outflowAdaptive_print
+
+      end select
+
+    end subroutine outflowAdaptive_init
+
+    ! Apply
+    module subroutine outflowAdaptive_apply(self, edge, part)
       use moduleSpecies
       use moduleRefParam, only: v_ref
       implicit none
@@ -497,25 +521,106 @@ submodule(moduleMesh) boundaryParticle
       class(boundaryOutflowAdaptive), intent(inout):: self
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
-      real(8):: v_cut
 
-      v_cut = 2.d0 * 40.d3/v_ref !This will be the drag velocity of the ions in the future
-
-      if (dot_product(part%v,-edge%normal) > v_cut) then
-        ! print *,part%v(1), v_cut
-        ! part%v = part%v + v_cut*edge%normal
-        part%v(1) = part%v(1) - v_cut
-        ! print *,part%v(1)
       call genericReflection(edge, part)
-        ! print *,part%v(1)
-        ! print *
 
-      else
+      part%v(1) = part%v(1) - self%velocity_shift(edge%n)
+
+      if (dot_product(part%v,edge%normal) <= 0.d0) then
+
         call genericTransparent(edge, part)
 
       end if
 
-    end subroutine outflowAdaptive
+    end subroutine outflowAdaptive_apply
+
+    ! Update
+    subroutine outflowAdaptive_update(self)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+      integer:: e, n, s
+      integer, allocatable:: nodes(:)
+      class(meshEdge), pointer:: edge
+      real(8), allocatable:: den_nodes(:)
+      real(8), allocatable:: mom_nodes(:)
+      class(meshNode), pointer:: node
+      real(8):: den_center, mom_center, vel_center
+      real(8):: velocity_incident, velocity_rest
+      real(8):: alpha
+      
+      select type(self)
+      type is(boundaryOutflowAdaptive)
+        do e = 1, size(self%edges)
+          edge => self%edges(e)%obj
+
+          velocity_incident = 0.d0
+          velocity_rest     = 0.d0
+
+          nodes = edge%getNodes(edge%nNodes)
+          allocate(den_nodes(1:edge%nNodes))
+          allocate(mom_nodes(1:edge%nNodes))
+          do s = 1, nSpecies
+            do n = 1, edge%nNodes
+              node => mesh%nodes(nodes(n))%obj
+
+              den_nodes(n) = node%output(s)%den
+              mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
+
+            end do
+
+            den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
+            mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
+
+            if (den_center > 1.0d-10) then
+              vel_center = mom_center / den_center
+
+            else
+              vel_center = 0.0d0
+            end if
+
+            if (s == self%s_incident) then
+              velocity_incident = vel_center
+            else
+              velocity_rest = velocity_rest + vel_center
+            end if
+
+          end do
+
+          alpha = velocity_rest - velocity_incident
+
+          ! Apply correction with a factor of 0.1 to avoid fast changes
+          self%velocity_shift(edge%n) = self%velocity_shift(edge%n) + 1.0d-2 * alpha
+
+          deallocate(den_nodes)
+          deallocate(mom_nodes)
+
+        end do
+
+      end select
+
+    end subroutine outflowAdaptive_update
+
+    ! Print
+    subroutine outflowAdaptive_print(self, fileID)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+      integer, intent(in):: fileID
+      integer:: e
+
+      write(fileID, '(A)') self%name
+      select type(self)
+      type is(boundaryOutflowAdaptive)
+        write(fileID, '(A,",",A)') '"Edge id"', '"vel shift non-dimensional"'
+        do e = 1, size(self%edges)
+          write(fileID, '('//fmtColInt//','//fmtColReal//')') self%edges(e)%obj%n, self%velocity_shift(self%edges(e)%obj%n)
+        end do
+
+      end select
+
+    end subroutine outflowAdaptive_print
+
 
     ! Generic boundary conditions for internal use
     ! reflection