diff --git a/src/makefile b/src/makefile index 83cf8db..de27c82 100644 --- a/src/makefile +++ b/src/makefile @@ -1,7 +1,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/moduleCompTime.o \ $(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \ $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ - $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.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)/moduleAverage.o $(OBJDIR)/moduleCoulomb.o \ @@ -9,7 +9,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \ $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ $(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \ - $(OBJDIR)/moduleMeshInputText.o $(OBJDIR)/moduleMeshOutputText.o \ + $(OBJDIR)/moduleMeshInputText.o $(OBJDIR)/moduleMeshOutputText.o \ $(OBJDIR)/moduleMesh3DCart.o \ $(OBJDIR)/moduleMesh2DCyl.o \ $(OBJDIR)/moduleMesh2DCart.o \ diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 3682adf..e421639 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -8,7 +8,6 @@ MODULE moduleInput SUBROUTINE readConfig(inputFile) USE json_module USE moduleErrors - USE moduleBoundary USE moduleOutput USE moduleMesh IMPLICIT NONE @@ -793,7 +792,7 @@ MODULE moduleInput !Reads boundary conditions for the mesh SUBROUTINE readBoundary(config) - USE moduleBoundary + use moduleMesh USE moduleErrors USE moduleSpecies USE moduleRefParam @@ -817,22 +816,22 @@ MODULE moduleInput INTEGER:: nTypes CALL config%info('boundary', found, n_children = nBoundary) - ALLOCATE(boundary(1:nBoundary)) + ALLOCATE(boundaries(1:nBoundary)) DO i = 1, nBoundary WRITE(iString, '(i2)') i object = 'boundary(' // TRIM(iString) // ')' - boundary(i)%n = i - CALL config%get(object // '.name', boundary(i)%name, found) - CALL config%get(object // '.physicalSurface', boundary(i)%physicalSurface, found) + boundaries(i)%n = i + CALL config%get(object // '.name', boundaries(i)%name, found) + CALL config%get(object // '.physicalSurface', boundaries(i)%physicalSurface, found) CALL config%info(object // '.bTypes', found, n_children = nTypes) IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary') - ALLOCATE(boundary(i)%bTypes(1:nSpecies)) + ALLOCATE(boundaries(i)%bTypes(1:nSpecies)) DO s = 1, nSpecies WRITE(sString,'(i2)') s object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')' CALL config%get(object // '.type', bType, found) - associate(bound => boundary(i)%bTypes(s)%obj) + associate(bound => boundaries(i)%bTypes(s)%obj) SELECT CASE(bType) CASE('reflection') ALLOCATE(boundaryReflection:: bound) diff --git a/src/modules/makefile b/src/modules/makefile index 508bb56..bcb96f6 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -1,5 +1,5 @@ OBJS = common.o output.o mesh.o solver.o init.o \ - moduleBoundary.o moduleCollisions.o moduleInject.o \ + moduleCollisions.o moduleInject.o \ moduleList.o moduleProbe.o moduleCoulomb.o \ moduleSpecies.o @@ -11,7 +11,7 @@ common.o: output.o: moduleSpecies.o common.o $(MAKE) -C output all -mesh.o: moduleCollisions.o moduleCoulomb.o moduleBoundary.o output.o common.o +mesh.o: moduleList.o moduleSpecies.o moduleCollisions.o moduleCoulomb.o output.o common.o $(MAKE) -C mesh all solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o @@ -20,9 +20,6 @@ solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o init.o: common.o solver.o moduleInject.o $(MAKE) -C init all -moduleBoundary.o: common.o moduleBoundary.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - moduleCollisions.o: moduleList.o moduleSpecies.o common.o moduleCollisions.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index dd86109..f7e91e2 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -4,7 +4,6 @@ ! z == unused MODULE moduleMesh1DCart USE moduleMesh - USE moduleMeshBoundary IMPLICIT NONE REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) @@ -103,7 +102,6 @@ MODULE moduleMesh1DCart !Init edge element SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface) USE moduleSpecies - USE moduleBoundary USE moduleErrors IMPLICIT NONE @@ -128,7 +126,7 @@ MODULE moduleMesh1DCart self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index - self%boundary => boundary(bt) + self%boundary => boundaries(bt) ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 3e36711..983743e 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -4,7 +4,6 @@ ! z == unused MODULE moduleMesh1DRad USE moduleMesh - USE moduleMeshBoundary IMPLICIT NONE REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) @@ -103,7 +102,6 @@ MODULE moduleMesh1DRad !Init edge element SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface) USE moduleSpecies - USE moduleBoundary USE moduleErrors IMPLICIT NONE @@ -128,7 +126,7 @@ MODULE moduleMesh1DRad self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index - self%boundary => boundary(bt) + self%boundary => boundaries(bt) ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 092c0d6..5886e85 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -4,7 +4,6 @@ ! z == unused MODULE moduleMesh2DCart USE moduleMesh - USE moduleMeshBoundary IMPLICIT NONE !Values for Gauss integral @@ -143,7 +142,6 @@ MODULE moduleMesh2DCart !Init edge element SUBROUTINE initEdge2DCart(self, n, p, bt, physicalSurface) USE moduleSpecies - USE moduleBoundary USE moduleErrors IMPLICIT NONE @@ -172,7 +170,7 @@ MODULE moduleMesh2DCart self%normal = self%normal/NORM2(self%normal) !Boundary index - self%boundary => boundary(bt) + self%boundary => boundaries(bt) ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 109512e..9634204 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -4,7 +4,6 @@ ! z == theta (unused) MODULE moduleMesh2DCyl USE moduleMesh - USE moduleMeshBoundary IMPLICIT NONE !Values for Gauss integral @@ -143,7 +142,6 @@ MODULE moduleMesh2DCyl !Init edge element SUBROUTINE initEdge2DCyl(self, n, p, bt, physicalSurface) USE moduleSpecies - USE moduleBoundary USE moduleErrors USE moduleConstParam, ONLY: PI IMPLICIT NONE @@ -181,7 +179,7 @@ MODULE moduleMesh2DCyl self%normal = self%normal/NORM2(self%normal) !Boundary index - self%boundary => boundary(bt) + self%boundary => boundaries(bt) ALLOCATE(self%fboundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 0a4ba79..1cbcac2 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -4,7 +4,6 @@ ! z == z MODULE moduleMesh3DCart USE moduleMesh - USE moduleMeshBoundary IMPLICIT NONE TYPE, PUBLIC, EXTENDS(meshNode):: meshNode3DCart @@ -106,7 +105,6 @@ MODULE moduleMesh3DCart !Init surface element SUBROUTINE initEdge3DCartTria(self, n, p, bt, physicalSurface) USE moduleSpecies - USE moduleBoundary USE moduleErrors USE moduleMath USE moduleRefParam, ONLY: L_ref @@ -146,7 +144,7 @@ MODULE moduleMesh3DCart self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D !Boundary index - self%boundary => boundary(bt) + self%boundary => boundaries(bt) ALLOCATE(self%fBoundary(1:nSpecies)) !Assign functions to boundary DO s = 1, nSpecies diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 1ad6a36..17961c4 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -30,7 +30,6 @@ MODULE moduleMeshInputGmsh2 USE moduleMesh2DCart USE moduleMesh1DRad USE moduleMesh1DCart - USE moduleBoundary IMPLICIT NONE CLASS(meshGeneric), INTENT(inout):: self diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index c7b89b5..f7bd0ac 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -161,7 +161,6 @@ MODULE moduleMeshInputVTU USE moduleMesh2DCart USE moduleMesh1DRad USE moduleMesh1DCart - USE moduleBoundary IMPLICIT NONE CLASS(meshGeneric), INTENT(inout):: self diff --git a/src/modules/mesh/makefile b/src/modules/mesh/makefile index 7c66c3b..e21187f 100644 --- a/src/modules/mesh/makefile +++ b/src/modules/mesh/makefile @@ -1,4 +1,4 @@ -all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o +all: moduleMesh.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o 3DCart.o: moduleMesh.o $(MAKE) -C 3DCart all @@ -20,9 +20,8 @@ all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o moduleMesh.o: moduleMesh.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ - -moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90 - $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + $(FC) $(FCFLAGS) -c moduleMesh@elements.f90 -o $(OBJDIR)/$@ + $(FC) $(FCFLAGS) -c moduleMesh@boundary.f90 -o $(OBJDIR)/$@ inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o $(MAKE) -C inout all diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 8a88e78..da7a43d 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -2,10 +2,10 @@ MODULE moduleMesh USE moduleList USE moduleOutput - USE moduleBoundary USE moduleCollisions IMPLICIT NONE + ! Declarations for elements !Generic mesh element TYPE, PUBLIC, ABSTRACT:: meshElement !Index @@ -14,7 +14,7 @@ MODULE moduleMesh END TYPE meshElement - !Parent of Node element + !Generic Node element TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshNode !Node volume REAL(8):: v = 0.D0 @@ -32,6 +32,14 @@ MODULE moduleMesh END TYPE meshNode + interface + pure module subroutine resetOutput(self) + class(meshNode), intent(inout):: self + + end subroutine resetOutput + + end interface + ABSTRACT INTERFACE !Interface of init a node (3D generic coordinates) SUBROUTINE initNode_interface(self, n, r) @@ -104,6 +112,18 @@ MODULE moduleMesh END TYPE meshEdge + interface + pure module function gatherF_edge_scalar(self, Xi, nNodes, valNodes) result(f) + class(meshEdge), intent(in):: self + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8), intent(in):: valNodes(1:nNodes) + real(8):: f + + end function gatherF_edge_scalar + + end interface + ABSTRACT INTERFACE !Inits the edge parameters SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface) @@ -210,14 +230,63 @@ MODULE moduleMesh !Subroutine to find in which cell a particle is located PROCEDURE, PASS:: findCell !Gather value and spatial derivative on the nodes at position Xi - PROCEDURE, PASS, PRIVATE:: gatherF_cell_scalar - PROCEDURE, PASS, PRIVATE:: gatherF_cell_array - PROCEDURE, PASS, PRIVATE:: gatherDF_cell_scalar - GENERIC:: gatherF => gatherF_cell_scalar, gatherF_cell_array - GENERIC:: gatherDF => gatherDF_cell_scalar + procedure, pass, private:: gatherF_cell_scalar + procedure, pass, private:: gatherF_cell_array + procedure, pass, private:: gatherDF_cell_scalar + generic:: gatherF => gatherF_cell_scalar, gatherF_cell_array + generic:: gatherDF => gatherDF_cell_scalar END TYPE meshCell + interface + recursive module subroutine findCell(self, part, oldCell) + use moduleSpecies, only: particle + + class(meshCell), intent(inout):: self + class(particle), intent(inout), target:: part + class(meshCell), optional, intent(in):: oldCell + + end subroutine findCell + + module subroutine scatter(self, nNodes, part) + use moduleSpecies, only: particle + + class(meshCell), intent(inout):: self + integer, intent(in):: nNodes + class(particle), intent(in):: part + + end subroutine scatter + + pure module function gatherF_cell_scalar(self, Xi, nNodes, valNodes) result(f) + class(meshCell), intent(in):: self + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8), intent(in):: valNodes(1:nNodes) + real(8):: f + + end function gatherF_cell_scalar + + pure module function gatherF_cell_array(self, Xi, nNodes, valNodes) result(f) + class(meshCell), intent(in):: self + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8), intent(in):: valNodes(1:nNodes, 1:3) + real(8):: f(1:3) + + end function gatherF_cell_array + + pure module function gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df) + + class(meshCell), intent(in):: self + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8), intent(in):: valNodes(1:nNodes) + real(8):: df(1:3) + + end function gatherDF_cell_scalar + + end interface + ABSTRACT INTERFACE SUBROUTINE initCell_interface(self, n, p, nodes) IMPORT:: meshCell @@ -358,6 +427,21 @@ MODULE moduleMesh END TYPE meshGeneric + interface + module subroutine doCollisions(self) + class(meshGeneric), intent(inout), target:: self + + end subroutine doCollisions + + module function findCellBrute(self, r) result(nVol) + class(meshGeneric), intent(in):: self + real(8), dimension(1:3), intent(in):: r + integer:: nVol + + end function findCellBrute + + end interface + ABSTRACT INTERFACE !Reads the mesh from a file SUBROUTINE readMesh_interface(self, filename) @@ -405,12 +489,25 @@ MODULE moduleMesh PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL() CONTAINS - !GENERIC PROCEDURES - PROCEDURE, PASS:: constructGlobalK - PROCEDURE, PASS:: doCoulomb + !GENERIC PROCEDURES + PROCEDURE, PASS:: constructGlobalK + PROCEDURE, PASS:: doCoulomb END TYPE meshParticles + interface + pure module subroutine constructGlobalK(self) + class(meshParticles), intent(inout):: self + + end subroutine constructGlobalK + + module subroutine doCoulomb(self) + class(meshParticles), intent(in), target:: self + + end subroutine doCoulomb + + end interface + ABSTRACT INTERFACE !Prints Species data SUBROUTINE printOutput_interface(self) @@ -427,12 +524,6 @@ MODULE moduleMesh END SUBROUTINE printEM_interface !Perform Coulomb Scattering - SUBROUTINE doCoulomb_interface(self) - IMPORT meshParticles - CLASS(meshParticles), INTENT(inout):: self - - END SUBROUTINE doCoulomb_interface - !Prints average values SUBROUTINE printAverage_interface(self) IMPORT meshParticles @@ -476,6 +567,23 @@ MODULE moduleMesh !Procedure to find a cell for a particle in meshColl PROCEDURE(findCellColl_interface), POINTER:: findCellColl => NULL() + interface + module subroutine findCellCollMesh(part) + use moduleSpecies + implicit none + + type(particle), intent(inout):: part + + end subroutine findCellCollMesh + + module subroutine findCellSameMesh(part) + use moduleSpecies + + type(particle), intent(inout):: part + + end subroutine findCellSameMesh + end interface + ABSTRACT INTERFACE SUBROUTINE findCellColl_interface(part) USE moduleSpecies @@ -497,797 +605,106 @@ MODULE moduleMesh !Complete path for the two meshes CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle - CONTAINS - !Constructs the global K matrix - PURE SUBROUTINE constructGlobalK(self) - IMPLICIT NONE + ! Boundary Particle Definitions + !Generic type for boundaries + TYPE, PUBLIC:: boundaryGeneric + CONTAINS - CLASS(meshParticles), INTENT(inout):: self - INTEGER:: c - INTEGER, ALLOCATABLE:: n(:) - REAL(8), ALLOCATABLE:: localK(:,:) - INTEGER:: i, j + END TYPE boundaryGeneric - DO c = 1, self%numCells - associate(nNodes => self%cells(c)%obj%nNodes) - ALLOCATE(n(1:nNodes)) - ALLOCATE(localK(1:nNodes, 1:nNodes)) - n = self%cells(c)%obj%getNodes(nNodes) - localK = self%cells(c)%obj%elemK(nNodes) + !Reflecting boundary + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryReflection + CONTAINS - DO i = 1, nNodes - DO j = 1, nNodes - self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j) + END TYPE boundaryReflection - END DO + !Absorption boundary + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAbsorption + CONTAINS - END DO + END TYPE boundaryAbsorption - DEALLOCATE(n, localK) + !Transparent boundary + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryTransparent + CONTAINS - end associate + END TYPE boundaryTransparent - END DO + !Symmetry axis + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis + CONTAINS - END SUBROUTINE constructGlobalK + END TYPE boundaryAxis - !Reset the output of node - PURE SUBROUTINE resetOutput(self) - USE moduleSpecies - USE moduleOutput - IMPLICIT NONE + !Wall Temperature boundary + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature + !Thermal velocity of the wall: square root(Wall temperature X specific heat) + REAL(8):: vTh + CONTAINS - CLASS(meshNode), INTENT(inout):: self - INTEGER:: k + END TYPE boundaryWallTemperature - DO k = 1, nSpecies - self%output(k)%den = 0.D0 - self%output(k)%mom = 0.D0 - self%output(k)%tensorS = 0.D0 + !Ionization boundary + TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization + REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals. + CLASS(speciesGeneric), POINTER:: species !Ion species + CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron + TYPE(table1D):: crossSection + REAL(8):: effectiveTime + REAL(8):: eThreshold + REAL(8):: deltaV + CONTAINS - END DO + END TYPE boundaryIonization - END SUBROUTINE resetOutput + ! Ensures quasi-neutrality by changing the reflection coefficient + type, public, extends(boundaryGeneric):: boundaryQuasiNeutrality + real(8):: alpha ! Reflection parameter + integer, allocatable:: edges(:) !Array with edges - !Gather the value of valNodes (scalar) at position Xi - PURE FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f) - IMPLICIT NONE + end type boundaryQuasiNeutrality + !Wrapper for boundary types (one per species) + TYPE:: bTypesCont + CLASS(boundaryGeneric), ALLOCATABLE:: obj + CONTAINS - CLASS(meshCell), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in):: valNodes(1:nNodes) - REAL(8):: f - REAL(8):: fPsi(1:nNodes) + END TYPE bTypesCont - fPsi = self%fPsi(Xi, nNodes) - f = DOT_PRODUCT(fPsi, valNodes) + !Wrapper for boundary conditions + TYPE:: boundaryCont + INTEGER:: n = 0 + CHARACTER(:), ALLOCATABLE:: name + INTEGER:: physicalSurface = 0 !Physical surface as defined in the mesh file + CLASS(bTypesCont), ALLOCATABLE:: bTypes(:) !Array for boundary per species + CONTAINS - END FUNCTION gatherF_cell_scalar + END TYPE boundaryCont - ! Gather the value of valNodes at position Xi of an edge - pure function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f) - implicit none + interface + module subroutine pointBoundaryFunction(edge, s) + class(meshEdge), intent(inout):: edge + integer, intent(in):: s !Species index - class(meshEdge), intent(in):: self - real(8), intent(in):: Xi(1:3) - integer, intent(in):: nNodes - real(8), intent(in):: valNodes(1:nNodes) - real(8):: f - real(8):: fPsi(1:nNodes) + end subroutine pointBoundaryFunction - fPsi = self%fPsi(Xi, nNodes) - f = dot_product(fPsi, valNodes) + module function getBoundaryId(physicalSurface) result(id) + integer:: physicalSurface - end function gatherF_edge_scalar + end function getBoundaryId - !Gather the value of valNodes (array) at position Xi - PURE FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f) - IMPLICIT NONE + module subroutine reflection(edge, part) + use moduleSpecies - CLASS(meshCell), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in):: valNodes(1:nNodes, 1:3) - REAL(8):: f(1:3) - REAL(8):: fPsi(1:nNodes) + class(meshEdge), intent(inout):: edge + class(particle), intent(inout):: part - fPsi = self%fPsi(Xi, nNodes) - f = MATMUL(fPsi, valNodes) + end subroutine reflection - END FUNCTION gatherF_cell_array + end interface - !Gather the spatial derivative of valNodes (scalar) at position Xi - PURE FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df) - IMPLICIT NONE - - CLASS(meshCell), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in):: valNodes(1:nNodes) - REAL(8):: df(1:3) - REAL(8):: dPsi(1:3, 1:nNodes) - REAL(8):: pDer(1:3,1:3) - REAL(8):: dPsiR(1:3, 1:nNodes) - REAL(8):: invJ(1:3, 1:3), detJ - - dPsi = self%dPsi(Xi, nNodes) - pDer = self%partialDer(nNodes, dPsi) - detJ = self%detJac(pDer) - invJ = self%invJac(pDer) - dPsiR = MATMUL(invJ, dPsi)/detJ - df = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), & - DOT_PRODUCT(dPsiR(2,:), valNodes), & - DOT_PRODUCT(dPsiR(3,:), valNodes) /) - - END FUNCTION gatherDF_cell_scalar - - !Scatters particle properties into cell nodes - SUBROUTINE scatter(self, nNodes, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshCell), INTENT(inout):: self - INTEGER, INTENT(in):: nNodes - CLASS(particle), INTENT(in):: part - REAL(8):: fPsi(1:nNodes) - INTEGER:: cellNodes(1:nNodes) - REAL(8):: tensorS(1:3, 1:3) - INTEGER:: sp - INTEGER:: i - CLASS(meshNode), POINTER:: node - REAL(8):: pFraction !Particle fraction - - cellNodes = self%getNodes(nNodes) - fPsi = self%fPsi(part%Xi, nNodes) - - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - - DO i = 1, nNodes - node => mesh%nodes(cellNodes(i))%obj - pFraction = fPsi(i)*part%weight - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + pFraction - node%output(sp)%mom(:) = node%output(sp)%mom(:) + pFraction*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END DO - - END SUBROUTINE scatter - - !Find next cell for particle - RECURSIVE SUBROUTINE findCell(self, part, oldCell) - USE moduleSpecies - USE moduleErrors - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshCell), INTENT(inout):: self - CLASS(particle), INTENT(inout), TARGET:: part - CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell - REAL(8):: Xi(1:3) - CLASS(meshElement), POINTER:: neighbourElement - INTEGER:: sp - - Xi = self%phy2log(part%r) - !Checks if particle is inside 'self' cell - IF (self%inside(Xi)) THEN - part%cell = self%n - part%Xi = Xi - part%n_in = .TRUE. - !Assign particle to listPart_in - IF (listInCells) THEN - CALL OMP_SET_LOCK(self%lock) - sp = part%species%n - CALL self%listPart_in(sp)%add(part) - self%totalWeight(sp) = self%totalWeight(sp) + part%weight - CALL OMP_UNSET_LOCK(self%lock) - - END IF - - ELSE - !If not, searches for a neighbour and repeats the process. - CALL self%neighbourElement(Xi, neighbourElement) - !Defines the next step - SELECT TYPE(neighbourElement) - CLASS IS(meshCell) - !Particle moved to new cell, repeat find procedure - CALL neighbourElement%findCell(part, self) - - CLASS IS (meshEdge) - !Particle encountered a surface, apply boundary - CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part) - - !If particle is still inside the domain, call findCell - IF (part%n_in) THEN - IF(PRESENT(oldCell)) THEN - CALL self%findCell(part, oldCell) - - ELSE - CALL self%findCell(part) - - END IF - END IF - - CLASS DEFAULT - WRITE (*, "(A, I6)") "Element = ", self%n - CALL criticalError("No connectivity found for element", "findCell") - - END SELECT - - END IF - - END SUBROUTINE findCell - - !If Coll and Particle are the same, simply copy the part%cell into part%cellColl - SUBROUTINE findCellSameMesh(part) - USE moduleSpecies - IMPLICIT NONE - - TYPE(particle), INTENT(inout):: part - - part%cellColl = part%cell - - END SUBROUTINE findCellSameMesh - - !TODO: try to combine this with the findCell for a regular mesh - !Find the volume in which particle reside in the mesh for collisions - !No boundary interaction taken into account - SUBROUTINE findCellCollMesh(part) - USE moduleSpecies - IMPLICIT NONE - - TYPE(particle), INTENT(inout):: part - LOGICAL:: found - CLASS(meshCell), POINTER:: cell - REAL(8), DIMENSION(1:3):: Xi - CLASS(meshElement), POINTER:: neighbourElement - INTEGER:: sp - - found = .FALSE. - - cell => meshColl%cells(part%cellColl)%obj - DO WHILE(.NOT. found) - Xi = cell%phy2log(part%r) - IF (cell%inside(Xi)) THEN - part%cellColl = cell%n - IF (listInCells) THEN - CALL OMP_SET_LOCK(cell%lock) - sp = part%species%n - CALL cell%listPart_in(sp)%add(part) - cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight - CALL OMP_UNSET_LOCK(cell%lock) - - END IF - found = .TRUE. - - ELSE - CALL cell%neighbourElement(Xi, neighbourElement) - SELECT TYPE(neighbourElement) - CLASS IS(meshCell) - !Try next element - cell => neighbourElement - - CLASS DEFAULT - !Should never happend, but just in case, stops loops - found = .TRUE. - - END SELECT - - END IF - - END DO - - END SUBROUTINE findCellCollMesh - - !Returns index of volume associated to a position (if any) - !If no voulme is found, returns 0 - !WARNING: This function is slow and should only be used in initialization phase - FUNCTION findCellBrute(self, r) RESULT(nVol) - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshGeneric), INTENT(in):: self - REAL(8), DIMENSION(1:3), INTENT(in):: r - INTEGER:: nVol - INTEGER:: e - REAL(8), DIMENSION(1:3):: Xi - - !Inits RESULT - nVol = 0 - - DO e = 1, self%numCells - Xi = self%cells(e)%obj%phy2log(r) - IF(self%cells(e)%obj%inside(Xi)) THEN - nVol = self%cells(e)%obj%n - EXIT - - END IF - - END DO - - END FUNCTION findCellBrute - - !Computes collisions in element - SUBROUTINE doCollisions(self) - USE moduleCollisions - USE moduleSpecies - USE moduleList - use moduleRefParam - USE moduleRandom - USE moduleOutput - USE moduleMath - USE moduleCaseParam, ONLY: timeStep - IMPLICIT NONE - - CLASS(meshGeneric), INTENT(inout), TARGET:: self - INTEGER:: e - CLASS(meshCell), POINTER:: cell - INTEGER:: k, i, j - INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell - REAL(8):: pMax !Maximum probability of collision - INTEGER:: nColl - TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:) - TYPE(particle), POINTER:: part_i, part_j - INTEGER:: n, c - REAL(8):: vRel, rMass, eRel - REAL(8):: sigmaVrelTotal - REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:) - REAL(8):: rnd_real !Random number for collision - INTEGER:: rnd_int !Random number for collision - - IF (MOD(timeStep, everyColl) == 0) THEN - !Collisions need to be performed in this iteration - !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) - DO e=1, self%numCells - - cell => self%cells(e)%obj - - !TODO: Simplify this, to many sublevels - !Iterate over the number of pairs - DO k = 1, nCollPairs - !Reset tally of collisions - IF (collOutput) THEN - cell%tallyColl(k)%tally = 0 - - END IF - - IF (interactionMatrix(k)%amount > 0) THEN - !Select the species for the collision pair - i = interactionMatrix(k)%sp_i%n - j = interactionMatrix(k)%sp_j%n - - !Number of particles per species in the collision pair - nPart_i = cell%listPart_in(i)%amount - nPart_j = cell%listPart_in(j)%amount - - IF (nPart_i > 0 .AND. nPart_j > 0) THEN - !Total number of particles for the collision pair - nPart = nPart_i + nPart_j - - !Resets the number of collisions in the cell - nColl = 0 - - !Probability of collision for pair i-j - pMax = (cell%totalWeight(i) + cell%totalWeight(j))*cell%sigmaVrelMax(k)*tauColl/cell%volume - - !Number of collisions in the cell - nColl = NINT(REAL(nPart)*pMax*0.5D0) - - !Converts the list of particles to an array for easy access - IF (nColl > 0) THEN - partTemp_i = cell%listPart_in(i)%convert2Array() - partTemp_j = cell%listPart_in(j)%convert2Array() - - END IF - - DO n = 1, nColl - !Select random particles - part_i => NULL() - part_j => NULL() - rnd_int = random(1, nPart_i) - part_i => partTemp_i(rnd_int)%part - rnd_int = random(1, nPart_j) - part_j => partTemp_j(rnd_int)%part - !If they are the same particle, skip - !TODO: Maybe try to improve this - IF (ASSOCIATED(part_i, part_j)) THEN - CYCLE - - END IF - - !If particles do not belong to the species, skip collision - !This can happen, for example, if particle has been previously ionized or removed - !TODO: Try to find a way to not lose these collisions. Maybe check new 'k' and use that for the collision? - IF (part_i%species%n /= i .OR. & - part_j%species%n /= j) THEN - CYCLE - - END IF - !Obtain the cross sections for the different processes - !TODO: From here it might be a procedure in interactionMatrix - vRel = NORM2(part_i%v-part_j%v) - rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m) - eRel = rMass*vRel**2 - CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) - - !Update maximum sigma*v_rel - IF (sigmaVrelTotal > cell%sigmaVrelMax(k)) THEN - cell%sigmaVrelMax(k) = sigmaVrelTotal - - END IF - - ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount)) - probabilityColl = 0.0 - DO c = 1, interactionMatrix(k)%amount - probabilityColl(c) = sigmaVrel(c)/cell%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1)) - - END DO - - !Selects random number between 0 and 1 - rnd_real = random() - - !If the random number is below the total probability of collision, collide particles - IF (rnd_real < sigmaVrelTotal / cell%sigmaVrelMax(k)) THEN - - !Loop over collisions - DO c = 1, interactionMatrix(k)%amount - IF (rnd_real <= probabilityColl(c)) THEN - !$OMP CRITICAL - CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel) - !$OMP END CRITICAL - - !If collisions are gonna be output, count the collision - IF (collOutput) THEN - cell%tallyColl(k)%tally(c) = cell%tallyColl(k)%tally(c) + 1 - - END IF - - !A collision has ocurred, exit the loop - EXIT - - END IF - - END DO - - END IF - - !Deallocate arrays for next collision - DEALLOCATE(sigmaVrel, probabilityColl) - - !End loop collisions in cell - END DO - - END IF - - END IF - - !End loop collision pairs - END DO - - !End loop volumes - END DO - !$OMP END DO - - END IF - - END SUBROUTINE doCollisions - - SUBROUTINE doCoulomb(self) - USE moduleCoulomb - USE moduleRandom - USE moduleOutput - USE moduleList - USE moduleMath - USE moduleRefParam - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshParticles), INTENT(in), TARGET:: self - CLASS(meshCell), POINTER:: cell - TYPE(interactionsCoulomb):: pair - INTEGER:: e - INTEGER:: k - INTEGER:: i, j - INTEGER:: n - INTEGER:: p - TYPE(lNode), POINTER:: partTemp - INTEGER(8), ALLOCATABLE:: cellNodes(:) - CLASS(meshNode), POINTER:: node - TYPE(outputFormat):: output - REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node - REAL(8):: density, velocity(1:3), temperature!values at particle position - REAL(8):: C(1:3), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference - REAL(8):: l, lW, l2 - REAL(8):: GlW, HlW - REAL(8):: normC - REAL(8):: cosThe, sinThe - REAL(8):: cosPhi, sinPhi - REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame - REAL(8):: A, AW - REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W - REAL(8):: theta_per !Random angle for perpendicular direction - REAL(8):: eps = 1.D-12 - REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ij, p_ij - REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ij - REAL(8):: massSum_ij - REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ji, p_ji - REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ji - REAL(8):: massSum_ji - REAL(8):: alpha_num, alpha_den, alpha, beta(1:3) - - - !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) - DO e = 1, self%numCells - cell => self%cells(e)%obj - cellNodes = cell%getNodes(cell%nNodes) - - ALLOCATE(densityNodes(1:cell%nNodes), & - velocityNodes(1:cell%nNodes, 1:3), & - temperatureNodes(1:cell%nNodes)) - - DO k=1, nCoulombPairs - pair = coulombMatrix(k) - i = pair%sp_i%n - j = pair%sp_j%n - - !Do scattering of particles from species_i due to species j - !Compute background properties of species_j - DO n = 1, cell%nNodes - node => self%nodes(cellNodes(n))%obj - CALL calculateOutput(node%output(j), output, node%v, pair%sp_j) - densityNodes(n) = output%density/n_ref - velocityNodes(n,1:3) = output%velocity(1:3)/v_ref - temperatureNodes(n) = output%temperature/T_ref - - END DO - - ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3)) - ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3)) - ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount)) - deltaV_ij = 0.D0 - p_ij = 0.D0 - mass_ij = 0.D0 - !Loop over particles of species_i - partTemp => cell%listPart_in(i)%head - p = 1 - DO WHILE(ASSOCIATED(partTemp)) - density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) - velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) - temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) - - !If cell temperature is too low, skip particle to avoid division by zero - IF (temperature>eps) THEN - l2 = pair%l2_j/temperature - l = SQRT(l2) - - ELSE - partTemp => partTemp%next - - CYCLE - - END IF - - A = pair%A_i*density - - C = partTemp%part%v - velocity - normC = NORM2(C) - - !C_3 = z; C_1, C2 = x, y (per) - C_per = NORM2(C(1:2)) - cosPhi = C(1) / C_per - sinPhi = C(2) / C_per - cosThe = C(3) / normC - sinThe = C_per / normC - - !Rotation matrix to go from W to C - rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column - -sinPhi, cosPhi, 0.D0, & !Second column - sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column - (/ 3, 3 /)) - - !W at start is = (0, 0, normC), so normW = normC - lW = l * normC - GlW = G(lW) - HlW = H(lW) - AW = A / normC - - !Calculate changes in W due to collision process - deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin - deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() - deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() - - !Random angle to distribute perpendicular change in velocity - theta_per = PI2*random() - - !Change W - W(1) = deltaW_per_square * COS(theta_per) - W(2) = deltaW_per_square * SIN(theta_per) - W(3) = normC + deltaW_par + deltaW_par_square - - !Compute changes in velocity for each particle - deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v - mass_ij(p) = pair%sp_i%m*partTemp%part%weight - p_ij(p,1:3) = mass_ij(p)*partTemp%part%v - - !Move to the next particle in the list - partTemp => partTemp%next - p = p + 1 - - END DO - - !Do corresponding collisions - IF (i /= j) THEN - !Do scattering of particles from species_j due to species i - !Compute background properties of species_i - DO n = 1, cell%nNodes - node => self%nodes(cellNodes(n))%obj - CALL calculateOutput(node%output(i), output, node%v, pair%sp_i) - densityNodes(n) = output%density/n_ref - velocityNodes(n,1:3) = output%velocity(1:3)/v_ref - temperatureNodes(n) = output%temperature/T_ref - - END DO - - ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3)) - ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3)) - ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount)) - deltaV_ji = 0.D0 - p_ji = 0.D0 - mass_ji = 0.D0 - !Loop over particles of species_j - partTemp => cell%listPart_in(j)%head - p = 1 - DO WHILE(ASSOCIATED(partTemp)) - density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) - velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) - temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) - - !If cell temperature is too low, skip particle to avoid division by zero - IF (temperature>eps) THEN - l2 = pair%l2_i/temperature - l = SQRT(l2) - - ELSE - partTemp => partTemp%next - - CYCLE - - END IF - A = pair%A_j*density - - C = partTemp%part%v - velocity - normC = NORM2(C) - - !C_3 = z; C_1, C2 = x, y (per) - C_per = NORM2(C(1:2)) - cosPhi = C(1) / C_per - sinPhi = C(2) / C_per - cosThe = C(3) / normC - sinThe = C_per / normC - - !Rotation matrix to go from W to C - rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column - -sinPhi, cosPhi, 0.D0, & !Second column - sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column - (/ 3, 3 /)) - - !W at start is = (0, 0, normC), so normW = normC - lW = l * normC - GlW = G(lW) - HlW = H(lW) - AW = A / normC - - !Calculate changes in W due to collision process - deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin - deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() - deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() - - !Random angle to distribute perpendicular change in velocity - theta_per = PI2*random() - - !Change W - W(1) = deltaW_per_square * COS(theta_per) - W(2) = deltaW_per_square * SIN(theta_per) - W(3) = normC + deltaW_par + deltaW_par_square - - !Compute changes in velocity for each particle - deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v - mass_ji(p) = pair%sp_j%m*partTemp%part%weight - p_ji(p,1:3) = mass_ji(p)*partTemp%part%v - - !Move to the next particle in the list - partTemp => partTemp%next - p = p + 1 - - END DO - - END IF - - !Calculate correction - !Total mass - massSum_ij = SUM(mass_ij) - massSum_ji = 0.D0 - - !Beta - beta = 0.D0 - DO p = 1, cell%listPart_in(i)%amount - beta = beta + mass_ij(p) * deltaV_ij(p,1:3) - - END DO - - IF (i /= j) THEN - massSum_ji = SUM(mass_ji) - DO p = 1, cell%listPart_in(j)%amount - beta = beta + mass_ji(p) * deltaV_ji(p,1:3) - - END DO - - END IF - - beta = beta / (massSum_ij + massSum_ji) - - !Alpha - alpha_num = 0.D0 - alpha_den = 0.D0 - DO p =1, cell%listPart_in(i)%amount - alpha_num = alpha_num + DOT_PRODUCT(p_ij(p,1:3), deltav_ij(p,1:3) - beta(1:3)) - alpha_den = alpha_den + mass_ij(p) * NORM2(deltav_ij(p,1:3) - beta(1:3))**2 - - END DO - - IF (i /= j) THEN - DO p = 1, cell%listPart_in(j)%amount - alpha_num = alpha_num + DOT_PRODUCT(p_ji(p,1:3), deltav_ji(p,1:3) - beta(1:3)) - alpha_den = alpha_den + mass_ji(p) * NORM2(deltav_ji(p,1:3) - beta(1:3))**2 - - END DO - - END IF - - alpha = -2.D0*alpha_num / alpha_den - - !Apply correction to particles velocity - partTemp => cell%listPart_in(i)%head - p = 1 - DO WHILE(ASSOCIATED(partTemp)) - partTemp%part%v = partTemp%part%v + alpha * (deltaV_ij(p,1:3) - beta(1:3)) - partTemp => partTemp%next - p = p + 1 - - END DO - - IF (i /= j) THEN - partTemp => cell%listPart_in(j)%head - p = 1 - DO WHILE(ASSOCIATED(partTemp)) - partTemp%part%v = partTemp%part%v + alpha * (deltaV_ji(p,1:3) - beta(1:3)) - partTemp => partTemp%next - p = p + 1 - - END DO - - END IF - - DEALLOCATE(deltaV_ij, p_ij, mass_ij) - - IF (i /= j) THEN - DEALLOCATE(deltaV_ji, p_ji, mass_ji) - - END IF - - END DO - - DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes) - - END DO - !$OMP END DO - - END SUBROUTINE doCoulomb + !Number of boundaries + INTEGER:: nBoundary = 0 + !Array for boundaries + TYPE(boundaryCont), ALLOCATABLE, TARGET:: boundaries(:) END MODULE moduleMesh diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMesh@boundary.f90 similarity index 70% rename from src/modules/mesh/moduleMeshBoundary.f90 rename to src/modules/mesh/moduleMesh@boundary.f90 index 6b8be87..0a9cd9a 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMesh@boundary.f90 @@ -1,9 +1,102 @@ !moduleMeshBoundary: Boundary functions for the mesh edges -MODULE moduleMeshBoundary - USE moduleMesh - +submodule(moduleMesh) boundary CONTAINS - SUBROUTINE reflection(edge, part) + FUNCTION getBoundaryId(physicalSurface) RESULT(id) + IMPLICIT NONE + + INTEGER:: physicalSurface + INTEGER:: id + INTEGER:: i + + id = 0 + DO i = 1, nBoundary + IF (physicalSurface == boundaries(i)%physicalSurface) id = boundaries(i)%n + + END DO + + END FUNCTION getBoundaryId + + SUBROUTINE initWallTemperature(boundary, T, c) + USE moduleRefParam + IMPLICIT NONE + + CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary + REAL(8), INTENT(in):: T, c !Wall temperature and specific heat + REAL(8):: vTh + + vTh = DSQRT(c * T) / v_ref + boundary = boundaryWallTemperature(vTh = vTh) + + END SUBROUTINE initWallTemperature + + SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) + USE moduleRefParam + USE moduleSpecies + USE moduleCaseParam + USE moduleConstParam + USE moduleErrors + IMPLICIT NONE + + CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary + REAL(8), INTENT(in):: me !Electron mass + REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties + INTEGER, INTENT(in):: ion + INTEGER, OPTIONAL, INTENT(in):: electronSecondary + REAL(8):: effTime + CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection + REAL(8), INTENT(in):: eThreshold + + ALLOCATE(boundaryIonization:: boundary) + + SELECT TYPE(boundary) + TYPE IS(boundaryIonization) + boundary%m0 = m0 / m_ref + boundary%n0 = n0 * Vol_ref + boundary%v0 = v0 / v_ref + boundary%vTh = DSQRT(kb*T0/m0)/v_ref + boundary%species => species(ion)%obj + IF (PRESENT(electronSecondary)) THEN + SELECT TYPE(sp => species(electronSecondary)%obj) + TYPE IS(speciesCharged) + boundary%electronSecondary => sp + + CLASS DEFAULT + CALL criticalError("Species " // sp%name // " chosen for " // & + "secondary electron is not a charged species", 'initIonization') + + END SELECT + + ELSE + boundary%electronSecondary => NULL() + + END IF + boundary%effectiveTime = effTime / ti_ref + CALL boundary%crossSection%init(crossSection) + CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) + boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2) + boundary%deltaV = DSQRT(boundary%eThreshold/me) + + END SELECT + + END SUBROUTINE initIonization + + subroutine initQuasiNeutrality(boundary) + implicit none + + class(boundaryGeneric), allocatable, intent(out):: boundary + integer:: e, et + + allocate(boundaryQuasiNeutrality:: boundary) + + select type(boundary) + type is(boundaryQuasiNeutrality) + boundary%alpha = 0.d0 + + end select + + end subroutine initQuasiNeutrality + + module SUBROUTINE reflection(edge, part) USE moduleCaseParam USE moduleSpecies IMPLICIT NONE @@ -224,18 +317,23 @@ MODULE moduleMeshBoundary class(meshEdge), intent(inout):: edge class(particle), intent(inout):: part real(8), allocatable:: density(:) - integer:: nNodes - integer, allocatable:: nodes(:) + class(meshCell), pointer:: cell + real(8):: EF_dir + real(8):: alpha select type(bound => edge%boundary%bTypes(part%species%n)%obj) type is(boundaryQuasiNeutrality) - nNodes = edge%nNodes - nodes = edge%getNodes(nNodes) - bound%alpha = 1.0d-1 + if (associated(edge%e1)) then + cell => edge%e1 - if (random() <= bound%alpha) then + else + cell => edge%e2 + + end if + + if (random() <= alpha) then call reflection(edge, part) else @@ -248,7 +346,7 @@ MODULE moduleMeshBoundary end subroutine quasiNeutrality !Points the boundary function to specific type - SUBROUTINE pointBoundaryFunction(edge, s) + module SUBROUTINE pointBoundaryFunction(edge, s) USE moduleErrors IMPLICIT NONE @@ -284,4 +382,4 @@ MODULE moduleMeshBoundary END SUBROUTINE pointBoundaryFunction -END MODULE moduleMeshBoundary +end submodule boundary diff --git a/src/modules/mesh/moduleMesh@common.f90 b/src/modules/mesh/moduleMesh@common.f90 new file mode 100644 index 0000000..fe710a8 --- /dev/null +++ b/src/modules/mesh/moduleMesh@common.f90 @@ -0,0 +1,3 @@ +module moduleMeshCommon + +end module moduleMeshCommon diff --git a/src/modules/mesh/moduleMesh@elements.f90 b/src/modules/mesh/moduleMesh@elements.f90 new file mode 100644 index 0000000..5b11ca1 --- /dev/null +++ b/src/modules/mesh/moduleMesh@elements.f90 @@ -0,0 +1,795 @@ +submodule(moduleMesh) elements + CONTAINS + !Reset the output of node + PURE module SUBROUTINE resetOutput(self) + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(meshNode), INTENT(inout):: self + INTEGER:: k + + DO k = 1, nSpecies + self%output(k)%den = 0.D0 + self%output(k)%mom = 0.D0 + self%output(k)%tensorS = 0.D0 + + END DO + + END SUBROUTINE resetOutput + + !Constructs the global K matrix + PURE module SUBROUTINE constructGlobalK(self) + IMPLICIT NONE + + CLASS(meshParticles), INTENT(inout):: self + INTEGER:: c + INTEGER, ALLOCATABLE:: n(:) + REAL(8), ALLOCATABLE:: localK(:,:) + INTEGER:: i, j + + DO c = 1, self%numCells + associate(nNodes => self%cells(c)%obj%nNodes) + ALLOCATE(n(1:nNodes)) + ALLOCATE(localK(1:nNodes, 1:nNodes)) + n = self%cells(c)%obj%getNodes(nNodes) + localK = self%cells(c)%obj%elemK(nNodes) + + DO i = 1, nNodes + DO j = 1, nNodes + self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j) + + END DO + + END DO + + DEALLOCATE(n, localK) + + end associate + + END DO + + END SUBROUTINE constructGlobalK + + ! Gather the value of valNodes at position Xi of an edge + pure module function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f) + implicit none + + class(meshEdge), intent(in):: self + real(8), intent(in):: Xi(1:3) + integer, intent(in):: nNodes + real(8), intent(in):: valNodes(1:nNodes) + real(8):: f + real(8):: fPsi(1:nNodes) + + fPsi = self%fPsi(Xi, nNodes) + f = dot_product(fPsi, valNodes) + + end function gatherF_edge_scalar + + !Gather the value of valNodes (scalar) at position Xi + PURE module FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes) + REAL(8):: f + REAL(8):: fPsi(1:nNodes) + + fPsi = self%fPsi(Xi, nNodes) + f = DOT_PRODUCT(fPsi, valNodes) + + END FUNCTION gatherF_cell_scalar + + !Gather the value of valNodes (array) at position Xi + PURE module FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes, 1:3) + REAL(8):: f(1:3) + REAL(8):: fPsi(1:nNodes) + + fPsi = self%fPsi(Xi, nNodes) + f = MATMUL(fPsi, valNodes) + + END FUNCTION gatherF_cell_array + + !Gather the spatial derivative of valNodes (scalar) at position Xi + PURE module FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df) + IMPLICIT NONE + + CLASS(meshCell), INTENT(in):: self + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes) + REAL(8):: df(1:3) + REAL(8):: dPsi(1:3, 1:nNodes) + REAL(8):: pDer(1:3,1:3) + REAL(8):: dPsiR(1:3, 1:nNodes) + REAL(8):: invJ(1:3, 1:3), detJ + + dPsi = self%dPsi(Xi, nNodes) + pDer = self%partialDer(nNodes, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + dPsiR = MATMUL(invJ, dPsi)/detJ + df = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), & + DOT_PRODUCT(dPsiR(2,:), valNodes), & + DOT_PRODUCT(dPsiR(3,:), valNodes) /) + + END FUNCTION gatherDF_cell_scalar + + !Scatters particle properties into cell nodes + module SUBROUTINE scatter(self, nNodes, part) + USE moduleMath + USE moduleSpecies + USE OMP_LIB + IMPLICIT NONE + + CLASS(meshCell), INTENT(inout):: self + INTEGER, INTENT(in):: nNodes + CLASS(particle), INTENT(in):: part + REAL(8):: fPsi(1:nNodes) + INTEGER:: cellNodes(1:nNodes) + REAL(8):: tensorS(1:3, 1:3) + INTEGER:: sp + INTEGER:: i + CLASS(meshNode), POINTER:: node + REAL(8):: pFraction !Particle fraction + + cellNodes = self%getNodes(nNodes) + fPsi = self%fPsi(part%Xi, nNodes) + + tensorS = outerProduct(part%v, part%v) + + sp = part%species%n + + DO i = 1, nNodes + node => mesh%nodes(cellNodes(i))%obj + pFraction = fPsi(i)*part%weight + CALL OMP_SET_LOCK(node%lock) + node%output(sp)%den = node%output(sp)%den + pFraction + node%output(sp)%mom(:) = node%output(sp)%mom(:) + pFraction*part%v(:) + node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS + CALL OMP_UNSET_LOCK(node%lock) + + END DO + + END SUBROUTINE scatter + + !Find next cell for particle + RECURSIVE module SUBROUTINE findCell(self, part, oldCell) + USE moduleSpecies + USE moduleErrors + USE OMP_LIB + IMPLICIT NONE + + CLASS(meshCell), INTENT(inout):: self + CLASS(particle), INTENT(inout), TARGET:: part + CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell + REAL(8):: Xi(1:3) + CLASS(meshElement), POINTER:: neighbourElement + INTEGER:: sp + + Xi = self%phy2log(part%r) + !Checks if particle is inside 'self' cell + IF (self%inside(Xi)) THEN + part%cell = self%n + part%Xi = Xi + part%n_in = .TRUE. + !Assign particle to listPart_in + IF (listInCells) THEN + CALL OMP_SET_LOCK(self%lock) + sp = part%species%n + CALL self%listPart_in(sp)%add(part) + self%totalWeight(sp) = self%totalWeight(sp) + part%weight + CALL OMP_UNSET_LOCK(self%lock) + + END IF + + ELSE + !If not, searches for a neighbour and repeats the process. + CALL self%neighbourElement(Xi, neighbourElement) + !Defines the next step + SELECT TYPE(neighbourElement) + CLASS IS(meshCell) + !Particle moved to new cell, repeat find procedure + CALL neighbourElement%findCell(part, self) + + CLASS IS (meshEdge) + !Particle encountered a surface, apply boundary + CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part) + + !If particle is still inside the domain, call findCell + IF (part%n_in) THEN + IF(PRESENT(oldCell)) THEN + CALL self%findCell(part, oldCell) + + ELSE + CALL self%findCell(part) + + END IF + END IF + + CLASS DEFAULT + WRITE (*, "(A, I6)") "Element = ", self%n + CALL criticalError("No connectivity found for element", "findCell") + + END SELECT + + END IF + + END SUBROUTINE findCell + + !If Coll and Particle are the same, simply copy the part%cell into part%cellColl + module SUBROUTINE findCellSameMesh(part) + USE moduleSpecies + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + + part%cellColl = part%cell + + END SUBROUTINE findCellSameMesh + + !TODO: try to combine this with the findCell for a regular mesh + !Find the volume in which particle reside in the mesh for collisions + !No boundary interaction taken into account + module SUBROUTINE findCellCollMesh(part) + USE moduleSpecies + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + LOGICAL:: found + CLASS(meshCell), POINTER:: cell + REAL(8), DIMENSION(1:3):: Xi + CLASS(meshElement), POINTER:: neighbourElement + INTEGER:: sp + + found = .FALSE. + + cell => meshColl%cells(part%cellColl)%obj + DO WHILE(.NOT. found) + Xi = cell%phy2log(part%r) + IF (cell%inside(Xi)) THEN + part%cellColl = cell%n + IF (listInCells) THEN + CALL OMP_SET_LOCK(cell%lock) + sp = part%species%n + CALL cell%listPart_in(sp)%add(part) + cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight + CALL OMP_UNSET_LOCK(cell%lock) + + END IF + found = .TRUE. + + ELSE + CALL cell%neighbourElement(Xi, neighbourElement) + SELECT TYPE(neighbourElement) + CLASS IS(meshCell) + !Try next element + cell => neighbourElement + + CLASS DEFAULT + !Should never happend, but just in case, stops loops + found = .TRUE. + + END SELECT + + END IF + + END DO + + END SUBROUTINE findCellCollMesh + + !Returns index of volume associated to a position (if any) + !If no voulme is found, returns 0 + !WARNING: This function is slow and should only be used in initialization phase + module FUNCTION findCellBrute(self, r) RESULT(nVol) + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshGeneric), INTENT(in):: self + REAL(8), DIMENSION(1:3), INTENT(in):: r + INTEGER:: nVol + INTEGER:: e + REAL(8), DIMENSION(1:3):: Xi + + !Inits RESULT + nVol = 0 + + DO e = 1, self%numCells + Xi = self%cells(e)%obj%phy2log(r) + IF(self%cells(e)%obj%inside(Xi)) THEN + nVol = self%cells(e)%obj%n + EXIT + + END IF + + END DO + + END FUNCTION findCellBrute + + !Computes collisions in element + module SUBROUTINE doCollisions(self) + USE moduleCollisions + USE moduleSpecies + USE moduleList + use moduleRefParam + USE moduleRandom + USE moduleOutput + USE moduleMath + USE moduleCaseParam, ONLY: timeStep + IMPLICIT NONE + + CLASS(meshGeneric), INTENT(inout), TARGET:: self + INTEGER:: e + CLASS(meshCell), POINTER:: cell + INTEGER:: k, i, j + INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell + REAL(8):: pMax !Maximum probability of collision + INTEGER:: nColl + TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:) + TYPE(particle), POINTER:: part_i, part_j + INTEGER:: n, c + REAL(8):: vRel, rMass, eRel + REAL(8):: sigmaVrelTotal + REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:) + REAL(8):: rnd_real !Random number for collision + INTEGER:: rnd_int !Random number for collision + + IF (MOD(timeStep, everyColl) == 0) THEN + !Collisions need to be performed in this iteration + !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) + DO e=1, self%numCells + + cell => self%cells(e)%obj + + !TODO: Simplify this, to many sublevels + !Iterate over the number of pairs + DO k = 1, nCollPairs + !Reset tally of collisions + IF (collOutput) THEN + cell%tallyColl(k)%tally = 0 + + END IF + + IF (interactionMatrix(k)%amount > 0) THEN + !Select the species for the collision pair + i = interactionMatrix(k)%sp_i%n + j = interactionMatrix(k)%sp_j%n + + !Number of particles per species in the collision pair + nPart_i = cell%listPart_in(i)%amount + nPart_j = cell%listPart_in(j)%amount + + IF (nPart_i > 0 .AND. nPart_j > 0) THEN + !Total number of particles for the collision pair + nPart = nPart_i + nPart_j + + !Resets the number of collisions in the cell + nColl = 0 + + !Probability of collision for pair i-j + pMax = (cell%totalWeight(i) + cell%totalWeight(j))*cell%sigmaVrelMax(k)*tauColl/cell%volume + + !Number of collisions in the cell + nColl = NINT(REAL(nPart)*pMax*0.5D0) + + !Converts the list of particles to an array for easy access + IF (nColl > 0) THEN + partTemp_i = cell%listPart_in(i)%convert2Array() + partTemp_j = cell%listPart_in(j)%convert2Array() + + END IF + + DO n = 1, nColl + !Select random particles + part_i => NULL() + part_j => NULL() + rnd_int = random(1, nPart_i) + part_i => partTemp_i(rnd_int)%part + rnd_int = random(1, nPart_j) + part_j => partTemp_j(rnd_int)%part + !If they are the same particle, skip + !TODO: Maybe try to improve this + IF (ASSOCIATED(part_i, part_j)) THEN + CYCLE + + END IF + + !If particles do not belong to the species, skip collision + !This can happen, for example, if particle has been previously ionized or removed + !TODO: Try to find a way to not lose these collisions. Maybe check new 'k' and use that for the collision? + IF (part_i%species%n /= i .OR. & + part_j%species%n /= j) THEN + CYCLE + + END IF + !Obtain the cross sections for the different processes + !TODO: From here it might be a procedure in interactionMatrix + vRel = NORM2(part_i%v-part_j%v) + rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m) + eRel = rMass*vRel**2 + CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel) + + !Update maximum sigma*v_rel + IF (sigmaVrelTotal > cell%sigmaVrelMax(k)) THEN + cell%sigmaVrelMax(k) = sigmaVrelTotal + + END IF + + ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount)) + probabilityColl = 0.0 + DO c = 1, interactionMatrix(k)%amount + probabilityColl(c) = sigmaVrel(c)/cell%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1)) + + END DO + + !Selects random number between 0 and 1 + rnd_real = random() + + !If the random number is below the total probability of collision, collide particles + IF (rnd_real < sigmaVrelTotal / cell%sigmaVrelMax(k)) THEN + + !Loop over collisions + DO c = 1, interactionMatrix(k)%amount + IF (rnd_real <= probabilityColl(c)) THEN + !$OMP CRITICAL + CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel) + !$OMP END CRITICAL + + !If collisions are gonna be output, count the collision + IF (collOutput) THEN + cell%tallyColl(k)%tally(c) = cell%tallyColl(k)%tally(c) + 1 + + END IF + + !A collision has ocurred, exit the loop + EXIT + + END IF + + END DO + + END IF + + !Deallocate arrays for next collision + DEALLOCATE(sigmaVrel, probabilityColl) + + !End loop collisions in cell + END DO + + END IF + + END IF + + !End loop collision pairs + END DO + + !End loop volumes + END DO + !$OMP END DO + + END IF + + END SUBROUTINE doCollisions + + module SUBROUTINE doCoulomb(self) + USE moduleCoulomb + USE moduleRandom + USE moduleOutput + USE moduleList + USE moduleMath + USE moduleRefParam + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshParticles), INTENT(in), TARGET:: self + CLASS(meshCell), POINTER:: cell + TYPE(interactionsCoulomb):: pair + INTEGER:: e + INTEGER:: k + INTEGER:: i, j + INTEGER:: n + INTEGER:: p + TYPE(lNode), POINTER:: partTemp + INTEGER(8), ALLOCATABLE:: cellNodes(:) + CLASS(meshNode), POINTER:: node + TYPE(outputFormat):: output + REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node + REAL(8):: density, velocity(1:3), temperature!values at particle position + REAL(8):: C(1:3), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference + REAL(8):: l, lW, l2 + REAL(8):: GlW, HlW + REAL(8):: normC + REAL(8):: cosThe, sinThe + REAL(8):: cosPhi, sinPhi + REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame + REAL(8):: A, AW + REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W + REAL(8):: theta_per !Random angle for perpendicular direction + REAL(8):: eps = 1.D-12 + REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ij, p_ij + REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ij + REAL(8):: massSum_ij + REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ji, p_ji + REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ji + REAL(8):: massSum_ji + REAL(8):: alpha_num, alpha_den, alpha, beta(1:3) + + + !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) + DO e = 1, self%numCells + cell => self%cells(e)%obj + cellNodes = cell%getNodes(cell%nNodes) + + ALLOCATE(densityNodes(1:cell%nNodes), & + velocityNodes(1:cell%nNodes, 1:3), & + temperatureNodes(1:cell%nNodes)) + + DO k=1, nCoulombPairs + pair = coulombMatrix(k) + i = pair%sp_i%n + j = pair%sp_j%n + + !Do scattering of particles from species_i due to species j + !Compute background properties of species_j + DO n = 1, cell%nNodes + node => self%nodes(cellNodes(n))%obj + CALL calculateOutput(node%output(j), output, node%v, pair%sp_j) + densityNodes(n) = output%density/n_ref + velocityNodes(n,1:3) = output%velocity(1:3)/v_ref + temperatureNodes(n) = output%temperature/T_ref + + END DO + + ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3)) + ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3)) + ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount)) + deltaV_ij = 0.D0 + p_ij = 0.D0 + mass_ij = 0.D0 + !Loop over particles of species_i + partTemp => cell%listPart_in(i)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) + velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) + temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) + + !If cell temperature is too low, skip particle to avoid division by zero + IF (temperature>eps) THEN + l2 = pair%l2_j/temperature + l = SQRT(l2) + + ELSE + partTemp => partTemp%next + + CYCLE + + END IF + + A = pair%A_i*density + + C = partTemp%part%v - velocity + normC = NORM2(C) + + !C_3 = z; C_1, C2 = x, y (per) + C_per = NORM2(C(1:2)) + cosPhi = C(1) / C_per + sinPhi = C(2) / C_per + cosThe = C(3) / normC + sinThe = C_per / normC + + !Rotation matrix to go from W to C + rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column + -sinPhi, cosPhi, 0.D0, & !Second column + sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column + (/ 3, 3 /)) + + !W at start is = (0, 0, normC), so normW = normC + lW = l * normC + GlW = G(lW) + HlW = H(lW) + AW = A / normC + + !Calculate changes in W due to collision process + deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin + deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() + deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() + + !Random angle to distribute perpendicular change in velocity + theta_per = PI2*random() + + !Change W + W(1) = deltaW_per_square * COS(theta_per) + W(2) = deltaW_per_square * SIN(theta_per) + W(3) = normC + deltaW_par + deltaW_par_square + + !Compute changes in velocity for each particle + deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v + mass_ij(p) = pair%sp_i%m*partTemp%part%weight + p_ij(p,1:3) = mass_ij(p)*partTemp%part%v + + !Move to the next particle in the list + partTemp => partTemp%next + p = p + 1 + + END DO + + !Do corresponding collisions + IF (i /= j) THEN + !Do scattering of particles from species_j due to species i + !Compute background properties of species_i + DO n = 1, cell%nNodes + node => self%nodes(cellNodes(n))%obj + CALL calculateOutput(node%output(i), output, node%v, pair%sp_i) + densityNodes(n) = output%density/n_ref + velocityNodes(n,1:3) = output%velocity(1:3)/v_ref + temperatureNodes(n) = output%temperature/T_ref + + END DO + + ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3)) + ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3)) + ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount)) + deltaV_ji = 0.D0 + p_ji = 0.D0 + mass_ji = 0.D0 + !Loop over particles of species_j + partTemp => cell%listPart_in(j)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) + velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) + temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) + + !If cell temperature is too low, skip particle to avoid division by zero + IF (temperature>eps) THEN + l2 = pair%l2_i/temperature + l = SQRT(l2) + + ELSE + partTemp => partTemp%next + + CYCLE + + END IF + A = pair%A_j*density + + C = partTemp%part%v - velocity + normC = NORM2(C) + + !C_3 = z; C_1, C2 = x, y (per) + C_per = NORM2(C(1:2)) + cosPhi = C(1) / C_per + sinPhi = C(2) / C_per + cosThe = C(3) / normC + sinThe = C_per / normC + + !Rotation matrix to go from W to C + rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column + -sinPhi, cosPhi, 0.D0, & !Second column + sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column + (/ 3, 3 /)) + + !W at start is = (0, 0, normC), so normW = normC + lW = l * normC + GlW = G(lW) + HlW = H(lW) + AW = A / normC + + !Calculate changes in W due to collision process + deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin + deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() + deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() + + !Random angle to distribute perpendicular change in velocity + theta_per = PI2*random() + + !Change W + W(1) = deltaW_per_square * COS(theta_per) + W(2) = deltaW_per_square * SIN(theta_per) + W(3) = normC + deltaW_par + deltaW_par_square + + !Compute changes in velocity for each particle + deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v + mass_ji(p) = pair%sp_j%m*partTemp%part%weight + p_ji(p,1:3) = mass_ji(p)*partTemp%part%v + + !Move to the next particle in the list + partTemp => partTemp%next + p = p + 1 + + END DO + + END IF + + !Calculate correction + !Total mass + massSum_ij = SUM(mass_ij) + massSum_ji = 0.D0 + + !Beta + beta = 0.D0 + DO p = 1, cell%listPart_in(i)%amount + beta = beta + mass_ij(p) * deltaV_ij(p,1:3) + + END DO + + IF (i /= j) THEN + massSum_ji = SUM(mass_ji) + DO p = 1, cell%listPart_in(j)%amount + beta = beta + mass_ji(p) * deltaV_ji(p,1:3) + + END DO + + END IF + + beta = beta / (massSum_ij + massSum_ji) + + !Alpha + alpha_num = 0.D0 + alpha_den = 0.D0 + DO p =1, cell%listPart_in(i)%amount + alpha_num = alpha_num + DOT_PRODUCT(p_ij(p,1:3), deltav_ij(p,1:3) - beta(1:3)) + alpha_den = alpha_den + mass_ij(p) * NORM2(deltav_ij(p,1:3) - beta(1:3))**2 + + END DO + + IF (i /= j) THEN + DO p = 1, cell%listPart_in(j)%amount + alpha_num = alpha_num + DOT_PRODUCT(p_ji(p,1:3), deltav_ji(p,1:3) - beta(1:3)) + alpha_den = alpha_den + mass_ji(p) * NORM2(deltav_ji(p,1:3) - beta(1:3))**2 + + END DO + + END IF + + alpha = -2.D0*alpha_num / alpha_den + + !Apply correction to particles velocity + partTemp => cell%listPart_in(i)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + partTemp%part%v = partTemp%part%v + alpha * (deltaV_ij(p,1:3) - beta(1:3)) + partTemp => partTemp%next + p = p + 1 + + END DO + + IF (i /= j) THEN + partTemp => cell%listPart_in(j)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + partTemp%part%v = partTemp%part%v + alpha * (deltaV_ji(p,1:3) - beta(1:3)) + partTemp => partTemp%next + p = p + 1 + + END DO + + END IF + + DEALLOCATE(deltaV_ij, p_ij, mass_ij) + + IF (i /= j) THEN + DEALLOCATE(deltaV_ji, p_ji, mass_ji) + + END IF + + END DO + + DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes) + + END DO + !$OMP END DO + + END SUBROUTINE doCoulomb + +end submodule elements diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 deleted file mode 100644 index 79d7b3d..0000000 --- a/src/modules/moduleBoundary.f90 +++ /dev/null @@ -1,178 +0,0 @@ -MODULE moduleBoundary - USE moduleTable - USE moduleSpecies - - !Generic type for boundaries - TYPE, PUBLIC:: boundaryGeneric - CONTAINS - - END TYPE boundaryGeneric - - !Reflecting boundary - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryReflection - CONTAINS - - END TYPE boundaryReflection - - !Absorption boundary - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAbsorption - CONTAINS - - END TYPE boundaryAbsorption - - !Transparent boundary - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryTransparent - CONTAINS - - END TYPE boundaryTransparent - - !Symmetry axis - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis - CONTAINS - - END TYPE boundaryAxis - - !Wall Temperature boundary - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature - !Thermal velocity of the wall: square root(Wall temperature X specific heat) - REAL(8):: vTh - CONTAINS - - END TYPE boundaryWallTemperature - - !Ionization boundary - TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization - REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals. - CLASS(speciesGeneric), POINTER:: species !Ion species - CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron - TYPE(table1D):: crossSection - REAL(8):: effectiveTime - REAL(8):: eThreshold - REAL(8):: deltaV - CONTAINS - - END TYPE boundaryIonization - - ! Ensures quasi-neutrality by changing the reflection coefficient - type, public, extends(boundaryGeneric):: boundaryQuasiNeutrality - real(8):: alpha ! Reflection parameter - - end type boundaryQuasiNeutrality - !Wrapper for boundary types (one per species) - TYPE:: bTypesCont - CLASS(boundaryGeneric), ALLOCATABLE:: obj - CONTAINS - - END TYPE bTypesCont - - !Wrapper for boundary conditions - TYPE:: boundaryCont - INTEGER:: n = 0 - CHARACTER(:), ALLOCATABLE:: name - INTEGER:: physicalSurface = 0 !Physical surface as defined in the mesh file - CLASS(bTypesCont), ALLOCATABLE:: bTypes(:) !Array for boundary per species - CONTAINS - - END TYPE boundaryCont - - !Number of boundaries - INTEGER:: nBoundary = 0 - !Array for boundaries - TYPE(boundaryCont), ALLOCATABLE, TARGET:: boundary(:) - - CONTAINS - FUNCTION getBoundaryId(physicalSurface) RESULT(id) - IMPLICIT NONE - - INTEGER:: physicalSurface - INTEGER:: id - INTEGER:: i - - id = 0 - DO i = 1, nBoundary - IF (physicalSurface == boundary(i)%physicalSurface) id = boundary(i)%n - - END DO - - END FUNCTION getBoundaryId - - SUBROUTINE initWallTemperature(boundary, T, c) - USE moduleRefParam - IMPLICIT NONE - - CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary - REAL(8), INTENT(in):: T, c !Wall temperature and specific heat - REAL(8):: vTh - - vTh = DSQRT(c * T) / v_ref - boundary = boundaryWallTemperature(vTh = vTh) - - END SUBROUTINE initWallTemperature - - SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary) - USE moduleRefParam - USE moduleSpecies - USE moduleCaseParam - USE moduleConstParam - USE moduleErrors - IMPLICIT NONE - - CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary - REAL(8), INTENT(in):: me !Electron mass - REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties - INTEGER, INTENT(in):: ion - INTEGER, OPTIONAL, INTENT(in):: electronSecondary - REAL(8):: effTime - CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection - REAL(8), INTENT(in):: eThreshold - - ALLOCATE(boundaryIonization:: boundary) - - SELECT TYPE(boundary) - TYPE IS(boundaryIonization) - boundary%m0 = m0 / m_ref - boundary%n0 = n0 * Vol_ref - boundary%v0 = v0 / v_ref - boundary%vTh = DSQRT(kb*T0/m0)/v_ref - boundary%species => species(ion)%obj - IF (PRESENT(electronSecondary)) THEN - SELECT TYPE(sp => species(electronSecondary)%obj) - TYPE IS(speciesCharged) - boundary%electronSecondary => sp - - CLASS DEFAULT - CALL criticalError("Species " // sp%name // " chosen for " // & - "secondary electron is not a charged species", 'initIonization') - - END SELECT - - ELSE - boundary%electronSecondary => NULL() - - END IF - boundary%effectiveTime = effTime / ti_ref - CALL boundary%crossSection%init(crossSection) - CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) - boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2) - boundary%deltaV = DSQRT(boundary%eThreshold/me) - - END SELECT - - END SUBROUTINE initIonization - - subroutine initQuasiNeutrality(boundary) - implicit none - - class(boundaryGeneric), allocatable, intent(out):: boundary - - allocate(boundaryQuasiNeutrality:: boundary) - - select type(boundary) - type is(boundaryQuasiNeutrality) - boundary%alpha = 0.d0 - - end select - - end subroutine initQuasiNeutrality - -END MODULE moduleBoundary diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 3962986..d9a0d4e 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -324,7 +324,6 @@ MODULE moduleInject USE moduleMesh USE moduleRandom USE moduleErrors - use moduleMeshBoundary, only: reflection IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self