Collisions and EM field in .vtu

The collisions and EM field information is now available in .vtu files.
A collection file .pvd is provided per dataset for time-dependent
plotting.

Still to do:
  Write average quantities in .vtu
  Read mesh from .vtu
This commit is contained in:
Jorge Gonzalez 2023-02-04 12:31:33 +01:00
commit f1c0c5755f
2 changed files with 233 additions and 61 deletions

View file

@ -1,22 +1,44 @@
MODULE moduleMeshOutputVTK
CHARACTER(LEN=6):: prefix = 'OUTPUT'
CONTAINS
SUBROUTINE writeFileHeader(self, fileID)
PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName)
USE moduleOutput
IMPLICIT NONE
CHARACTER(*), INTENT(in):: prefix, suffix, extension
INTEGER, INTENT(in), OPTIONAL:: t
CHARACTER (LEN=iterationDigits):: tString
CHARACTER(:), ALLOCATABLE:: fileName
IF (PRESENT(t)) THEN
WRITE(tString, iterationFormat) t
fileName = prefix // '_' // tString // '_' // suffix // '.' // extension
ELSE
fileName = prefix // '_' // suffix // '.' // extension
END IF
END FUNCTION formatFileName
SUBROUTINE writeHeader(nNodes, nCells, fileID)
USE moduleMesh
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: nNodes, nCells
INTEGER, INTENT(in):: fileID
WRITE(fileID,"(A)") '<?xml version="1.0"?>'
WRITE(fileID,"(2X, A)") '<VTKFile type="UnstructuredGrid">'
WRITE(fileID,"(4X, A,ES20.6E3,A)") '<UnstructuredGrid>'
WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', self%numNodes, '" NumberOfCells="', self%numCells, '">'
WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', nNodes, '" NumberOfCells="', nCells, '">'
END SUBROUTINE writeFileHeader
END SUBROUTINE writeHeader
SUBROUTINE writeFileFooter(fileID)
SUBROUTINE writeFooter(fileID)
IMPLICIT NONE
INTEGER, INTENT(in):: fileID
@ -25,7 +47,7 @@ MODULE moduleMeshOutputVTK
WRITE(fileID,"(4X, A)") '</UnstructuredGrid>'
WRITE(fileID,"(2X, A)") '</VTKFile>'
END SUBROUTINE writeFileFooter
END SUBROUTINE writeFooter
FUNCTION getCellType(cell) RESULT(indexType)
USE moduleMesh3DCart
@ -74,12 +96,12 @@ MODULE moduleMeshOutputVTK
END FUNCTION getCellType
SUBROUTINE writeFileMesh(self, fileID)
SUBROUTINE writeMesh(self, fileID)
USE moduleMesh
USE moduleRefParam
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: fileID
CHARACTER(LEN=25):: nodeFormat
CHARACTER(LEN=25):: cellFormat
@ -123,104 +145,252 @@ MODULE moduleMeshOutputVTK
WRITE(fileID, "(10X, A)") '</DataArray>'
WRITE(fileID, "(8X, A)") '</Cells>'
END SUBROUTINE writeFileMesh
END SUBROUTINE writeMesh
SUBROUTINE writeNodeData(self, fileID, output)
SUBROUTINE writeSpeciesOutput(self, fileID, speciesIndex)
USE moduleMesh
USE moduleOutput
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: fileID
INTEGER, INTENT(in):: speciesIndex
TYPE(outputFormat):: output(1:self%numNodes)
CHARACTER(LEN=25):: nodeFormat
INTEGER:: n
WRITE(fileID,"(A)") '<PointData>'
WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
DO n = 1, self%numNodes
CALL calculateOutput(self%nodes(n)%obj%output(speciesIndex), output(n), self%nodes(n)%obj%v, species(speciesIndex)%obj)
END DO
WRITE(fileID,"(8X,A)") '<PointData>'
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
WRITE(fileID,nodeFormat) (output(n)%density, n = 1, self%numNodes)
WRITE(fileID,"(A)") '</DataArray>'
WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Velocity (m s^-1)" NumberOfComponents="3">'
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Velocity (m s^-1)" NumberOfComponents="3">'
WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
WRITE(fileID,nodeFormat) (output(n)%velocity(1:3), n = 1, self%numNodes)
WRITE(fileID,"(A)") '</DataArray>'
WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Pressure (Pa)" NumberOfComponents="1">'
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Pressure (Pa)" NumberOfComponents="1">'
WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
WRITE(fileID,nodeFormat) (output(n)%pressure, n = 1, self%numNodes)
WRITE(fileID,"(A)") '</DataArray>'
WRITE(fileID,"(A)") '<DataArray type="Float64" Name="Temperature (K)" NumberOfComponents="1">'
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Temperature (K)" NumberOfComponents="1">'
WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
WRITE(fileID,nodeFormat) (output(n)%temperature, n = 1, self%numNodes)
WRITE(fileID,"(A)") '</DataArray>'
WRITE(fileID,"(A)") '</PointData>'
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(8X,A)") '</PointData>'
END SUBROUTINE writeNodeData
END SUBROUTINE writeSpeciesOutput
SUBROUTINE writeCollOutput(self,fileID)
USE moduleMesh
USE moduleCollisions
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: fileID
INTEGER:: k, c, n
CHARACTER(:), ALLOCATABLE:: title
CHARACTER (LEN=2):: cString
CHARACTER(LEN=25):: cellFormat
WRITE(fileID,"(8X,A)") '<CellData>'
DO k = 1, nCollPairs
DO c = 1, interactionMatrix(k)%amount
WRITE(cString, "(I2)") c
title = 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="',title, '" NumberOfComponents="1">'
WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(I10))"
WRITE(fileID, cellFormat) (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells)
WRITE(fileID, "(10X, A)") '</DataArray>'
END DO
END DO
WRITE(fileID,"(8X,A)") '</CellData>'
END SUBROUTINE writeCollOutput
SUBROUTINE writeEM(self, fileID)
USE moduleMesh
USE moduleRefParam
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: fileID
CHARACTER(LEN=25):: nodeFormat
CHARACTER(LEN=25):: cellFormat
INTEGER:: n
REAL(8):: Xi(1:3)
Xi = (/ 0.D0, 0.D0, 0.D0 /)
!Points in nodes
WRITE(fileID,"(8X,A)") '<PointData>'
!Electric potential
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Potential (V)" NumberOfComponents="1">'
WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
!Magnetic Field
WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Magnetic Field (T)" NumberOfComponents="3">'
WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes)
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(8X,A)") '</PointData>'
!Cell Data
WRITE(fileID,"(8X,A)") '<CellData>'
!Electric field
WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(3(ES20.6E3)))"
WRITE(fileID, cellFormat) (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
WRITE(fileID,"(10X,A)") '</DataArray>'
WRITE(fileID,"(8X,A)") '</CellData>'
END SUBROUTINE writeEM
SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
USE moduleCaseParam
USE moduleOutput
USE moduleRefParam
IMPLICIT NONE
INTEGER:: fileID
INTEGER, INTENT(in):: t
CHARACTER(*):: fileNameStep, fileNameCollection
IF (t == tInitial) THEN
!Create collection file
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection)
WRITE (fileID + 1, "(A)") '<VTKFile type="Collection">'
WRITE (fileID + 1, "(2X, A)") '<Collection>'
CLOSE(fileID + 1)
END IF
!Write iteration file in collection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND')
WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'
!Close collection file
IF (t == tFinal) THEN
WRITE (fileID + 1, "(2X, A)") '</Collection>'
WRITE (fileID + 1, "(A)") '</VTKFile>'
END IF
CLOSE(fileID + 1)
END SUBROUTINE writeCollection
SUBROUTINE printOutputVTK(self,t)
USE moduleMesh
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
USE moduleCaseParam
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i, fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
CHARACTER (LEN=iterationDigits):: tstring
TYPE(outputFormat):: output(1:self%numNodes)
fileID = 60
DO i = 1, nSpecies
WRITE(tstring, iterationFormat) t
fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.vtu'
fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
fileNameCollection = 'OUTPUT_Collection_' // species(i)%obj%name // '.pvd'
IF (t == tInitial) THEN
!Create collection file
WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection)
WRITE (fileID + 1, "(A)") '<VTKFile type="Collection">'
WRITE (fileID + 1, "(2X, A)") '<Collection>'
CLOSE(fileID + 1)
CALL writeHeader(self%numNodes, self%numCells, fileID)
END IF
CALL writeMesh(self, fileID)
!Write iteration file in collection
OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND')
WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)"), '<DataSet timestep="', t*ti_ref,'" file="', fileName,'"/>'
CALL writeSpeciesOutput(self, fileID, i)
!Close collection file
IF (t == tFinal) THEN
WRITE (fileID + 1, "(2X, A)") '</Collection>'
WRITE (fileID + 1, "(A)") '</VTKFile>'
END IF
CLOSE(fileID + 1)
CALL writeFileHeader(self, fileID)
CALL writeFileMesh(self, fileID)
DO n = 1, self%numNodes
CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
END DO
CALL writeNodeData(self, fileID, output)
CALL writeFileFooter(fileID)
CALL writeFooter(fileID)
CLOSE(fileID)
!Write collection file for time plotting
fileNameCollection = formatFileName(prefix, 'Collection_' // species(i)%obj%name, 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
END DO
END SUBROUTINE printOutputVTK
SUBROUTINE printCollVTK(self,t)
USE moduleMesh
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i, fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
CHARACTER (LEN=iterationDigits):: tstring
TYPE(outputFormat):: output(1:self%numNodes)
fileID = 62
IF (collOutput) THEN
fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
CALL writeHeader(self%numNodes, self%numCells, fileID)
CALL writeMesh(self, fileID)
CALL writeCollOutput(self, fileID)
CALL writeFooter(fileID)
CLOSE(fileID)
!Write collection file for time plotting
fileNameCollection = formatFileName(prefix, 'Collection_Collisions', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
END IF
END SUBROUTINE printCollVTK
SUBROUTINE printEMVTK(self, t)
USE moduleMesh
IMPLICIT NONE
CLASS(meshParticles), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: fileID
CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
fileID = 64
IF (emOutput) THEN
fileName = formatFileName(prefix, 'EMField', 'vtu', t)
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
OPEN (fileID, file = path // folder // '/' // fileName)
CALL writeHeader(self%numNodes, self%numCells, fileID)
CALL writeMesh(self, fileID)
CALL writeEM(self, fileID)
CALL writeFooter(fileID)
CLOSE(fileID)
!Write collection file for time plotting
fileNameCollection = formatFileName(prefix, 'Collection_EMField', 'pvd')
CALL writeCollection(fileID, t, fileName, filenameCollection)
END IF
END SUBROUTINE printEMVTK
END MODULE moduleMeshOutputVTK

View file

@ -528,9 +528,11 @@ MODULE moduleSolver
CALL outputProbes(t)
CALL mesh%printOutput(t)
CALL printOutputVTK(mesh, t)
CALL printOutputVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT
IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
IF (ASSOCIATED(meshForMCC)) CALL printCollVTK(meshForMCC,t) !TEMPORARY TO TEST VTK OUTPUT
CALL mesh%printEM(t)
CALL printEMVTK(mesh, t) !TEMPORARY TO TEST VTK OUTPUT
WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld
IF (t == 0) THEN