From b2eb7c562281b361256063df4018e2aa63c743c8 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Wed, 14 Dec 2022 16:22:59 +0100 Subject: [PATCH] First commit for average scheme New module defined that will take care of averaging the output in the nodes. --- src/makefile | 2 +- src/modules/makefile | 2 +- src/modules/moduleAverage.f90 | 65 +++++++++++++++++++++++++++++++++++ src/modules/moduleMath.f90 | 11 ++++++ src/modules/moduleOutput.f90 | 64 ++++++++++++++++++++++++++++++---- src/modules/moduleSolver.f90 | 16 +++++++++ 6 files changed, 151 insertions(+), 9 deletions(-) create mode 100644 src/modules/moduleAverage.f90 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..d3931f4 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -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/moduleAverage.f90 b/src/modules/moduleAverage.f90 new file mode 100644 index 0000000..9b614c0 --- /dev/null +++ b/src/modules/moduleAverage.f90 @@ -0,0 +1,65 @@ +MODULE moduleAverage + USE moduleOutput + + TYPE:: averageData + TYPE(outputNode), ALLOCATABLE:: output(:) + TYPE(emNode):: emData + + END TYPE averageData + + !Generic type for average scheme + TYPE, PUBLIC:: averageGeneric + INTEGER:: tStart !Starting iteartion for average scheme + TYPE(averageData), ALLOCATABLE:: mean(:) + TYPE(averageData), ALLOCATABLE:: deviation(:) + CONTAINS + PROCEDURE, PASS:: updateAverage + + END TYPE averageGeneric + + TYPE(averageGeneric):: averageScheme + + !Logical to determine if average scheme must be used + LOGICAL:: useAverage + + CONTAINS + !Based on Welford's online algorithm + SUBROUTINE updateAverage(self, t) + USE moduleMesh + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(averageGeneric), INTENT(inout):: self + INTEGER, INTENT(in):: t + INTEGER:: tAverage + INTEGER:: n, s + TYPE(averageData):: newAverage + + tAverage = t - self%tStart + + IF (tAverage == 1) THEN + !First iteration in which average scheme is used + DO n = 1, mesh%numNodes + !Copy current data as mean + self%mean(n)%output(:) = mesh%nodes(n)%obj%output(:) + + END DO + + ELSEIF (tAverage > 1) THEN + !Normal average step + DO n = 1, mesh%numNodes + DO s = 1, nSpecies + newAverage%output(s) = self%mean(n)%output(s) + (mesh%nodes(n)%obj%output(s) - self%mean(n)%output(s))/tAverage + + END DO + + self%mean(n)%output(:) = newAverage%output(:) + + END DO + + END IF + + END SUBROUTINE updateAverage + +END MODULE moduleAverage 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..9cb23a7 100644 --- a/src/modules/moduleOutput.f90 +++ b/src/modules/moduleOutput.f90 @@ -1,9 +1,19 @@ !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:: copyOutputNode + PROCEDURE, PASS(self), PRIVATE:: addOutputNode + PROCEDURE, PASS(self), PRIVATE:: subOutputNode + PROCEDURE, PASS(self), PRIVATE:: divOutputNode_int + GENERIC, PUBLIC :: ASSIGNMENT(=) => copyOutputNode + GENERIC, PUBLIC :: OPERATOR(+) => addOutputNode + GENERIC, PUBLIC :: OPERATOR(-) => subOutputNode + GENERIC, PUBLIC :: OPERATOR(/) => divOutputNode_int END TYPE @@ -32,16 +42,56 @@ MODULE moduleOutput LOGICAL:: emOutput = .FALSE. CONTAINS - FUNCTION tensorTrace(a) RESULT(t) + PURE SUBROUTINE copyOutputNode(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 copyOutputNode + + PURE FUNCTION addOutputNode(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 addOutputNode + + PURE FUNCTION subOutputNode(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 subOutputNode + + PURE FUNCTION divOutputNode_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 divOutputNode_int SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn) USE moduleConstParam diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 017a321..7182184 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 @@ -818,5 +820,19 @@ MODULE moduleSolver END SUBROUTINE doOutput + SUBROUTINE doAverage(t) + USE moduleAverage + IMPLICIT NONE + + INTEGER, INTENT(in):: t + + IF (useAverage) THEN + CALL averageScheme%updateAverage(t) + + END IF + + END SUBROUTINE doAverage + + END MODULE moduleSolver