Merge branch 'development' into feature/0DGrid

Conflicts:
	src/modules/mesh/moduleMesh.f90

doc/user-manual/fpakc_UserManual.tex

@@ -618,25 +618,12 @@ make
                                   Required values are:
                                   \begin{itemize}
                                     \item \textbf{speciesName}: Character.
-                                                                Name of species.
+                                                                Name of species as defined in the object \textbf{species}.
                                     \item \textbf{initialState}: Character.
-                                                                 Plain text file that contains the information about the species macroscopic properties in the grid.
-                                                                 Initial particles are assumed to be Maxwellian.
+                                                                 Output file from previous run used as an initial state for the species.
+                                                                 The file format must be the same as in \textbf{geometry.meshType}
+                                                                 Initial particles are assumed to have a Maxwellian distribution.
                                                                  File must be located at \textbf{output.path}.
-                                                                 The file must contain the following columns in this specific order:
-                                                                 \begin{itemize}
-                                                                   \item \textit{Element}: Integer. 
-                                                                                           Identifier of the volume in the mesh.
-                                                                                           It is not required to specify empty volumes.
-                                                                   \item \textit{Density}: Real.
-                                                                                           Species density in $\unit{m^-3}$.
-                                                                   \item \textit{Velocity}: Real.
-                                                                                            Three colums representing the initial velocity in each direction.
-                                                                                            Units are $\unit{m s^-1}$.
-                                                                   \item \textit{Temperature}: Real.
-                                                                                               One column that represents the initial temperature in $\unit{K}$.
-                                                                                 
-                                                                 \end{itemize}
                                   \end{itemize}
         \end{itemize}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

src/modules/mesh/0D/moduleMesh0D.f90

@@ -17,6 +17,7 @@ MODULE moduleMesh0D
       PROCEDURE, PASS:: init => initVol0D
       PROCEDURE, PASS:: getNodes    => getNodes0D
       PROCEDURE, PASS:: randPos     => randPos0D
+      PROCEDURE, NOPASS:: fPsi      => fPsi0D
       PROCEDURE, PASS:: scatter     => scatter0D
       PROCEDURE, PASS:: gatherEF    => gatherEF0D
       PROCEDURE, PASS:: elemK       => elemK0D
@@ -100,6 +101,15 @@ MODULE moduleMesh0D
 
     END FUNCTION randPos0D
 
+    PURE FUNCTION fPsi0D(xi) RESULT(fPsi)
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+      ALLOCATE(fPsi(1:1))
+      fPsi = 1.D0
+
+    END FUNCTION fPsi0D
+
     SUBROUTINE scatter0D(self, part)
       USE moduleMath
       USE moduleSpecies

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -36,19 +36,12 @@ MODULE moduleMesh1DCart
     CONTAINS
       PROCEDURE, PASS:: detJac => detJ1DCart
       PROCEDURE, PASS:: invJac => invJ1DCart
-      PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
       PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
       PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
   
   END TYPE meshVol1DCart
 
   ABSTRACT INTERFACE
-    PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
-      REAL(8), INTENT(in):: xi(1:3)
-      REAL(8), ALLOCATABLE:: fPsi(:)
-
-    END FUNCTION fPsi_interface
-
     PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
       REAL(8), INTENT(in):: xi(1:3)
       REAL(8), ALLOCATABLE:: dPsi(:,:)

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -36,19 +36,12 @@ MODULE moduleMesh1DRad
     CONTAINS
       PROCEDURE, PASS:: detJac => detJ1DRad
       PROCEDURE, PASS:: invJac => invJ1DRad
-      PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
       PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
       PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
   
   END TYPE meshVol1DRad
 
   ABSTRACT INTERFACE
-    PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
-      REAL(8), INTENT(in):: xi(1:3)
-      REAL(8), ALLOCATABLE:: fPsi(:)
-
-    END FUNCTION fPsi_interface
-
     PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
       REAL(8), INTENT(in):: xi(1:3)
       REAL(8), ALLOCATABLE:: dPsi(:,:)

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -41,19 +41,12 @@ MODULE moduleMesh2DCart
     CONTAINS
       PROCEDURE, PASS:: detJac    => detJ2DCart
       PROCEDURE, PASS:: invJac    => invJ2DCart
-      PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
       PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
       PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
 
   END TYPE meshVol2DCart
 
   ABSTRACT INTERFACE
