diff --git a/doc/user-manual/bibliography.bib b/doc/user-manual/bibliography.bib index 0ef1811..878b7ee 100644 --- a/doc/user-manual/bibliography.bib +++ b/doc/user-manual/bibliography.bib @@ -51,4 +51,15 @@ howpublished = {\url{https://gmsh.info/}}, } +@Article{welford1962note, + author = {Welford, BP}, + journal = {Technometrics}, + title = {Note on a method for calculating corrected sums of squares and products}, + year = {1962}, + number = {3}, + pages = {419--420}, + volume = {4}, + publisher = {Taylor \& Francis}, +} + @Comment{jabref-meta: databaseType:bibtex;} diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index a84376f..66aba8b 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 b41f809..e5c0d4f 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -1,4 +1,4 @@ -\documentclass[10pt,a4paper,oneside]{book} +\documentclass[10pt,a4paper,twoside]{book} \usepackage[latin1]{inputenc} \usepackage{amsmath} \usepackage{amsfonts} @@ -22,6 +22,7 @@ Author = {Jorge Gonzalez} } } +\usepackage[margin=1in]{geometry} % Reduces margins of the document % Allows breaking of URL in bibliography. \setcounter{biburllcpenalty}{7000} \setcounter{biburlucpenalty}{8000} @@ -260,6 +261,18 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Electromagnetic field} WIP. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + \section{Average scheme} + Particle-in-cell codes has an intrinsic statistical noise associated with them. + Although this can be reduced by increasing the number of particles, this also increases the CPU requirements of the case. + + It is quite common that most cases reach a quasi-steady state after a number of iterations and time-average results can be obtained after to improve analysis, plotting and restarting the case using these time-average results as new species backgrounds. + Although this is possible to do once the simulation is finished with post-processing tools, this is limited to the amount of iterations printed. + + \Gls{fpakc} implements a simple average scheme that, after a start time provided by the user, scores a mean and standard deviation of all the main species properties, and the electromagnetic field. + This scheme is based on the Welford's online algorithm~\cite{welford1962note}. + The averaged data is written in the same format as the input mesh at the end of the simulation. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \chapter{Installation} @@ -642,7 +655,7 @@ make \item \textbf{initialTime}: Real. Units of $\unit{s}$. Initial simulation time. - If no value is provided, the initial time is set to $\unit[0]{s}$. + If no value is provided, the initial time is set to $\unit[0.0]{s}$. \item \textbf{pusher}: Character. Array dimension 'number of species'. Indicates the type of pusher used for each species: @@ -684,6 +697,17 @@ make File must be located at \textbf{output.path}. \end{itemize} \end{itemize} +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + \subsection{average} + This object determines the use of an average scheme. + If this object exists in the input file, average will be written at the end of the simulation. + Acceptable values are: + \begin{itemize} + \item \textbf{startTime}: Real. + Units in $\unit{s}$. + Simulation physical time in which average scheme will start to compute the mean and standard validation. + If no value is provided, the initial time is set to $\unit[0.0]{s}$. + \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{interactions}\label{ssec:input_interactions} This object determine the different interactions among species. diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 03d0bbb..1224f53 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -114,6 +114,8 @@ PROGRAM fpakc !$OMP SINGLE tEMField = omp_get_wtime() - tEMField + CALL doAverage(t) + tStep = omp_get_wtime() - tStep !Output data diff --git a/src/makefile b/src/makefile index 9022723..24ad089 100644 --- a/src/makefile +++ b/src/makefile @@ -4,7 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \ - $(OBJDIR)/moduleProbe.o \ + $(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \ $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ $(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \ $(OBJDIR)/moduleMesh3DCart.o \ diff --git a/src/modules/makefile b/src/modules/makefile index c83d436..f1545dc 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -7,7 +7,7 @@ OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \ all: $(OBJS) -mesh.o: moduleCollisions.o moduleBoundary.o +mesh.o: moduleCollisions.o moduleBoundary.o moduleAverage.o $(MAKE) -C mesh all moduleCollisions.o: moduleList.o moduleMath.o moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90 @@ -25,7 +25,7 @@ moduleList.o: moduleConstParam.o moduleErrors.o moduleCaseParam.o moduleSpecies. moduleOutput.o: moduleMath.o moduleRefParam.o moduleOutput.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleSolver.o: moduleProbe.o moduleEM.o moduleSolver.f90 +moduleSolver.o: moduleProbe.o moduleEM.o moduleAverage.o moduleSolver.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ moduleProbe.o: mesh.o moduleProbe.f90 diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 57619ee..660a313 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -260,10 +260,17 @@ MODULE moduleMesh1DRad Xii = 0.D0 fPsi = self%fPsi(Xii) detJ = self%detJac(Xii) + !Computes total volume of the cell r = DOT_PRODUCT(fPsi, self%r) l = 2.D0*detJ - self%volume = r*l - self%arNodes = fPsi*r*l + self%volume = r*l + !Computes volume per node + xi = (/-5.D-1, 0.D0, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(1) = fPsi(1)*r*l + xi = (/ 5.D-1, 0.D0, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(2) = fPsi(2)*r*l END SUBROUTINE areaRad diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 0a6b7f4..c576c71 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -249,7 +249,7 @@ MODULE moduleMesh2DCyl dr = self%r(2) - self%r(1) dz = self%z(2) - self%z(1) IF (dr /= 0.D0) THEN - r(2) = dr*DSQRT(rnd) + self%r(1) + r(2) = dr * DSQRT(rnd) + self%r(1) r(1) = dz * (r(2) - self%r(1))/dr + self%z(1) ELSE @@ -318,9 +318,22 @@ MODULE moduleMesh2DCyl xi = 0.D0 detJ = self%detJac(xi)*PI8 !4*2*pi fPsi = self%fPsi(xi) + !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) - self%volume = r*detJ - self%arNodes = fPsi*r*detJ + self%volume = r*detJ + !Computes volume per node + xi = (/-5.D-1, -5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(1) = fPsi(1)*r*detJ + xi = (/ 5.D-1, -5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(2) = fPsi(2)*r*detJ + xi = (/ 5.D-1, 5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(3) = fPsi(3)*r*detJ + xi = (/-5.D-1, 5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(4) = fPsi(4)*r*detJ END SUBROUTINE areaQuad @@ -701,8 +714,10 @@ MODULE moduleMesh2DCyl xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) detJ = self%detJac(xi)*PI !2PI*1/2 fPsi = self%fPsi(xi) + !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) self%volume = r*detJ + !Computes volume per node self%arNodes = fPsi*r*detJ END SUBROUTINE areaTria diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 0a3a557..7832843 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -13,9 +13,10 @@ MODULE moduleMeshInputGmsh2 IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollGmsh2 SELECT TYPE(self) TYPE IS(meshParticles) - self%printOutput => printOutputGmsh2 - self%printEM => printEMGmsh2 - self%readInitial => readInitialGmsh2 + self%printOutput => printOutputGmsh2 + self%printEM => printEMGmsh2 + self%readInitial => readInitialGmsh2 + self%printAverage => printAverageGmsh2 END SELECT self%readMesh => readGmsh2 diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 40eda4f..7eade3e 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -1,6 +1,84 @@ MODULE moduleMeshOutputGmsh2 CONTAINS + !Header for mesh format + SUBROUTINE writeGmsh2HeaderMesh(fileID) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + + WRITE(fileID, "(A)") '$MeshFormat' + WRITE(fileID, "(A)") '2.2 0 8' + WRITE(fileID, "(A)") '$EndMeshFormat' + + END SUBROUTINE writeGmsh2HeaderMesh + + !Node data subroutines + !Header + SUBROUTINE writeGmsh2HeaderNodeData(fileID, title, iteration, time, dimensions, nNodes) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + CHARACTER(*), INTENT(in):: title + INTEGER, INTENT(in):: iteration, dimensions, nNodes + REAL(8), INTENT(in):: time + + + WRITE(fileID, "(A)") '$NodeData' + WRITE(fileID, "(I10)") 1 + WRITE(fileID, "(A1, A, A1)") '"' , title , '"' + WRITE(fileID, "(I10)") 1 + WRITE(fileID, "(ES20.6E3)") time + WRITE(fileID, "(I10)") 3 + WRITE(fileID, "(I10)") iteration + WRITE(fileID, "(I10)") dimensions + WRITE(fileID, "(I10)") nNodes + + END SUBROUTINE writeGmsh2HeaderNodeData + + !Footer + SUBROUTINE writeGmsh2FooterNodeData(fileID) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + + WRITE(fileID, "(A)") '$EndNodeData' + + END SUBROUTINE writeGmsh2FooterNodeData + + !Element data subroutines + !Header + SUBROUTINE writeGmsh2HeaderElementData(fileID, title, iteration, time, dimensions, nVols) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + CHARACTER(*), INTENT(in):: title + INTEGER, INTENT(in):: iteration, dimensions, nVols + REAL(8), INTENT(in):: time + + + WRITE(fileID, "(A)") '$ElementData' + WRITE(fileID, "(I10)") 1 + WRITE(fileID, "(A1, A, A1)") '"' , title , '"' + WRITE(fileID, "(I10)") 1 + WRITE(fileID, "(ES20.6E3)") time + WRITE(fileID, "(I10)") 3 + WRITE(fileID, "(I10)") iteration + WRITE(fileID, "(I10)") dimensions + WRITE(fileID, "(I10)") nVols + + END SUBROUTINE writeGmsh2HeaderElementData + + !Footer + SUBROUTINE writeGmsh2FooterElementData(fileID) + IMPLICIT NONE + + INTEGER, INTENT(in):: fileID + + WRITE(fileID, "(A)") '$EndElementData' + + END SUBROUTINE writeGmsh2FooterElementData + !Prints the scattered properties of particles into the nodes SUBROUTINE printOutputGmsh2(self, t) USE moduleMesh @@ -21,65 +99,36 @@ MODULE moduleMeshOutputGmsh2 DO i = 1, nSpecies WRITE(tstring, iterationFormat) t - fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh' + fileName= 'OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh' WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) - WRITE(60, "(A)") '$MeshFormat' - WRITE(60, "(A)") '2.2 0 8' - WRITE(60, "(A)") '$EndMeshFormat' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"' // species(i)%obj%name // ' density (m^-3)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes + + CALL writeGmsh2HeaderMesh(60) + + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes) DO n=1, self%numNodes CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) WRITE(60, "(I6,ES20.6E3)") n, output(n)%density END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"' // species(i)%obj%name // ' velocity (m s^-1)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 3 - WRITE(60, *) self%numNodes + CALL writeGmsh2FooterNodeData(60) + + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"' // species(i)%obj%name // ' pressure (Pa)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes + CALL writeGmsh2FooterNodeData(60) + + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"' // species(i)%obj%name // ' temperature (K)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes + CALL writeGmsh2FooterNodeData(60) + + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature END DO - WRITE(60, "(A)") '$EndNodeData' + CALL writeGmsh2FooterNodeData(60) CLOSE (60) END DO @@ -102,7 +151,9 @@ MODULE moduleMeshOutputGmsh2 INTEGER:: n REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tstring + CHARACTER (LEN=iterationDigits):: tString + CHARACTER(:), ALLOCATABLE:: title + CHARACTER (LEN=2):: cString SELECT TYPE(self) TYPE IS(meshParticles) @@ -118,30 +169,26 @@ MODULE moduleMeshOutputGmsh2 IF (collOutput) THEN time = DBLE(t)*tauMin*ti_ref - WRITE(tstring, iterationFormat) t + WRITE(tString, iterationFormat) t - fileName='OUTPUT_' // tstring// '_Collisions.msh' + fileName='OUTPUT_' // tString// '_Collisions.msh' WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) - WRITE(60, "(A)") '$MeshFormat' - WRITE(60, "(A)") '2.2 0 8' - WRITE(60, "(A)") '$EndMeshFormat' + + CALL writeGmsh2HeaderMesh(60) + DO k = 1, nCollPairs DO c = 1, interactionMatrix(k)%amount - WRITE(60, "(A)") '$ElementData' - WRITE(60, "(A)") '1' - WRITE(60, "(5A,I2)") '"Pair ', interactionMatrix(k)%sp_i%name, '-', interactionMatrix(k)%sp_j%name, ' collision ', c - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numVols + WRITE(cString, "(I2)") c + title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString + CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numVols) DO n=1, self%numVols WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c) END DO - WRITE(60, "(A)") '$EndElementData' + CALL writeGmsh2FooterElementData(60) + END DO + END DO CLOSE(60) @@ -175,50 +222,26 @@ MODULE moduleMeshOutputGmsh2 fileName='OUTPUT_' // tstring// '_EMField.msh' WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (20, file = path // folder // '/' // fileName) - WRITE(20, "(A)") '$MeshFormat' - WRITE(20, "(A)") '2.2 0 8' - WRITE(20, "(A)") '$EndMeshFormat' - WRITE(20, "(A)") '$NodeData' - WRITE(20, "(A)") '1' - WRITE(20, "(A)") '"Potential (V)"' - WRITE(20, *) 1 - WRITE(20, *) time - WRITE(20, *) 3 - WRITE(20, *) t - WRITE(20, *) 1 - WRITE(20, *) self%numNodes + + CALL writeGmsh2HeaderMesh(20) + + CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref END DO - WRITE(20, "(A)") '$EndNodeData' + CALL writeGmsh2FooterNodeData(20) - WRITE(20, "(A)") '$ElementData' - WRITE(20, "(A)") '1' - WRITE(20, "(A)") '"Electric Field (V m^-1)"' - WRITE(20, *) 1 - WRITE(20, *) time - WRITE(20, *) 3 - WRITE(20, *) t - WRITE(20, *) 3 - WRITE(20, *) self%numVols + CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numVols) DO e=1, self%numVols WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref END DO - WRITE(20, "(A)") '$EndElementData' + CALL writeGmsh2FooterElementData(20) - WRITE(20, "(A)") '$NodeData' - WRITE(20, "(A)") '1' - WRITE(20, "(A)") '"Magnetic Field (T)"' - WRITE(20, *) 1 - WRITE(20, *) time - WRITE(20, *) 3 - WRITE(20, *) t - WRITE(20, *) 3 - WRITE(20, *) self%numNodes + CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref END DO - WRITE(20, "(A)") '$EndNodeData' + CALL writeGmsh2FooterNodeData(20) CLOSE(20) @@ -226,4 +249,76 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printEMGmsh2 + !Prints the average properties of particles into the nodes + SUBROUTINE printAverageGmsh2(self) + USE moduleMesh + USE moduleRefParam + USE moduleSpecies + USE moduleOutput + USE moduleAverage + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in):: self + INTEGER:: n, i + TYPE(outputFormat), DIMENSION(1:self%numNodes):: outputMean, outputDeviation + CHARACTER(:), ALLOCATABLE:: fileName + INTEGER:: fileMean=10, fileDeviation=20 + + DO i = 1, nSpecies + fileName= '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' + WRITE(*, "(6X,A15,A)") "Creating file: ", fileName + OPEN (filedeviation, file = path // folder // '/' // fileName) + + CALL writeGmsh2HeaderMesh(fileMean) + CALL writeGmsh2HeaderMesh(fileDeviation) + + CALL writeGmsh2HeaderNodeData(fileMean, species(i)%obj%name // ' density, mean (m^-3)', 0, 0.D0, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' density, sd (m^-3)', 0, 0.D0, 1, self%numNodes) + DO n=1, self%numNodes + CALL calculateOutput(averageScheme(n)%mean%output(i), outputMean(n), self%nodes(n)%obj%v, species(i)%obj) + WRITE(fileMean, "(I6,ES20.6E3)") n, outputMean(n)%density + CALL calculateOutput(averageScheme(n)%deviation%output(i), outputDeviation(n), self%nodes(n)%obj%v, species(i)%obj) + WRITE(fileDeviation, "(I6,ES20.6E3)") n, outputDeviation(n)%density + END DO + CALL writeGmsh2FooterNodeData(fileMean) + CALL writeGmsh2FooterNodeData(fileDeviation) + + CALL writeGmsh2HeaderNodeData(fileMean, species(i)%obj%name // ' velocity, mean (m s^-1)', 0, 0.D0, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' velocity, sd (m s^-1)', 0, 0.D0, 3, self%numNodes) + DO n=1, self%numNodes + WRITE(fileMean, "(I6,3(ES20.6E3))") n, outputMean(n)%velocity + WRITE(fileDeviation, "(I6,3(ES20.6E3))") n, outputDeviation(n)%velocity + END DO + CALL writeGmsh2FooterNodeData(fileMean) + CALL writeGmsh2FooterNodeData(fileDeviation) + + CALL writeGmsh2HeaderNodeData(fileMean, species(i)%obj%name // ' Pressure, mean (Pa)', 0, 0.D0, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' Pressure, sd (Pa)', 0, 0.D0, 1, self%numNodes) + DO n=1, self%numNodes + WRITE(fileMean, "(I6,3(ES20.6E3))") n, outputMean(n)%pressure + WRITE(fileDeviation, "(I6,3(ES20.6E3))") n, outputDeviation(n)%pressure + END DO + CALL writeGmsh2FooterNodeData(fileMean) + CALL writeGmsh2FooterNodeData(fileDeviation) + + CALL writeGmsh2HeaderNodeData(fileMean, species(i)%obj%name // ' Temperature, mean (K)', 0, 0.D0, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(fileDeviation, species(i)%obj%name // ' Temperature, sd (K)', 0, 0.D0, 1, self%numNodes) + DO n=1, self%numNodes + WRITE(fileMean, "(I6,3(ES20.6E3))") n, outputMean(n)%temperature + WRITE(fileDeviation, "(I6,3(ES20.6E3))") n, outputDeviation(n)%temperature + END DO + CALL writeGmsh2FooterNodeData(fileMean) + CALL writeGmsh2FooterNodeData(fileDeviation) + + CLOSE (fileMean) + CLOSE(fileDeviation) + + END DO + + END SUBROUTINE printAverageGmsh2 + END MODULE moduleMeshOutputGmsh2 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index ebbbed4..05c4b85 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -289,7 +289,7 @@ MODULE moduleMesh CONTAINS PROCEDURE, PASS:: doCollisions - END TYPE + END TYPE meshGeneric ABSTRACT INTERFACE !Reads the mesh from a file @@ -338,9 +338,10 @@ MODULE moduleMesh REAL(8), ALLOCATABLE, DIMENSION(:,:):: K !Permutation matrix for P L U factorization INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV - PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() - PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() - PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL() + PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() + PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() + PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL() + PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL() CONTAINS PROCEDURE, PASS:: constructGlobalK @@ -373,6 +374,14 @@ MODULE moduleMesh END SUBROUTINE printEM_interface + !Prints average values + SUBROUTINE printAverage_interface(self) + IMPORT meshParticles + + CLASS(meshParticles), INTENT(in):: self + + END SUBROUTINE printAverage_interface + END INTERFACE diff --git a/src/modules/moduleAverage.f90 b/src/modules/moduleAverage.f90 new file mode 100644 index 0000000..8d705d5 --- /dev/null +++ b/src/modules/moduleAverage.f90 @@ -0,0 +1,72 @@ +MODULE moduleAverage + USE moduleOutput + + TYPE:: averageData + TYPE(outputNode), ALLOCATABLE:: output(:) + TYPE(emNode):: emData + + END TYPE averageData + + !Generic type for average scheme + TYPE, PUBLIC:: averageGeneric + TYPE(averageData):: mean + TYPE(averageData):: deviation + CONTAINS + PROCEDURE, PASS:: firstAverage + PROCEDURE, PASS:: updateAverage + + END TYPE averageGeneric + + TYPE(averageGeneric), ALLOCATABLE:: averageScheme(:) + + !Logical to determine if average scheme must be used + LOGICAL:: useAverage = .FALSE. + INTEGER:: tAverageStart = 1 !Starting iteartion for average scheme + + CONTAINS + PURE SUBROUTINE firstAverage(self, output) + USE moduleOutput + IMPLICIT NONE + + CLASS(averageGeneric), INTENT(inout):: self + TYPE(outputNode), INTENT(in):: output(:) + + !Copy current data as mean + self%mean%output(:) = output(:) + + !Set deviation at 0 + self%deviation%output(:) = 0.D0 + + END SUBROUTINE firstAverage + + + !Based on Welford's online algorithm + SUBROUTINE updateAverage(self, output, tAverage) + USE moduleOutput + USE moduleSpecies, ONLY: nSpecies + IMPLICIT NONE + + CLASS(averageGeneric), INTENT(inout):: self + TYPE(outputNode), INTENT(in):: output(:) + INTEGER, INTENT(in):: tAverage + TYPE(averageData):: newMean, newDeviation + + ALLOCATE(newMean%output(1:nSpecies)) + ALLOCATE(newDeviation%output(1:nSpecies)) + + newMean%output(:) = self%mean%output(:) + (output(:) - self%mean%output(:))/tAverage + + newDeviation%output(:) = self%deviation%output(:) + ((output(:) - self%mean%output(:)) * & + (output(:) - newMean%output(:)) - & + self%deviation%output(:))/tAverage + + self%mean%output(:) = newMean%output(:) + + self%deviation%output(:) = newDeviation%output(:) + + DEALLOCATE(newMean%output) + DEALLOCATE(newDeviation%output) + + END SUBROUTINE updateAverage + +END MODULE moduleAverage diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index de6cc65..348edf3 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -75,6 +75,11 @@ MODULE moduleInput CALL readInject(config) CALL checkStatus(config, "readInject") + !Read average scheme + CALL verboseError('Reading average scheme...') + CALL readAverage(config) + CALL checkStatus(config, "readAverage") + !Read parallel parameters CALL verboseError('Reading Parallel configuration...') CALL readParallel(config) @@ -1180,6 +1185,40 @@ MODULE moduleInput END SUBROUTINE readInject + SUBROUTINE readAverage(config) + USE moduleAverage + USE moduleCaseParam, ONLY: tauMin + USE moduleMesh, ONLY: mesh + USE moduleSpecies, ONLY: nSpecies + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + LOGICAL:: found + REAL(8):: tStart + INTEGER:: n + + CALL config%info('average', found) + IF (found) THEN + useAverage = .TRUE. + CALL config%get('average.startTime', tStart, found) + + IF (found) THEN + tAverageStart = INT(tStart / tauMin) + + END IF + + ALLOCATE(averageScheme(1:mesh%numNodes)) + + DO n = 1, mesh%numNodes + ALLOCATE(averageScheme(n)%mean%output(1:nSpecies)) + ALLOCATE(averageScheme(n)%deviation%output(1:nSpecies)) + + END DO + + END IF + + END SUBROUTINE readAverage + !Reads the velocity distribution functions for each inject SUBROUTINE readVelDistr(config, inj, object) USE moduleErrors diff --git a/src/modules/moduleMath.f90 b/src/modules/moduleMath.f90 index ca8d780..adb0d24 100644 --- a/src/modules/moduleMath.f90 +++ b/src/modules/moduleMath.f90 @@ -59,4 +59,15 @@ MODULE moduleMath END FUNCTION + FUNCTION tensorTrace(a) RESULT(t) + IMPLICIT NONE + + REAL(8), DIMENSION(1:3,1:3):: a + REAL(8):: t + + t = 0.D0 + t = a(1,1)+a(2,2)+a(3,3) + + END FUNCTION tensorTrace + END MODULE moduleMath diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90 index ecc5462..18fbb7f 100644 --- a/src/modules/moduleOutput.f90 +++ b/src/modules/moduleOutput.f90 @@ -1,9 +1,22 @@ !Contains information about output MODULE moduleOutput IMPLICIT NONE + !Output for each node - TYPE outputNode + TYPE, PUBLIC:: outputNode REAL(8):: den = 0.D0, mom(1:3) = 0.D0, tensorS(1:3,1:3) = 0.D0 + CONTAINS + PROCEDURE, PASS(self), PRIVATE:: outputNode_equal_outputNode + PROCEDURE, PASS(self), PRIVATE:: outputNode_equal_real + PROCEDURE, PASS(self), PRIVATE:: outputNode_add_outputNode + PROCEDURE, PASS(self), PRIVATE:: outputNode_sub_outputNode + PROCEDURE, PASS(self), PRIVATE:: outputNode_mul_outputNode + PROCEDURE, PASS(self), PRIVATE:: outputNode_div_int + GENERIC, PUBLIC :: ASSIGNMENT(=) => outputNode_equal_outputNode, outputNode_equal_real + GENERIC, PUBLIC :: OPERATOR(+) => outputNode_add_outputNode + GENERIC, PUBLIC :: OPERATOR(-) => outputNode_sub_outputNode + GENERIC, PUBLIC :: OPERATOR(*) => outputNode_mul_outputNode + GENERIC, PUBLIC :: OPERATOR(/) => outputNode_div_int END TYPE @@ -32,16 +45,81 @@ MODULE moduleOutput LOGICAL:: emOutput = .FALSE. CONTAINS - FUNCTION tensorTrace(a) RESULT(t) + PURE SUBROUTINE outputNode_equal_outputNode(self, from) IMPLICIT NONE - REAL(8), DIMENSION(1:3,1:3):: a - REAL(8):: t + CLASS(outputNode), INTENT(inout):: self + CLASS(outputNode), INTENT(in):: from - t = 0.D0 - t = a(1,1)+a(2,2)+a(3,3) + self%den = from%den + self%mom = from%mom + self%tensorS = from%tensorS - END FUNCTION tensorTrace + END SUBROUTINE outputNode_equal_outputNode + + PURE ELEMENTAL SUBROUTINE outputNode_equal_real(self, from) + IMPLICIT NONE + + CLASS(outputNode), INTENT(inout):: self + REAL(8), INTENT(in):: from + + self%den = from + self%mom = from + self%tensorS = from + + END SUBROUTINE outputNode_equal_real + + PURE ELEMENTAL FUNCTION outputNode_add_outputNode(self, that) RESULT(total) + IMPLICIT NONE + + CLASS(outputNode), INTENT(in):: self + CLASS(outputNode), INTENT(in):: that + TYPE(outputNode):: total + + total%den = self%den + that%den + total%mom = self%mom + that%mom + total%tensorS = self%tensorS + that%tensorS + + END FUNCTION outputNode_add_outputNode + + PURE ELEMENTAL FUNCTION outputNode_sub_outputNode(self, that) RESULT(total) + IMPLICIT NONE + + CLASS(outputNode), INTENT(in):: self + CLASS(outputNode), INTENT(in):: that + TYPE(outputNode):: total + + total%den = self%den - that%den + total%mom = self%mom - that%mom + total%tensorS = self%tensorS - that%tensorS + + END FUNCTION outputNode_sub_outputNode + + PURE ELEMENTAL FUNCTION outputNode_mul_outputNode(self, that) RESULT(total) + IMPLICIT NONE + + CLASS(outputNode), INTENT(in):: self + CLASS(outputNode), INTENT(in):: that + TYPE(outputNode):: total + + total%den = self%den * that%den + total%mom = self%mom * that%mom + total%tensorS = self%tensorS * that%tensorS + + END FUNCTION outputNode_mul_outputNode + + PURE ELEMENTAL FUNCTION outputNode_div_int(self, that) RESULT(total) + IMPLICIT NONE + + CLASS(outputNode), INTENT(in):: self + INTEGER, INTENT(in):: that + TYPE(outputNode):: total + + total%den = self%den / REAL(that) + total%mom = self%mom / REAL(that) + total%tensorS = self%tensorS / REAL(that) + + END FUNCTION outputNode_div_int SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn) USE moduleConstParam diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 017a321..d0b2d01 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -1,4 +1,5 @@ MODULE moduleSolver + USE moduleAverage !Generic type for pusher of particles TYPE, PUBLIC:: pusherGeneric @@ -15,6 +16,7 @@ MODULE moduleSolver !Generic type for solver TYPE, PUBLIC:: solverGeneric TYPE(pusherGeneric), ALLOCATABLE:: pusher(:) + TYPE(averageGeneric), ALLOCATABLE:: averageScheme PROCEDURE(solveEM_interface), POINTER, NOPASS:: solveEM => NULL() PROCEDURE(weightingScheme_interface), POINTER, NOPASS:: weightingScheme => NULL() CONTAINS @@ -658,24 +660,18 @@ MODULE moduleSolver TYPE(particle), INTENT(inout):: part CLASS(meshVol), POINTER, INTENT(in):: volOld CLASS(meshVol), POINTER, INTENT(inout):: volNew - REAL(8):: fractionVolume, fractionWeight - INTEGER:: nSplit + REAL(8):: fractionVolume, pSplit - !If particle has change cell, call Weighting scheme - IF (volOld%n /= volNew%n) THEN + !If particle changes volume to smaller cell + IF (volOld%volume > volNew%volume) THEN fractionVolume = volOld%volume/volNew%volume - part%weight = part%weight * fractionVolume + !Calculate probability of splitting particle + pSplit = 1.D0 - DEXP(-fractionVolume) - fractionWeight = part%weight / part%species%weight - - IF (fractionWeight >= 2.D0) THEN - nSplit = FLOOR(fractionWeight) - CALL splitParticle(part, nSplit, volNew) - - ELSEIF (part%weight < 1.D0) THEN - !Particle has lost statistical meaning and will be terminated - part%n_in = .FALSE. + IF (random() < pSplit THEN + !Split particle in two + CALL splitParticle(part, 2, volNew) END IF @@ -816,7 +812,44 @@ MODULE moduleSolver END IF + !Output average values + IF (useAverage .AND. t == tFinal) THEN + CALL mesh%printAverage() + + END IF + END SUBROUTINE doOutput + SUBROUTINE doAverage(t) + USE moduleAverage + USE moduleMesh + IMPLICIT NONE + + INTEGER, INTENT(in):: t + INTEGER:: tAverage, n + + + IF (useAverage) THEN + tAverage = t - tAverageStart + + IF (tAverage == 1) THEN + !First iteration in which average scheme is used + DO n = 1, mesh%numNodes + CALL averageScheme(n)%firstAverage(mesh%nodes(n)%obj%output) + + END DO + ELSEIF (tAverage > 1) THEN + DO n = 1, mesh%numNodes + CALL averageScheme(n)%updateAverage(mesh%nodes(n)%obj%output, tAverage) + + END DO + + END IF + + END IF + + END SUBROUTINE doAverage + + END MODULE moduleSolver