diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 49536bf..da9471c 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 8897869..a55c7dd 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -38,6 +38,7 @@ \newacronym{io}{I/O}{input/output} \newacronym{mcc}{MCC}{Monte-Carlo Collisions} \newacronym{cs}{CS}{Coulomb Scattering} +\newacronym{vtu}{VTU}{Visualization Toolkit for unstructured grids} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newglossaryentry{openmp}{name={OpenMP},description={Shared-memory parallelization}} @@ -437,25 +438,27 @@ make \begin{itemize} \item \textbf{Cart}: Cartesian coordinates. Available for \textbf{geometry.dimension} $1$, $2$ and $3$. - For \Gls{gmsh} mesh format, the coordinates $x$, $y$ and $z$ correspond to $x$, $y$ and $z$ respectively. + The coordinates $x$, $y$ and $z$ correspond to $x$, $y$ and $z$ respectively. \item \textbf{Cyl}: Cylindrical coordinates ($z \hyphen r$) with symmetry axis at $r = 0$. Only available for \textbf{geometry.dimension} $2$. - For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $z$ and $r$ respectively. + The coordinates $x$ and $y$ correspond to $z$ and $r$ respectively. \item \textbf{Rad}: One-dimensional radial space ($r$). Only available for \textbf{geometry.dimension} $1$. - For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$. + The coordinates $x$ corresponds to $r$. \end{itemize} \item \textbf{meshType}: Character. Format of mesh file. + The output will be written in the same format as specified here. Accepted formats are: \begin{itemize} - \item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0. + \item \textbf{gmsh2}: \Gls{gmsh} file format in version 2.0. This has to be in ASCII format. + \item \textbf{vtu}: \Gls{vtu} file format. This has to be in ASCII format. \end{itemize} \item \textbf{meshFile}: Character. Mesh filename. This file is searched in the path \textbf{output.path} and must contain the file extension. \item \textbf{volume}: Real. - Units of $\unit{m^-3}$. + Units of $\unit{m^{-3}}$. Used to set a fictitious volume for the $0$ dimension. Ignored in the other cases. \end{itemize} @@ -536,7 +539,7 @@ make Density of neutral background. Required parameter. \item \textbf{velocity}: Real. - Units in $\unit{m^{-3}}$. + Units in $\unit{m s^{-1}}$. Array of dimension $3$. Mean velocity of neutral background. Required parameter. @@ -604,7 +607,8 @@ make \item \textbf{part/s}: Particles (real) per second. \end{itemize} \item \textbf{v}: Real. - Module of velocity vector, in $\unitfrac{m}{s}$. + Units of $\unit{m s^{-1}}$. + Module of velocity vector. \item \textbf{n}: Real. Array dimension $3$. Direction of injection. @@ -618,6 +622,7 @@ make \item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction. \end{itemize} \item \textbf{T}: Real. + Units of $\unit{K}$ Array dimension $3$. Temperature in each direction. \item \textbf{physicalSurface}: Integer. diff --git a/src/makefile b/src/makefile index 2a421f8..39a39a7 100644 --- a/src/makefile +++ b/src/makefile @@ -5,6 +5,8 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \ $(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \ + $(OBJDIR)/moduleMeshInoutCommon.o \ + $(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \ $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ $(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \ $(OBJDIR)/moduleMesh3DCart.o \ diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 0a33903..da43e98 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -860,7 +860,13 @@ MODULE moduleInput SUBROUTINE readGeometry(config) USE moduleMesh USE moduleMeshInputGmsh2, ONLY: initGmsh2 + USE moduleMeshInputVTU, ONLY: initVTU USE moduleMeshInput0D, ONLY: init0D + USE moduleMesh3DCart + USE moduleMesh2DCyl + USE moduleMesh2DCart + USE moduleMesh1DRad + USE moduleMesh1DCart USE moduleErrors USE moduleOutput USE moduleRefParam @@ -873,6 +879,7 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile REAL(8):: volume + CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT object = 'geometry' @@ -957,7 +964,39 @@ MODULE moduleInput END SELECT - !Get the type of mesh + !Link the procedure to connect meshes + SELECT CASE(mesh%dimen) + CASE(3) + mesh%connectMesh => connectMesh3DCart + + CASE(2) + SELECT CASE(mesh%geometry) + CASE("Cyl") + mesh%connectMesh => connectMesh2DCyl + + CASE("Cart") + mesh%connectMesh => connectMesh2DCart + + END SELECT + + CASE(1) + SELECT CASE(mesh%geometry) + CASE("Rad") + mesh%connectMesh => connectMesh1DRad + + CASE("Cart") + mesh%connectMesh => connectMesh1DCart + + END SELECT + + END SELECT + + IF (doubleMesh) THEN + meshColl%connectMesh => mesh%connectMesh + + END IF + + !Get the format of mesh CALL config%get(object // '.meshType', meshFormat, found) SELECT CASE(meshFormat) CASE ("gmsh2") @@ -967,6 +1006,16 @@ MODULE moduleInput END IF + CASE ("vtu") + CALL initVTU(mesh) + IF (doubleMesh) THEN + CALL initVTU(meshColl) + + END IF + + CASE DEFAULT + CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry') + END SELECT !Reads the mesh file diff --git a/src/modules/makefile b/src/modules/makefile index 254d015..743144f 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -11,7 +11,7 @@ common.o: output.o: moduleSpecies.o common.o $(MAKE) -C output all -mesh.o: moduleCollisions.o moduleBoundary.o output.o +mesh.o: moduleCollisions.o moduleBoundary.o output.o common.o $(MAKE) -C mesh all solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index dfe9605..97ec729 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -42,13 +42,13 @@ MODULE moduleMeshOutput0D USE moduleOutput IMPLICIT NONE - CLASS(meshGeneric), INTENT(inout):: self + CLASS(meshGeneric), INTENT(in):: self INTEGER, INTENT(in):: t CHARACTER(:), ALLOCATABLE:: fileName INTEGER:: k fileName='OUTPUT_Collisions.dat' - IF (t == 0) THEN + IF (t == tInitial) THEN OPEN(20, file = path // folder // '/' // fileName, action = 'write') WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions" WRITE(*, "(6X,A15,A)") "Creating file: ", fileName diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index ae1fe05..0df1289 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -2,7 +2,7 @@ MODULE moduleMeshInputGmsh2 !Reads a mesh in the Gmsh v2.0 format CONTAINS - !Inits a mesh to use Gmsh2 format + !Init a mesh to use Gmsh2 format SUBROUTINE initGmsh2(self) USE moduleMesh USE moduleMeshOutputGmsh2 @@ -23,7 +23,7 @@ MODULE moduleMeshInputGmsh2 END SUBROUTINE initGmsh2 - !Reads a Gmsh 2 format + !Read a Gmsh 2 format SUBROUTINE readGmsh2(self, filename) USE moduleMesh3DCart USE moduleMesh2DCyl @@ -71,30 +71,26 @@ MODULE moduleMeshInputGmsh2 SELECT CASE(self%dimen) CASE(3) ALLOCATE(meshNode3Dcart::self%nodes(n)%obj) - self%connectMesh => connectMesh3DCart CASE(2) SELECT CASE(self%geometry) CASE("Cyl") ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj) - self%connectMesh => connectMesh2DCyl CASE("Cart") ALLOCATE(meshNode2DCart:: self%nodes(n)%obj) - self%connectMesh => connectMesh2DCart END SELECT + r(3) = 0.D0 CASE(1) SELECT CASE(self%geometry) CASE("Rad") ALLOCATE(meshNode1DRad:: self%nodes(n)%obj) - self%connectMesh => connectMesh1DRad CASE("Cart") ALLOCATE(meshNode1DCart:: self%nodes(n)%obj) - self%connectMesh => connectMesh1DCart END SELECT r(2:3) = 0.D0 @@ -111,7 +107,7 @@ MODULE moduleMeshInputGmsh2 !Reads total number of elements (no nodes) READ(10, *) totalNumElem - !conts edges and volume elements + !Count edges and volume elements SELECT TYPE(self) TYPE IS(meshParticles) self%numEdges = 0 @@ -151,7 +147,7 @@ MODULE moduleMeshInputGmsh2 END SELECT - !Allocates arrays + !Allocates array of cells ALLOCATE(self%cells(1:self%numCells)) SELECT TYPE(self) @@ -179,45 +175,30 @@ MODULE moduleMeshInputGmsh2 END SELECT CASE (2) + ALLOCATE(p(1:2)) + READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) + !Associate boundary condition procedure. + bt = getBoundaryId(boundaryType) + SELECT CASE(self%geometry) CASE("Cyl") - ALLOCATE(p(1:2)) - - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) - ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj) CASE("Cart") - ALLOCATE(p(1:2)) - - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) - ALLOCATE(meshEdge2DCart:: self%edges(e)%obj) END SELECT CASE(1) + ALLOCATE(p(1:1)) + READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) + !Associate boundary condition + bt = getBoundaryId(boundaryType) SELECT CASE(self%geometry) CASE("Rad") - ALLOCATE(p(1:1)) - - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - ALLOCATE(meshEdge1DRad:: self%edges(e)%obj) CASE("Cart") - ALLOCATE(p(1:1)) - - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - ALLOCATE(meshEdge1DCart:: self%edges(e)%obj) END SELECT @@ -231,9 +212,9 @@ MODULE moduleMeshInputGmsh2 END SELECT - !Read and initialize volumes + !Read and initialize cells DO e = 1, self%numCells - !Reads the volume according to the geometry + !Read the cell according to the geometry SELECT CASE(self%dimen) CASE(3) READ(10, *) n, elemType diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 53485f4..2217daf 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -85,6 +85,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleRefParam USE moduleSpecies USE moduleOutput + USE moduleMeshInoutCommon IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -93,13 +94,11 @@ MODULE moduleMeshOutputGmsh2 TYPE(outputFormat):: output(1:self%numNodes) REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tstring time = DBLE(t)*tauMin*ti_ref DO i = 1, nSpecies - WRITE(tstring, iterationFormat) t - fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh' + fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) @@ -142,16 +141,16 @@ MODULE moduleMeshOutputGmsh2 USE moduleCaseParam USE moduleCollisions USE moduleOutput + USE moduleMeshInoutCommon IMPLICIT NONE - CLASS(meshGeneric), INTENT(inout):: self + CLASS(meshGeneric), INTENT(in):: self INTEGER, INTENT(in):: t INTEGER:: numEdges INTEGER:: k, c INTEGER:: n REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tString CHARACTER(:), ALLOCATABLE:: title CHARACTER (LEN=2):: cString @@ -169,9 +168,8 @@ MODULE moduleMeshOutputGmsh2 IF (collOutput) THEN time = DBLE(t)*tauMin*ti_ref - WRITE(tString, iterationFormat) t - fileName='OUTPUT_' // tString// '_Collisions.msh' + fileName = formatFileName(prefix, 'Collisions', 'msh', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) @@ -203,6 +201,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleRefParam USE moduleCaseParam USE moduleOutput + USE moduleMeshInoutCommon IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -210,16 +209,14 @@ MODULE moduleMeshOutputGmsh2 INTEGER:: n, e REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tstring REAL(8):: Xi(1:3) Xi = (/ 0.D0, 0.D0, 0.D0 /) IF (emOutput) THEN time = DBLE(t)*tauMin*ti_ref - WRITE(tstring, iterationFormat) t - fileName='OUTPUT_' // tstring// '_EMField.msh' + fileName = formatFileName(prefix, 'Collisions', 'msh', t) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (20, file = path // folder // '/' // fileName) @@ -256,6 +253,7 @@ MODULE moduleMeshOutputGmsh2 USE moduleSpecies USE moduleOutput USE moduleAverage + USE moduleMeshInoutCommon IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self @@ -265,11 +263,11 @@ MODULE moduleMeshOutputGmsh2 INTEGER:: fileMean=10, fileDeviation=20 DO i = 1, nSpecies - fileName= 'Average_mean_' // species(i)%obj%name // '.msh' + fileName = formatFileName('Average_mean', species(i)%obj%name, 'msh') WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileMean, file = path // folder // '/' // fileName) - fileName= 'Average_deviation_' // species(i)%obj%name // '.msh' + fileName = formatFileName('Average_deviation', species(i)%obj%name, 'msh') WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (filedeviation, file = path // folder // '/' // fileName) diff --git a/src/modules/mesh/inout/makefile b/src/modules/mesh/inout/makefile index 2f73e73..a93161d 100644 --- a/src/modules/mesh/inout/makefile +++ b/src/modules/mesh/inout/makefile @@ -1,7 +1,13 @@ -all: gmsh2.o 0D.o +all: vtu.o gmsh2.o 0D.o + +vtu.o: moduleMeshInoutCommon.o + $(MAKE) -C vtu all gmsh2.o: $(MAKE) -C gmsh2 all 0D.o: $(MAKE) -C 0D all + +%.o: %.f90 + $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 new file mode 100644 index 0000000..e4a6c72 --- /dev/null +++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 @@ -0,0 +1,27 @@ +MODULE moduleMeshInoutCommon + + CHARACTER(LEN=4):: prefix = 'Step' + + CONTAINS + 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 + + +END MODULE moduleMeshInoutCommon diff --git a/src/modules/mesh/inout/vtu/makefile b/src/modules/mesh/inout/vtu/makefile new file mode 100644 index 0000000..07b471d --- /dev/null +++ b/src/modules/mesh/inout/vtu/makefile @@ -0,0 +1,7 @@ +all: moduleMeshInputVTU.o moduleMeshOutputVTU.o + +moduleMeshInputVTU.o: moduleMeshOutputVTU.o moduleMeshInputVTU.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +%.o: %.f90 + $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 new file mode 100644 index 0000000..763517f --- /dev/null +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -0,0 +1,552 @@ +MODULE moduleMeshInputVTU + !Reads mesh in the VTU format + + INTERFACE getValueFromLine + MODULE PROCEDURE getIntegerFromLine, getRealFromLine + + END INTERFACE + + INTERFACE readDataBlock + MODULE PROCEDURE readIntegerBlock, readRealBlock + + END INTERFACE + + CONTAINS + FUNCTION findLine(fileID, text) RESULT(line) + USE moduleErrors + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + CHARACTER(*):: text + CHARACTER(LEN=256):: line + INTEGER:: error, found + + error = 0 + found = 0 + !Reads all the file to find the 'text' string in a line + DO WHILE(error == 0 .AND. found == 0) + READ(fileID, "(A)", IOSTAT=error) line + found = INDEX(line, text) + IF (found > 0) THEN + EXIT + + END IF + END DO + + !If no line is found, return an error + IF (found == 0) THEN + CALL criticalError('String ' // text // ' not found in file.', 'findLine') + + END IF + + END FUNCTION findLine + + SUBROUTINE getIntegerFromLine(line, label, valueInteger) + USE moduleErrors + IMPLICIT NONE + + CHARACTER(LEN=256), INTENT(in):: line + CHARACTER(*), INTENT(in):: label + INTEGER, INTENT(out):: valueInteger + INTEGER:: labelStart, valueStart, valueEnd + + labelStart = 0 + labelStart = INDEX(line, label) + IF (labelStart == 0) THEN + CALL criticalError('Label '// label // ' not found in line', 'getValueFromLine') + + END IF + valueStart = INDEX(line(labelStart:256), '"') + labelStart + valueEnd = INDEX(line(valueStart:256), '"') + valueStart - 2 + READ(line(valueStart:valueEnd), '(I10)') valueInteger + + END SUBROUTINE getIntegerFromLine + + SUBROUTINE getRealFromLine(line, label, valueReal) + USE moduleErrors + IMPLICIT NONE + + CHARACTER(LEN=256), INTENT(in):: line + CHARACTER(*), INTENT(in):: label + REAL(8), INTENT(out):: valueReal + INTEGER:: labelStart, valueStart, valueEnd + + labelStart = 0 + labelStart = INDEX(line, label) + IF (labelStart == 0) THEN + CALL criticalError('Label '// label // ' not found in line', 'getValueFromLine') + + END IF + valueStart = INDEX(line(labelStart:256), '"') + labelStart + valueEnd = INDEX(line(valueStart:256), '"') + valueStart - 2 + READ(line(valueStart:valueEnd), '(F16.14)') valueReal + + END SUBROUTINE getRealFromLine + + SUBROUTINE readIntegerBlock(fileID, nData, array) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + INTEGER, INTENT(in):: nData + INTEGER, INTENT(out):: array(1:nData) + INTEGER:: iStart, iEnd, block + + iStart = 0 + iEnd = 0 + block = 6 !Assumes block of data in 6 columns + + DO WHILE (iStart < nData) + iStart = iStart + 1 + iEnd = iStart - 1 + block + IF (iEnd > nData) THEN + iEnd = nData + + END IF + READ(fileID, *) array(iStart:iEnd) + iStart = iEnd + + END DO + + END SUBROUTINE readIntegerBlock + + SUBROUTINE readRealBlock(fileID, nData, array) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + INTEGER, INTENT(in):: nData + REAL(8), INTENT(out):: array(1:nData) + INTEGER:: iStart, iEnd, block + + iStart = 0 + iEnd = 0 + block = 6 !Assumes block of data in 6 columns + + DO WHILE (iStart < nData) + iStart = iStart + 1 + iEnd = iStart - 1 + block + IF (iEnd > nData) THEN + iEnd = nData + + END IF + READ(fileID, *) array(iStart:iEnd) + iStart = iEnd + + END DO + + END SUBROUTINE readRealBlock + + SUBROUTINE initVTU(self) + USE moduleMesh + USE moduleMeshOutputVTU + IMPLICIT NONE + + CLASS(meshGeneric), INTENT(inout), TARGET:: self + + IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollVTU + SELECT TYPE(self) + TYPE IS(meshParticles) + self%printOutput => printOutputVTU + self%printEM => printEMVTU + self%readInitial => readInitialVTU + self%printAverage => printAverageVTU + + END SELECT + self%readMesh => readVTU + + END SUBROUTINE initVTU + + SUBROUTINE readVTU(self, filename) + USE moduleMesh3DCart + USE moduleMesh2DCyl + USE moduleMesh2DCart + USE moduleMesh1DRad + USE moduleMesh1DCart + USE moduleBoundary + IMPLICIT NONE + + CLASS(meshGeneric), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + REAL(8):: r(1:3) !3 generic coordinates + INTEGER:: fileID, error, found + CHARACTER(LEN=256):: line + INTEGER:: numNodes, numElements, numEdges + INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types + REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates + INTEGER:: n, e, c + INTEGER, ALLOCATABLE:: p(:) + INTEGER:: bt + + fileID = 10 + + OPEN(fileID, FILE=TRIM(filename)) + + !Find the number of nodes and elements (edges+cells) in the mesh. + line = findLine(fileID, '' + WRITE(fileID,"(2X, A)") '' + WRITE(fileID,"(4X, A,ES20.6E3,A)") '' + WRITE(fileID,"(6X, A, I10, A, I10, A)") '' + + END SUBROUTINE writeHeader + + SUBROUTINE writeFooter(fileID) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + + WRITE(fileID,"(6X, A)") '' + WRITE(fileID,"(4X, A)") '' + WRITE(fileID,"(2X, A)") '' + + END SUBROUTINE writeFooter + + FUNCTION getCellType(cell) RESULT(indexType) + USE moduleMesh3DCart + USE moduleMesh2DCyl + USE moduleMesh2DCart + USE moduleMesh1DRad + USE moduleMesh1DCart + USE moduleMesh0D + USE moduleErrors + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: cell + INTEGER:: indexType + + indexType = 0 + + SELECT TYPE(cell) + TYPE IS(meshCell3DCartTetra) + indexType = 10 + + TYPE IS(meshCell2DCylQuad) + indexType = 9 + + TYPE IS(meshCell2DCartQuad) + indexType = 9 + + TYPE IS(meshCell2DCylTria) + indexType = 5 + + TYPE IS(meshCell2DCartTria) + indexType = 5 + + TYPE IS(meshCell1DRadSegm) + indexType = 3 + + TYPE IS(meshCell1DCartSegm) + indexType = 3 + + TYPE IS(meshCell0D) + indexType = 1 + + CLASS DEFAULT + CALL criticalError('Cell not valid for VTU output', 'getCellType') + + END SELECT + + END FUNCTION getCellType + + SUBROUTINE writeMesh(self, fileID) + USE moduleMesh + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: fileID + INTEGER:: e + INTEGER:: offset(1:self%numCells), types(1:self%numCells) + + !Write nodes coordinates + WRITE(fileID, "(8X, A)") '' + WRITE(fileID, "(10X,A)") '' + WRITE(fileID, "(6(ES20.6E3))") (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes) + WRITE(fileID, "(10X, A)") '' + WRITE(fileID, "(8X, A)") '' + + WRITE(fileID, "(8X, A)") '' + !Write nodes connectivity of each cell + WRITE(fileID, "(10X,A)") '' + WRITE(fileID, "(6(I10))") (self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1, e = 1, self%numCells) !Array starts on 0 + WRITE(fileID, "(10X, A)") '' + !Write offset of each cell + offset(1) = self%cells(1)%obj%nNodes + WRITE(fileID, "(10X,A)") '' + DO e = 2, self%numCells + offset(e) = offset(e - 1) + self%cells(e)%obj%nNodes + + END DO + WRITE(fileID, "(6(I10))") offset + WRITE(fileID, "(10X, A)") '' + !Write type of each cell + WRITE(fileID, "(10X,A)") '' + DO e = 1, self%numCells + types(e) = getCellType(self%cells(e)%obj) + + END DO + WRITE(fileID, "(6(I10))") types + WRITE(fileID, "(10X, A)") '' + WRITE(fileID, "(8X, A)") '' + + END SUBROUTINE writeMesh + + 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) + INTEGER:: n + + 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 node data + WRITE(fileID,"(8X,A)") '' + !Write density + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(6(ES20.6E3))") (output(n)%density, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + !Write velocity + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(6(ES20.6E3))") (output(n)%velocity, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + !Write pressure + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(6(ES20.6E3))") (output(n)%pressure, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + !Write temperature + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(6(ES20.6E3))") (output(n)%temperature, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + !End node data + WRITE(fileID,"(8X,A)") '' + + 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 + + WRITE(fileID,"(8X,A)") '' + 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)") '' + WRITE(fileID, "(6(I10))") (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells) + WRITE(fileID, "(10X, A)") '' + + END DO + END DO + WRITE(fileID,"(8X,A)") '' + + END SUBROUTINE writeCollOutput + + SUBROUTINE writeEM(self, fileID) + USE moduleMesh + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: fileID + INTEGER:: n + REAL(8):: Xi(1:3) + + !Points in nodes + WRITE(fileID,"(8X,A)") '' + !Electric potential + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + !Magnetic Field + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes) + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(8X,A)") '' + + !Cell Data + Xi = (/ 0.D0, 0.D0, 0.D0 /) + WRITE(fileID,"(8X,A)") '' + !Electric field + WRITE(fileID,"(10X,A, A, A)") '' + WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) + WRITE(fileID,"(10X,A)") '' + WRITE(fileID,"(8X,A)") '' + + 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)") '' + WRITE (fileID + 1, "(2X, A)") '' + 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)") '' + + !Close collection file + IF (t == tFinal) THEN + WRITE (fileID + 1, "(2X, A)") '' + WRITE (fileID + 1, "(A)") '' + + END IF + CLOSE(fileID + 1) + + END SUBROUTINE writeCollection + + SUBROUTINE writeAverage(self, fileIDMean, fileIDDeviation, speciesIndex) + USE moduleMesh + USE moduleOutput + USE moduleAverage + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: fileIDMean, fileIDDeviation + INTEGER, INTENT(in):: speciesIndex + TYPE(outputFormat):: outputMean(1:self%numNodes) + TYPE(outputFormat):: outputDeviation(1:self%numNodes) + INTEGER:: n + + DO n = 1, self%numNodes + CALL calculateOutput(averageScheme(n)%mean%output(speciesIndex), outputMean(n), & + self%nodes(n)%obj%v, species(speciesIndex)%obj) + CALL calculateOutput(averageScheme(n)%deviation%output(speciesIndex), outputDeviation(n), & + self%nodes(n)%obj%v, species(speciesIndex)%obj) + + END DO + + !Write node data + WRITE(fileIDMean, "(8X,A)") '' + WRITE(fileIDDeviation,"(8X,A)") '' + !Write density + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%density, n = 1, self%numNodes) + WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%density, n = 1, self%numNodes) + WRITE(fileIDMean,"(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !Write velocity + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%velocity, n = 1, self%numNodes) + WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%velocity, n = 1, self%numNodes) + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !Write pressure + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%pressure, n = 1, self%numNodes) + WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%pressure, n = 1, self%numNodes) + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !Write temperature + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + WRITE(fileIDMean, "(6(ES20.6E3))") (outputMean(n)%temperature, n = 1, self%numNodes) + WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%temperature, n = 1, self%numNodes) + WRITE(fileIDMean, "(10X,A)") '' + WRITE(fileIDDeviation,"(10X,A)") '' + !End node data + WRITE(fileIDMean, "(8X,A)") '' + WRITE(fileIDDeviation,"(8X,A)") '' + + END SUBROUTINE writeAverage + + SUBROUTINE printOutputVTU(self,t) + USE moduleMesh + USE moduleSpecies + USE moduleMeshInoutCommon + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER, INTENT(in):: t + INTEGER:: n, i, fileID + CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection + TYPE(outputFormat):: output(1:self%numNodes) + + fileID = 60 + + DO i = 1, nSpecies + fileName = formatFileName(prefix, species(i)%obj%name, '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 writeSpeciesOutput(self, fileID, i) + + CALL writeFooter(fileID) + + CLOSE(fileID) + + !Write collection file for time plotting + fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd') + CALL writeCollection(fileID, t, fileName, filenameCollection) + + END DO + + END SUBROUTINE printOutputVTU + + SUBROUTINE printCollVTU(self,t) + USE moduleMesh + USE moduleOutput + USE moduleMeshInoutCommon + 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('Collection', 'Collisions', 'pvd') + CALL writeCollection(fileID, t, fileName, filenameCollection) + + END IF + + END SUBROUTINE printCollVTU + + SUBROUTINE printEMVTU(self, t) + USE moduleMesh + USE moduleMeshInoutCommon + 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('Collection', 'EMField', 'pvd') + CALL writeCollection(fileID, t, fileName, filenameCollection) + + END IF + + END SUBROUTINE printEMVTU + + SUBROUTINE printAverageVTU(self) + USE moduleMesh + USE moduleSpecies + USE moduleMeshInoutCommon + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER:: n, i, fileIDMean, fileIDDeviation + CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation + TYPE(outputFormat):: output(1:self%numNodes) + + fileIDMean = 66 + fileIDDeviation = 67 + + DO i = 1, nSpecies + fileNameMean = formatFileName('Average_mean', species(i)%obj%name, 'vtu') + WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameMean + OPEN (fileIDMean, file = path // folder // '/' // fileNameMean) + + fileNameDeviation = formatFileName('Average_deviation', species(i)%obj%name, 'vtu') + WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation + OPEN (fileIDDeviation, file = path // folder // '/' // fileNameDeviation) + + CALL writeHeader(self%numNodes, self%numCells, fileIDMean) + CALL writeHeader(self%numNodes, self%numCells, fileIDDeviation) + + CALL writeMesh(self, fileIDMean) + CALL writeMesh(self, fileIDDeviation) + + CALL writeAverage(self, fileIDMean, fileIDDeviation, i) + + CALL writeFooter(fileIDMean) + CALL writeFooter(fileIDDeviation) + + CLOSE(fileIDMean) + CLOSE(fileIDDeviation) + + END DO + + END SUBROUTINE printAverageVTU + +END MODULE moduleMeshOutputVTU diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 7482842..cb252e8 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -374,7 +374,7 @@ MODULE moduleMesh !Prints number of collisions in each cell SUBROUTINE printColl_interface(self, t) IMPORT meshGeneric - CLASS(meshGeneric), INTENT(inout):: self + CLASS(meshGeneric), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printColl_interface