-    PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
-      REAL(8), INTENT(in):: xi(1:3)
-      REAL(8), ALLOCATABLE:: fPsi(:)
-
-    END FUNCTION fPsi_interface
-
     PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
       REAL(8), INTENT(in):: xi(1:3)
       REAL(8), ALLOCATABLE:: dPsi(:,:)

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -41,19 +41,12 @@ MODULE moduleMesh2DCyl
     CONTAINS
       PROCEDURE, PASS:: detJac    => detJ2DCyl
       PROCEDURE, PASS:: invJac    => invJ2DCyl
-      PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
       PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
       PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
 
   END TYPE meshVol2DCyl
 
   ABSTRACT INTERFACE
-    PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
-      REAL(8), INTENT(in):: xi(1:3)
-      REAL(8), ALLOCATABLE:: fPsi(:)
-
-    END FUNCTION fPsi_interface
-
     PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
       REAL(8), INTENT(in):: xi(1:3)
       REAL(8), ALLOCATABLE:: dPsi(:,:)

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -35,19 +35,12 @@ MODULE moduleMesh3DCart
     CONTAINS
       PROCEDURE, PASS:: detJac => detJ3DCart
       PROCEDURE, PASS:: invJac => invJ3DCart
-      PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
       PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
       PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
   
   END TYPE meshVol3DCart
 
   ABSTRACT INTERFACE
-    PURE FUNCTION fPsi_interface(xii) RESULT(fPsi)
-      REAL(8), INTENT(in):: xii(1:3)
-      REAL(8), ALLOCATABLE:: fPsi(:)
-
-    END FUNCTION fPsi_interface
-
     PURE FUNCTION dPsi_interface(xii) RESULT(dPsi)
       REAL(8), INTENT(in):: xii(1:3)
       REAL(8), ALLOCATABLE:: dPsi(:,:)
@@ -327,18 +320,18 @@ MODULE moduleMesh3DCart
     END SUBROUTINE volumeTetra
 
     !Computes element functions in point xii
-    PURE FUNCTION fPsiTetra(xii) RESULT(fPsi)
+    PURE FUNCTION fPsiTetra(xi) RESULT(fPsi)
       IMPLICIT NONE
 
-      REAL(8), INTENT(in):: xii(1:3)
+      REAL(8), INTENT(in):: xi(1:3)
       REAL(8), ALLOCATABLE:: fPsi(:)
 
       ALLOCATE(fPsi(1:4))
 
-      fPsi(1) = 1.D0 - xii(1) - xii(2) - xii(3)
-      fPsi(2) =        xii(1)
-      fPsi(3) =                 xii(2)
-      fPsi(4) =                          xii(3)
+      fPsi(1) = 1.D0 - xi(1) - xi(2) - xi(3)
+      fPsi(2) =        xi(1)
+      fPsi(3) =                xi(2)
+      fPsi(4) =                        xi(3)
 
     END FUNCTION fPsiTetra
 

src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90

@@ -15,6 +15,7 @@ MODULE moduleMeshInputGmsh2
       TYPE IS(meshParticles)
         self%printOutput => printOutputGmsh2
         self%printEM     => printEMGmsh2
+        self%readInitial => readInitialGmsh2
 
       END SELECT
       self%readMesh    => readGmsh2
@@ -289,4 +290,76 @@ MODULE moduleMeshInputGmsh2
 
     END SUBROUTINE readGmsh2
 
+    !Reads the initial information from an output file for an species
+    SUBROUTINE readInitialGmsh2(sp, filename, density, velocity, temperature)
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      INTEGER, INTENT(in):: sp
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
+      INTEGER:: i, e
+      INTEGER:: numNodes
+
+      OPEN(10, file = TRIM(filename))
+
+      !Skip first lines
+      DO i = 1, 11
+        READ(10, *)
+
+      END DO
+
+      !Reads number of nodes in file
+      READ(10, *) numNodes
+      ALLOCATE(density(1:numNodes))
+      ALLOCATE(velocity(1:numNodes, 1:3))
+      ALLOCATE(temperature(1:numNodes))
+
+      DO i = 1, numNodes
+        !Reads the density
+        READ(10, *), e, density(i)
+
+      END DO
+
+      DO i = 1, 10
+        READ(10, *)
+
+      END DO
+
+      DO i = 1, numNodes
+        !Reads the velocity
+        READ(10, *), e, velocity(i, 1:3)
+
+      END DO
+
+      !Skip uneccessary lines
+      DO i = 1, 10
+        READ(10, *)
+
+      END DO
+
+      !Assign density to nodes
+      DO i = 1, numNodes
+        !Skips pressure
+        READ(10, *)
+
+      END DO
+
+      !Skip uneccessary lines
+      DO i = 1, 10
+        READ(10, *)
+
+      END DO
+
+      !Assign density to nodes
+      DO i = 1, numNodes
+        !Skips pressure
+        READ(10, *) e, temperature(i)
+
+      END DO
+
+    END SUBROUTINE readInitialGmsh2
+
 END MODULE moduleMeshInputGmsh2

