diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 5f8e76a..8e8309b 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index a485c95..3249d02 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/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} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index bc86719..e945f25 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/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 diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index a6e4248..73a4c1d 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/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(:,:) diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 2bd9d5b..d8e5c51 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/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(:,:) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 89704d6..c830f8d 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/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(:,:) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 2784414..80c0665 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/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(:,:) diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index ad5c65d..e65ff25 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/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 diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 629a6f2..b446d76 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/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 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 53fa341..304dd86 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/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 diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 5272206..fed9228 100644 --- a/src/modules/moduleInput.f90 +++ b/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