src/modules/mesh/moduleMesh.f90

@@ -161,6 +161,7 @@ MODULE moduleMesh
       PROCEDURE(initVol_interface),     DEFERRED, PASS::   init
       PROCEDURE(getNodesVol_interface), DEFERRED, PASS::   getNodes
       PROCEDURE(randPosVol_interface),  DEFERRED, PASS::   randPos
+      PROCEDURE(fPsi_interface),        DEFERRED, NOPASS:: fPsi
       PROCEDURE(scatter_interface),     DEFERRED, PASS::   scatter
       PROCEDURE(gatherEF_interface),    DEFERRED, PASS::   gatherEF
       PROCEDURE(elemK_interface),       DEFERRED, PASS::   elemK
@@ -207,6 +208,12 @@ MODULE moduleMesh
 
     END FUNCTION getNodesVol_interface
 
+    PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
+      REAL(8), INTENT(in):: xi(1:3)
+      REAL(8), ALLOCATABLE:: fPsi(:)
+
+    END FUNCTION fPsi_interface
+
     PURE FUNCTION elemK_interface(self) RESULT(localK)
       IMPORT:: meshVol
       CLASS(meshVol), INTENT(in):: self
@@ -270,9 +277,10 @@ MODULE moduleMesh
     TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
     !Array of volume elements
     TYPE(meshVolCont),  ALLOCATABLE:: vols(:)
-    PROCEDURE(readMesh_interface),     POINTER, PASS:: readMesh     => NULL()
-    PROCEDURE(connectMesh_interface),  POINTER, PASS:: connectMesh  => NULL()
-    PROCEDURE(printColl_interface),    POINTER, PASS:: printColl    => NULL()
+    PROCEDURE(readMesh_interface),     POINTER, PASS::   readMesh     => NULL()
+    PROCEDURE(readInitial_interface),  POINTER, NOPASS:: readInitial  => NULL()
+    PROCEDURE(connectMesh_interface),  POINTER, PASS::   connectMesh  => NULL()
+    PROCEDURE(printColl_interface),    POINTER, PASS::   printColl    => NULL()
     CONTAINS
       PROCEDURE, PASS:: doCollisions
 
@@ -288,6 +296,15 @@ MODULE moduleMesh
 
     END SUBROUTINE readMesh_interface
 
+    SUBROUTINE readInitial_interface(sp, filename, density, velocity, temperature)
+      INTEGER, INTENT(in):: sp
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
+      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
+
+    END SUBROUTINE readInitial_interface
+
     !Connects volume and edges to the mesh
     SUBROUTINE connectMesh_interface(self)
       IMPORT meshGeneric

src/modules/moduleInput.f90

@@ -206,7 +206,6 @@ MODULE moduleInput
       CALL config%get(object // '.WeightingScheme', WSType, found)
       CALL solver%initWS(WSType)
 
-
       !Makes tau(s) non-dimensional
       tau = tau / ti_ref
       tauMin = tauMin / ti_ref
@@ -219,7 +218,6 @@ MODULE moduleInput
       !Read initial state for species
       CALL verboseError('Reading Initial state...')
       CALL readInitial(config)
-      CALL checkStatus(config, "readInitial")
 
     END SUBROUTINE readCase
 
@@ -237,14 +235,21 @@ MODULE moduleInput
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: object
       INTEGER:: nInitial
-      INTEGER:: i, p, e
+      INTEGER:: i, j, p, e
       CHARACTER(LEN=2):: iString
       CHARACTER(:), ALLOCATABLE:: spName
       INTEGER:: sp
       CHARACTER(:), ALLOCATABLE:: spFile
-      INTEGER:: stat
-      CHARACTER(100):: dummy
-      REAL(8):: density, velocity(1:3), temperature
+      CHARACTER(:), ALLOCATABLE:: filename
+      REAL(8), ALLOCATABLE, DIMENSION(:):: density, temperature
+      REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity
+      INTEGER, ALLOCATABLE, DIMENSION(:):: nodes
+      INTEGER:: nNodes
+      REAL(8), ALLOCATABLE, DIMENSION(:):: source, fPsi
+      !Density at the volume centroid
+      REAL(8):: densityCen
+      !Mean velocity and temperature at particle position
+      REAL(8):: velocityXi(1:3), temperatureXi
       INTEGER:: nNewPart = 0.D0
       TYPE(particle), POINTER:: partNew
       REAL(8):: vTh
@@ -260,38 +265,76 @@ MODULE moduleInput
           CALL config%get(object // '.speciesName', spName, found)
           sp = speciesName2Index(spName)
           CALL config%get(object // '.initialState', spFile, found)
-          OPEN (10, FILE = path // spFile, ACTION = 'READ')
-          DO
-            READ(10, '(A)', IOSTAT = stat) dummy
-            !If EoF, exit reading
-            IF (stat /= 0) EXIT
-            !If comment, skip
-            IF (INDEX(dummy,'#') /= 0) CYCLE
-            !Go up one line
-            BACKSPACE(10)
-            !Read information
-            READ(10, *) e, density, velocity, temperature
+          !Reads node values at the nodes
+          filename = path // spFile
+          CALL mesh%readInitial(sp, filename, density, velocity, temperature)
+          !For each volume in the node, create corresponding particles
+          DO e = 1, mesh%numVols
             !Scale variables
-            !Particles in cell volume
-            nNewPart = INT(density * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
-            !Non-dimensional velocity
-            velocity = velocity / v_ref
-            !Non-dimensional temperature
-            temperature = temperature / T_ref
-            !Non-dimensional thermal temperature
-            vTh = DSQRT(temperature/species(sp)%obj%m)
+            !Density at centroid of cell
+            nodes = mesh%vols(e)%obj%getNodes()
+            nNodes = SIZE(nodes)
+            !TODO: Procedure to obtain centroid from element (also for printing Electric Field)
+            fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/))
+            ALLOCATE(source(1:nNodes))
+            DO j = 1, nNodes
+              source(j) = density(nodes(j))
+
+            END DO
+            densityCen = DOT_PRODUCT(fPsi, source)
+            DEALLOCATE(fPsi)
+
+            !Calculate number of particles
+            nNewPart = INT(densityCen * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight)
+
             !Allocate new particles
             DO p = 1, nNewPart
               ALLOCATE(partNew)
               partNew%species => species(sp)%obj
-              partNew%v(1)    =  velocity(1) + vTh*randomMaxwellian()
-              partNew%v(2)    =  velocity(2) + vTh*randomMaxwellian()
-              partNew%v(3)    =  velocity(3) + vTh*randomMaxwellian()
-              partNew%vol     =  e
               partNew%r       =  mesh%vols(e)%obj%randPos()
               partNew%xi      =  mesh%vols(e)%obj%phy2log(partNew%r)
-              partNew%n_in    =  .TRUE.
-              partNew%weight  =  species(sp)%obj%weight
+              !Get mean velocity at particle position
+              fPsi            =  mesh%vols(e)%obj%fPsi(partNew%xi)
+              DO j = 1, nNodes
+                source(j) = velocity(nodes(j), 1)
+
+              END DO
+              velocityXi(1) = DOT_PRODUCT(fPsi, source)
+              DO j = 1, nNodes
+                source(j) = velocity(nodes(j), 2)
+
+              END DO
+              velocityXi(2) = DOT_PRODUCT(fPsi, source)
+              DO j = 1, nNodes
+                source(j) = velocity(nodes(j), 3)
+
+              END DO
+              velocityXi(3) = DOT_PRODUCT(fPsi, source)
+
+              velocityXi = velocityXi / v_ref
+
+              !Get temperature at particle position
+              DO j = 1, nNodes
+                source(j) = temperature(nodes(j))
+
+              END DO
+              temperatureXi = DOT_PRODUCT(fPsi, source)
+              temperatureXi = temperatureXi / T_ref
+
+              vTh = DSQRT(temperatureXi / species(sp)%obj%m)
+              partNew%v(1)    =  velocityXi(1) + vTh*randomMaxwellian()
+              partNew%v(2)    =  velocityXi(2) + vTh*randomMaxwellian()
+              partNew%v(3)    =  velocityXi(3) + vTh*randomMaxwellian()
+              partNew%vol     =  e
+              IF (ASSOCIATED(meshForMCC, mesh)) THEN
+                partNew%volColl = partNew%vol
+
+              ELSE
+                partNew%volColl = findCellBrute(meshColl, partNew%r)
+
+              END IF
+              partNew%n_in   =  .TRUE.
+              partNew%weight =  species(sp)%obj%weight
               !If charged species, add qm to particle
               SELECT TYPE(sp => species(sp)%obj)
               TYPE IS (speciesCharged)
@@ -305,8 +348,11 @@ MODULE moduleInput
               !Assign particle to temporal list of particles
               CALL partInitial%add(partNew)
 
+
             END DO
 
+            DEALLOCATE(source)
+
           END DO
 
         END DO