JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Implementation of different distribution functions for velocities.

Browse source 
Maxwellian and Diract Delta distributions have been implemented.

The input for injection of particles should be rewritten to allow more
clear input file.
This commit is contained in:
Jorge Gonzalez 2020-12-13 13:56:48 +01:00
commit 37b0139b1f
37 changed files with 252 additions and 5497 deletions
Download patch file Download diff file Expand all files Collapse all files

103
src/fpakc.f95
View file

@ -1,103 +0,0 @@
! FPAKC main program
PROGRAM fpakc
USE moduleInput
USE moduleErrors
USE OMP_LIB
USE moduleInject
USE moduleSolver
USE moduleOutput
USE moduleCompTime
USE moduleMesh
IMPLICIT NONE
! t = time step
INTEGER:: t = 0
! arg1 = Input argument 1 (input file)
CHARACTER(200):: arg1
! inputFile = path+name of input file
CHARACTER(:), ALLOCATABLE:: inputFile
tStep = omp_get_wtime()
!Gets the input file
CALL getarg(1, arg1)
inputFile = TRIM(arg1)
!If no input file is declared, a critical error is called
IF (inputFile == "") CALL criticalError("No input file provided", "fpakc")
!Reads the json configuration file
CALL readConfig(inputFile)
CALL verboseError("Calculating initial EM field...")
CALL doEMField()
tStep = omp_get_wtime() - tStep
!Output initial state
CALL doOutput(t)
CALL verboseError('Starting main loop...')
!$OMP PARALLEL DEFAULT(SHARED)
DO t = 1, tmax
!Insert new particles and push them
!$OMP SINGLE
tStep = omp_get_wtime()
!Checks if a species needs to me moved in this iteration
CALL solver%updatePushSpecies(t)
tPush = omp_get_wtime()
!$OMP END SINGLE
!Injects new particles
CALL doInjects()
!Push old particles
CALL doPushes()
!$OMP SINGLE
tPush = omp_get_wtime() - tPush
!Collisions
tColl = omp_get_wtime()
!$OMP END SINGLE
CALL doCollisions()
!$OMP SINGLE
tColl = omp_get_wtime() - tColl
!Reset particles
tReset = omp_get_wtime()
!$OMP END SINGLE
CALL doReset()
!$OMP SINGLE
tReset = omp_get_wtime() - tReset
!Weight
tWeight = omp_get_wtime()
!$OMP END SINGLE
CALL doScatter()
!$OMP SINGLE
tWeight = omp_get_wtime() - tWeight
tEMField = omp_get_wtime()
!$OMP END SINGLE
CALL doEMField()
!$OMP SINGLE
tEMField = omp_get_wtime() - tEMField
tStep = omp_get_wtime() - tStep
!Output data
CALL doOutput(t)
!$OMP END SINGLE
END DO
!$OMP END PARALLEL
END PROGRAM fpakc

4
src/makefile
View file

@ -9,8 +9,8 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolv
all: $(OUTPUT)
$(OUTPUT): modules.o $(OUTPUT).f95
$(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f95
$(OUTPUT): modules.o $(OUTPUT).f90
$(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f90
$(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(OBJDIR)/$(OUTPUT).o $(JSONLIB) -L/usr/local/lib -llapack -lopenblas
modules.o:

42
src/modules/makefile
View file

@ -9,35 +9,35 @@ all: $(OBJS) mesh.o
mesh.o: moduleCollisions.o moduleBoundary.o
$(MAKE) -C mesh all
moduleCollisions.o: moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleCollisions.o: moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleInput.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleInput.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleInject.o: moduleSpecies.o moduleSolver.o moduleInject.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleInject.o: moduleSpecies.o moduleSolver.o moduleInject.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleRefParam.o: moduleConstParam.o moduleRefParam.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleRefParam.o: moduleConstParam.o moduleRefParam.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleSpecies.o: moduleErrors.o moduleCaseParam.o moduleSpecies.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleSpecies.o: moduleErrors.o moduleCaseParam.o moduleSpecies.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleSolver.o: mesh.o moduleList.o moduleEM.o moduleSpecies.o moduleRefParam.o moduleOutput.o moduleSolver.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleSolver.o: mesh.o moduleList.o moduleEM.o moduleSpecies.o moduleRefParam.o moduleOutput.o moduleSolver.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleTable.o: moduleErrors.o moduleTable.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleTable.o: moduleErrors.o moduleTable.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleEM.o: mesh.o moduleSpecies.o moduleEM.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleEM.o: mesh.o moduleSpecies.o moduleEM.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
%.o: %.f95
%.o: %.f90
$(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@

12
src/modules/mesh/1DCart/makefile
View file

@ -1,11 +1,11 @@
all: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.o
moduleMesh1D.o: moduleMesh1D.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleMesh1D.o: moduleMesh1D.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh1DBoundary.o: moduleMesh1D.o moduleMesh1DBoundary.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleMesh1DBoundary.o: moduleMesh1D.o moduleMesh1DBoundary.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh1DRead.o: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleMesh1DRead.o: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

489
src/modules/mesh/1DCart/moduleMesh1D.f95
View file

@ -1,489 +0,0 @@
!moduleMesh1D: 1D cartesian module
! x == x
! y == unused
! z == unused
MODULE moduleMesh1D
USE moduleMesh
IMPLICIT NONE
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1D
!Element coordinates
REAL(8):: x = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initNode1D
PROCEDURE, PASS:: getCoordinates => getCoord1D
END TYPE meshNode1D
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdge1D
!Element coordinates
REAL(8):: x = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL()
CONTAINS
PROCEDURE, PASS:: init => initEdge1D
PROCEDURE, PASS:: getNodes => getNodes1D
PROCEDURE, PASS:: randPos => randPos1D
END TYPE meshEdge1D
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1D
CONTAINS
PROCEDURE, PASS:: detJac => detJ1D
PROCEDURE, PASS:: invJac => invJ1D
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol1D
ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xi) RESULT(fPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
END FUNCTION dPsi_interface
PURE SUBROUTINE partialDer_interface(self, dPsi, dx)
IMPORT meshVol1D
CLASS(meshVol1D), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1):: dx
END SUBROUTINE partialDer_interface
END INTERFACE
TYPE, PUBLIC, EXTENDS(meshVol1D):: meshVol1DSegm
!Element coordinates
REAL(8):: x(1:2)
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
!Connectivity to adjacent elements
CLASS(*), POINTER:: e1 => NULL(), e2 => NULL()
REAL(8):: arNodes(1:2)
CONTAINS
PROCEDURE, PASS:: init => initVol1DSegm
PROCEDURE, PASS:: area => areaSegm
PROCEDURE, NOPASS:: fPsi => fPsiSegm
PROCEDURE, NOPASS:: dPsi => dPsiSegm
PROCEDURE, PASS:: partialDer => partialDerSegm
PROCEDURE, PASS:: elemK => elemKSegm
PROCEDURE, PASS:: elemF => elemFSegm
PROCEDURE, NOPASS:: weight => weightSegm
PROCEDURE, NOPASS:: inside => insideSegm
PROCEDURE, PASS:: scatter => scatterSegm
PROCEDURE, PASS:: gatherEF => gatherEFSegm
PROCEDURE, PASS:: getNodes => getNodesSegm
PROCEDURE, PASS:: phy2log => phy2logSegm
PROCEDURE, PASS:: nextElement => nextElementSegm
PROCEDURE, PASS:: resetOutput => resetOutputSegm
END TYPE meshVol1DSegm
CONTAINS
!NODE FUNCTIONS
!Init node element
SUBROUTINE initNode1D(self, n, r)
USE moduleSpecies
USE moduleRefParam
IMPLICIT NONE
CLASS(meshNode1D), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
self%n = n
self%x = r(1)/L_ref
!Node volume, to be determined in mesh
self%v = 0.D0
!Allocates output
ALLOCATE(self%output(1:nSpecies))
END SUBROUTINE initNode1D
PURE FUNCTION getCoord1D(self) RESULT(r)
IMPLICIT NONE
CLASS(meshNode1D), INTENT(in):: self
REAL(8):: r(1:3)
r = (/ self%x, 0.D0, 0.D0 /)
END FUNCTION getCoord1D
!EDGE FUNCTIONS
!Inits edge element
SUBROUTINE initEdge1D(self, n, p, bt, physicalSurface)
IMPLICIT NONE
CLASS(meshEdge1D), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
REAL(8), DIMENSION(1:3):: r1
self%n = n
self%n1 => mesh%nodes(p(1))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
self%x = r1(1)
self%normal = (/ 1.D0, 0.D0, 0.D0 /)
!Boundary index
self%bt = bt
!Physical Surface
self%physicalSurface = physicalSurface
END SUBROUTINE initEdge1D
!Get nodes from edge
PURE FUNCTION getNodes1D(self) RESULT(n)
IMPLICIT NONE
CLASS(meshEdge1D), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1))
n = (/ self%n1%n /)
END FUNCTION getNodes1D
!Calculates a 'random' position in edge
FUNCTION randPos1D(self) RESULT(r)
CLASS(meshEdge1D), INTENT(in):: self
REAL(8):: r(1:3)
r = (/ self%x, 0.D0, 0.D0 /)
END FUNCTION randPos1D
!VOLUME FUNCTIONS
!SEGMENT FUNCTIONS
!Init segment element
SUBROUTINE initVol1DSegm(self, n, p)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
REAL(8), DIMENSION(1:3):: r1, r2
self%n = n
self%n1 => mesh%nodes(p(1))%obj
self%n2 => mesh%nodes(p(2))%obj
!Get element coordinates
r1 = self%n1%getCoordinates()
r2 = self%n2%getCoordinates()
self%x = (/ r1(1), r2(1) /)
!Assign node volume
CALL self%area()
self%n1%v = self%n1%v + self%arNodes(1)
self%n2%v = self%n2%v + self%arNodes(2)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
END SUBROUTINE initVol1DSegm
!Computes element area
PURE SUBROUTINE areaSegm(self)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout):: self
REAL(8):: l !element length
REAL(8):: fPsi(1:2)
REAL(8):: detJ
REAL(8):: Xii(1:3)
self%volume = 0.D0
self%arNodes = 0.D0
!1 point Gauss integral
Xii = 0.D0
fPsi = self%fPsi(Xii)
detJ = self%detJac(Xii)
l = 2.D0*detJ
self%volume = l
self%arNodes = fPsi*l
END SUBROUTINE areaSegm
!Computes element functions at point xii
PURE FUNCTION fPsiSegm(xi) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:2))
fPsi(1) = 1.D0 - xi(1)
fPsi(2) = 1.D0 + xi(1)
fPsi = fPsi * 5.D-1
END FUNCTION fPsiSegm
!Computes element derivative shape function at Xii
PURE FUNCTION dPsiSegm(xi) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
ALLOCATE(dPsi(1:1, 1:2))
dPsi(1, 1) = -5.D-1
dPsi(1, 2) = 5.D-1
END FUNCTION dPsiSegm
!Computes partial derivatives of coordinates
PURE SUBROUTINE partialDerSegm(self, dPsi, dx)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1):: dx
dx(1) = DOT_PRODUCT(dPsi(1,:), self%x)
END SUBROUTINE partialDerSegm
!Computes local stiffness matrix
PURE FUNCTION elemKSegm(self) RESULT(ke)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8):: ke(1:2,1:2)
REAL(8):: Xii(1:3)
REAL(8):: dPsi(1:1, 1:2)
REAL(8):: invJ
ke = 0.D0
Xii = (/ 0.D0, 0.D0, 0.D0 /)
dPsi = self%dPsi(Xii)
invJ = self%invJac(Xii, dPsi)
ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /)
ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /)
ke = 2.D0*ke*invJ
END FUNCTION elemKSegm
PURE FUNCTION elemFSegm(self, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:2)
REAL(8):: detJ
REAL(8):: Xii(1:3)
Xii = 0.D0
fPsi = self%fPsi(Xii)
detJ = self%detJac(Xii)
ALLOCATE(localF(1:2))
localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ
END FUNCTION elemFSegm
PURE FUNCTION weightSegm(xi) RESULT(w)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: w(1:3)
w = fPsiSegm(xi)
END FUNCTION weightSegm
PURE FUNCTION insideSegm(xi) RESULT(ins)
IMPLICIT NONE
REAL(8), INTENT(in):: xi(1:3)
LOGICAL:: ins
ins = xi(1) >=-1.D0 .AND. &
xi(1) <= 1.D0
END FUNCTION insideSegm
SUBROUTINE scatterSegm(self, part)
USE moduleOutput
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
CLASS(particle), INTENT(in):: part
TYPE(outputNode), POINTER:: vertex
REAL(8):: w_p(1:2)
REAL(8):: tensorS(1:3,1:3)
w_p = self%weight(part%xi)
tensorS = outerProduct(part%v, part%v)
vertex => self%n1%output(part%sp)
vertex%den = vertex%den + part%weight*w_p(1)
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
vertex => self%n2%output(part%sp)
vertex%den = vertex%den + part%weight*w_p(2)
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:)
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
END SUBROUTINE scatterSegm
!Gathers EF at position Xii
PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1, 1:2)
REAL(8):: phi(1:2)
REAL(8):: EF(1:3)
REAL(8):: invJ
phi = (/ self%n1%emData%phi, &
self%n2%emData%phi /)
dPsi = self%dPsi(xi)
invJ = self%invJac(xi, dPsi)
EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ
EF(2) = 0.D0
EF(3) = 0.D0
END FUNCTION gatherEFSegm
!Get nodes from 1D volume
PURE FUNCTION getNodesSegm(self) RESULT(n)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1:2))
n = (/ self%n1%n, self%n2%n /)
END FUNCTION getNodesSegm
PURE FUNCTION phy2logSegm(self, r) RESULT(xN)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
xN = 0.D0
xN(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0
END FUNCTION phy2logSegm
!Get next element for a logical position xi
SUBROUTINE nextElementSegm(self, xi, nextElement)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), POINTER, INTENT(out):: nextElement
NULLIFY(nextElement)
IF (xi(1) < -1.D0) THEN
nextElement => self%e2
ELSEIF (xi(1) > 1.D0) THEN
nextElement => self%e1
END IF
END SUBROUTINE nextElementSegm
!Reset the output of nodes in element
PURE SUBROUTINE resetOutputSegm(self)
USE moduleSpecies
USE moduleOutput
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout):: self
INTEGER:: k
DO k = 1, nSpecies
self%n1%output(k)%den = 0.D0
self%n1%output(k)%mom = 0.D0
self%n1%output(k)%tensorS = 0.D0
self%n2%output(k)%den = 0.D0
self%n2%output(k)%mom = 0.D0
self%n2%output(k)%tensorS = 0.D0
END DO
END SUBROUTINE resetOutputSegm
!COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS
!Computes the element Jacobian determinant
PURE FUNCTION detJ1D(self, xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE
CLASS(meshVol1D), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dJ
REAL(8):: dx(1)
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
dJ = dx(1)
END FUNCTION detJ1D
!Computes the invers Jacobian
PURE FUNCTION invJ1D(self, xi, dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshVol1D), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1)
REAL(8):: invJ
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx)
invJ = 1.D0/dx(1)
END FUNCTION invJ1D
END MODULE moduleMesh1D

40
src/modules/mesh/1DCart/moduleMesh1DBoundary.f95
View file

@ -1,40 +0,0 @@
MODULE moduleMesh1DBoundary
USE moduleMesh1D
TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DRef
CONTAINS
PROCEDURE, PASS:: fBoundary => reflection
END TYPE meshEdge1DRef
TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DAbs
CONTAINS
PROCEDURE, PASS:: fBoundary => absorption
END TYPE meshEdge1DAbs
CONTAINS
SUBROUTINE reflection(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge1DRef), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
part%v(1) = -part%v(1)
part%r(1) = 2.D0*self%x - part%r(1)
END SUBROUTINE reflection
SUBROUTINE absorption(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge1DAbs), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
part%n_in = .FALSE.
END SUBROUTINE absorption
END MODULE moduleMesh1DBoundary

268
src/modules/mesh/1DCart/moduleMesh1DRead.f95
View file

@ -1,268 +0,0 @@
MODULE moduleMesh1DRead
USE moduleMesh
USE moduleMesh1D
USE moduleMesh1DBoundary
!TODO: make this abstract to allow different mesh formats
TYPE, EXTENDS(meshGeneric):: mesh1DGeneric
CONTAINS
PROCEDURE, PASS:: init => init1DMesh
PROCEDURE, PASS:: readMesh => readMesh1D
END TYPE
INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connected
CONTAINS
!Init 1D mesh
SUBROUTINE init1DMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(mesh1DGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
SELECT CASE(meshFormat)
CASE ("gmsh")
self%printOutput => printOutputGmsh
self%printColl => printCollGmsh
self%printEM => printEMGmsh
CASE DEFAULT
CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1D")
END SELECT
END SUBROUTINE init1DMesh
!Reads 1D mesh
SUBROUTINE readMesh1D(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(mesh1DGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: x
INTEGER:: p(1:2)
INTEGER:: e, et, n, eTemp, elemType, bt
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Open file mesh
OPEN(10, FILE=TRIM(filename))
!Skip header
READ(10, *)
READ(10, *)
READ(10, *)
READ(10, *)
!Read number of nodes
READ(10, *) self%numNodes
!Allocate required matrices and vectors
ALLOCATE(self%nodes(1:self%numNodes))
ALLOCATE(self%K(1:self%numNodes, 1:self%numNodes))
ALLOCATE(self%IPIV(1:self%numNodes, 1:self%numNodes))
self%K = 0.D0
self%IPIV = 0
!Read nodes coordinates. Only relevant for x
DO e = 1, self%numNodes
READ(10, *) n, x
ALLOCATE(meshNode1D:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /))
END DO
!Skips comments
READ(10, *)
READ(10, *)
!Reads the total number of elements (edges+vol)
READ(10, *) totalNumElem
self%numEdges = 0
DO e = 1, totalNumElem
READ(10, *) eTemp, elemType
IF (elemType == 15) THEN !15 is physical node in GMSH
self%numEdges = e
END IF
END DO
!Substract the number of edges to the total number of elements
!to obtain the number of volume elements
self%numVols = totalNumelem - self%numEdges
!Allocates arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beginning of reading elements
DO e = 1, totalNumelem
BACKSPACE(10)
END DO
!Reads edges
DO e = 1, self%numEdges
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
!Associate boundary condition
bt = getBoundaryId(boundaryType)
SELECT CASE(boundary(bt)%obj%boundaryType)
CASE ('reflection')
ALLOCATE(meshEdge1DRef:: self%edges(e)%obj)
CASE ('absorption')
ALLOCATE(meshEdge1DAbs:: self%edges(e)%obj)
END SELECT
CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType)
END DO
!Read and initialize volumes
DO e = 1, self%numVols
READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2)
ALLOCATE(meshVol1DSegm:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2))
END DO
CLOSE(10)
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
DO et = 1, self%numVols
IF (e /= et) THEN
CALL connected(self%vols(e)%obj, self%vols(et)%obj)
END IF
END DO
!Connectivity betwen vols and edges
DO et = 1, self%numEdges
CALL connected(self%vols(e)%obj, self%edges(et)%obj)
END DO
!Constructs the global K matrix
CALL constructGlobalK(self%K, self%vols(e)%obj)
END DO
END SUBROUTINE readMesh1D
SUBROUTINE connectedVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshVol1DSegm)
SELECT TYPE(elemB)
TYPE IS(meshVol1DSegm)
CALL connectedSegmSegm(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolVol
SUBROUTINE connectedSegmSegm(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA
CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END SUBROUTINE connectedSegmSegm
SUBROUTINE connectedVolEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS (meshVol1DSegm)
SELECT TYPE(elemB)
CLASS IS(meshEdge1D)
CALL connectedSegmEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
SUBROUTINE connectedSegmEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA
CLASS(meshEdge1D), INTENT(inout), TARGET:: elemB
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END SUBROUTINE connectedSegmEdge
SUBROUTINE constructGlobalK(K, elem)
IMPLICIT NONE
REAL(8), INTENT(inout):: K(1:,1:)
CLASS(meshVol), INTENT(in):: elem
REAL(8):: localK(1:2,1:2)
INTEGER:: i, j
INTEGER:: n(1:2)
SELECT TYPE(elem)
TYPE IS(meshVol1DSegm)
localK = elem%elemK()
n = (/ elem%n1%n, elem%n2%n /)
CLASS DEFAULT
n = 0
localK = 0.D0
END SELECT
DO i = 1, 2
DO j = 1, 2
K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
END DO
END DO
END SUBROUTINE constructGlobalK
END MODULE moduleMesh1DRead

11
src/modules/mesh/Cyl/makefile
View file

@ -1,11 +0,0 @@
all : moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.o
moduleMeshCyl.o: moduleMeshCyl.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleMeshCylRead.o: moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@

1022
src/modules/mesh/Cyl/moduleMeshCyl.f95
View file

File diff suppressed because it is too large Load diff

80
src/modules/mesh/Cyl/moduleMeshCylBoundary.f95
View file

@ -1,80 +0,0 @@
!moduleMeshCylBoundary: Edge elements for Cylindrical mesh.
MODULE moduleMeshCylBoundary
USE moduleMeshCyl
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylRef
CONTAINS
PROCEDURE, PASS:: fBoundary => reflection
END TYPE meshEdgeCylRef
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAbs
CONTAINS
PROCEDURE, PASS:: fBoundary => absorption
END TYPE meshEdgeCylAbs
TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAxis
CONTAINS
PROCEDURE, PASS:: fBoundary => symmetryAxis
END TYPE meshEdgeCylAxis
CONTAINS
SUBROUTINE reflection(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylRef), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
edgeNorm = DSQRT((self%r(2)-self%r(1))**2 + (self%z(2)-self%z(1))**2)
cosT = (self%z(2)-self%z(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - self%z(1);
rp(2) = part%r(2) - self%r(1);
rpp(1) = cosT*rp(1) - sinT*rp(2)
rpp(2) = sinT*rp(1) + cosT*rp(2)
rpp(2) = -rpp(2)
vpp(1) = cosT*part%v(1) - sinT*part%v(2)
vpp(2) = sinT*part%v(1) + cosT*part%v(2)
vpp(2) = -vpp(2)
part%r(1) = cosT*rpp(1) + sinT*rpp(2) + self%z(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + self%r(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
part%n_in = .TRUE.
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylAbs), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
!TODO: Add scatter to mesh nodes
part%n_in = .FALSE.
END SUBROUTINE absorption
SUBROUTINE symmetryAxis(self, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdgeCylAxis), INTENT(inout):: self
CLASS(particle), INTENT(inout):: part
END SUBROUTINE symmetryAxis
END MODULE moduleMeshCylBoundary

610
src/modules/mesh/Cyl/moduleMeshCylRead.f95
View file

@ -1,610 +0,0 @@
MODULE moduleMeshCylRead
USE moduleMesh
USE moduleMeshCyl
USE moduleMeshCylBoundary
!TODO: make this abstract to allow different mesh formats
TYPE, EXTENDS(meshGeneric):: meshCylGeneric
CONTAINS
PROCEDURE, PASS:: init => initCylMesh
PROCEDURE, PASS:: readMesh => readMeshCylGmsh
END TYPE
INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connected
CONTAINS
!Init mesh
SUBROUTINE initCylMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
SELECT CASE(meshFormat)
CASE ("gmsh")
self%printOutput => printOutputGmsh
self%printColl => printCollGmsh
self%printEM => printEMGmsh
CASE DEFAULT
CALL criticalError("Mesh type " // meshFormat // " not supported.", "initCylMesh")
END SELECT
END SUBROUTINE initCylMesh
!Read mesh from gmsh file
SUBROUTINE readMeshCylGmsh(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(meshCylGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r, z
INTEGER:: p(1:4)
INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Read msh
OPEN(10, FILE=TRIM(filename))
!Skip header
READ(10, *)
READ(10, *)
READ(10, *)
READ(10, *)
!Read number of nodes
READ(10, *) self%numNodes
!Allocate required matrices and vectors
ALLOCATE(self%nodes(1:self%numNodes))
ALLOCATE(self%K(1:self%numNodes,1:self%numNodes))
ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
self%K = 0.D0
self%IPIV = 0
!Read nodes cartesian coordinates (x=z, y=r, z=null)
DO e=1, self%numNodes
READ(10, *) n, z, r
ALLOCATE(meshNodeCyl:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/r, z, 0.D0 /))
END DO
!Skips comments
READ(10, *)
READ(10, *)
!Reads Totalnumber of elements
READ(10, *) TotalnumElem
!counts edges and volume elements
self%numEdges = 0
DO e=1, TotalnumElem
READ(10,*) eTemp, elemType
IF (elemType==1) THEN
self%numEdges=e
END IF
END DO
!Substract the number of edges to the total number of elements
!to obtain the number of volume elements
self%numVols = TotalnumElem - self%numEdges
!Allocates arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beggining to read elements
DO e=1, totalNumElem
BACKSPACE(10)
END DO
!Reads edges
DO e=1, self%numEdges
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
SELECT CASE(boundary(bt)%obj%boundaryType)
CASE ('reflection')
ALLOCATE(meshEdgeCylRef:: self%edges(e)%obj)
CASE ('absorption')
ALLOCATE(meshEdgeCylAbs:: self%edges(e)%obj)
CASE ('axis')
ALLOCATE(meshEdgeCylAxis:: self%edges(e)%obj)
END SELECT
CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType)
END DO
!Read and initialize volumes
DO e=1, self%numVols
READ(10,*) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE (2)
!Triangular element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVolCylTria:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3))
CASE (3)
!Quadrilateral element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVolCylQuad:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
END SELECT
END DO
CLOSE(10)
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
DO et = 1, self%numVols
IF (e /= et) THEN
CALL connected(self%vols(e)%obj, self%vols(et)%obj)
END IF
END DO
!Connectivity between vols and edges
DO et = 1, self%numEdges
CALL connected(self%vols(e)%obj, self%edges(et)%obj)
END DO
!Constructs the global K matrix
CALL constructGlobalK(self%K, self%vols(e)%obj)
END DO
END SUBROUTINE readMeshCylGmsh
!Selects type of elements to build connection
SUBROUTINE connectedVolVol(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshVol), INTENT(inout):: elemB
SELECT TYPE(elemA)
TYPE IS(meshVolCylQuad)
!Element A is a quadrilateral
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
!Element B is a quadrilateral
CALL connectedQuadQuad(elemA, elemB)
TYPE IS(meshVolCylTria)
!Element B is a triangle
CALL connectedQuadTria(elemA, elemB)
END SELECT
TYPE IS(meshVolCylTria)
!Element A is a Triangle
SELECT TYPE(elemB)
TYPE IS(meshVolCylQuad)
!Element B is a quadrilateral
CALL connectedQuadTria(elemB, elemA)
TYPE IS(meshVolCylTria)
!Element B is a triangle
CALL connectedTriaTria(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolVol
SUBROUTINE connectedVolEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: elemA
CLASS(meshEdge), INTENT(inout):: elemB
SELECT TYPE(elemB)
CLASS IS(meshEdgeCyl)
SELECT TYPE(elemA)
TYPE IS(meshVolCylQuad)
!Element A is a quadrilateral
CALL connectedQuadEdge(elemA, elemB)
TYPE IS(meshVolCylTria)
!Element A is a triangle
CALL connectedTriaEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
SUBROUTINE connectedQuadQuad(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n4%n) THEN
elemA%e2 => elemB
elemB%e4 => elemA
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
END IF
END SUBROUTINE connectedQuadQuad
SUBROUTINE connectedQuadTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e4 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedQuadTria
SUBROUTINE connectedTriaTria(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n3%n) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
END IF
END IF
END SUBROUTINE connectedTriaTria
SUBROUTINE connectedQuadEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 4
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
IF (elemA%n4%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e4 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedQuadEdge
SUBROUTINE connectedTriaEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB
!Check direction 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF (elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n1%n) THEN
elemA%e1 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF (elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n1%n) THEN
elemA%e2 => elemB
elemB%e1 => elemA
END IF
END IF
!Check direction 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF (elemA%n3%n == elemB%n1%n .AND. &
elemA%n1%n == elemB%n2%n) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
elemA%n1%n == elemB%n1%n) THEN
elemA%e3 => elemB
elemB%e1 => elemA
END IF
END IF
END SUBROUTINE connectedTriaEdge
SUBROUTINE constructGlobalK(K, elem)
IMPLICIT NONE
REAL(8), INTENT(inout):: K(1:,1:)
CLASS(meshVol), INTENT(in):: elem
REAL(8), ALLOCATABLE:: localK(:,:)
INTEGER:: nNodes, i, j
INTEGER, ALLOCATABLE:: n(:)
SELECT TYPE(elem)
TYPE IS(meshVolCylQuad)
nNodes = 4
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, &
elem%n3%n, elem%n4%n /)
TYPE IS(meshVolCylTria)
nNodes = 3
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, elem%n3%n /)
CLASS DEFAULT
nNodes = 0
ALLOCATE(localK(1:1, 1:1))
localK = 0.D0
ALLOCATE(n(1:1))
n = 0
END SELECT
DO i = 1, nNodes
DO j = 1, nNodes
K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
END DO
END DO
END SUBROUTINE constructGlobalK
END MODULE moduleMeshCylRead

10
src/modules/mesh/makefile
View file

@ -1,11 +1,11 @@
all: moduleMesh.o Cyl.o 1DCart.o
all: moduleMesh.o 2DCyl.o 1DCart.o
Cyl.o:
$(MAKE) -C Cyl all
2DCyl.o:
$(MAKE) -C 2DCyl all
1DCart.o:
$(MAKE) -C 1DCart all
moduleMesh.o: moduleMesh.f95
$(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@
moduleMesh.o: moduleMesh.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

1
src/modules/mesh/moduleMesh.f90
View file

@ -240,7 +240,6 @@ MODULE moduleMesh
PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL()
PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL()
CONTAINS
PROCEDURE(initMesh_interface), DEFERRED, PASS:: init
PROCEDURE(readMesh_interface), DEFERRED, PASS:: readMesh

606
src/modules/mesh/moduleMesh.f95
View file

@ -1,606 +0,0 @@
!moduleMesh: General module for Finite Element mesh
MODULE moduleMesh
USE moduleList
USE moduleOutput
IMPLICIT NONE
!Parent of Node element
TYPE, PUBLIC, ABSTRACT:: meshNode
!Node index
INTEGER:: n = 0
!Node volume
REAL(8):: v = 0.D0
!Output values
TYPE(outputNode), ALLOCATABLE:: output(:)
TYPE(emNode):: emData
CONTAINS
PROCEDURE(initNode_interface), DEFERRED, PASS:: init
PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates
END TYPE meshNode
ABSTRACT INTERFACE
!Interface of init a node (3D generic coordinates)
SUBROUTINE initNode_interface(self, n, r)
IMPORT:: meshNode
CLASS(meshNode), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
END SUBROUTINE initNode_interface
!Interface to get coordinates from node
PURE FUNCTION getCoord_interface(self) RESULT(r)
IMPORT:: meshNode
CLASS(meshNode), INTENT(in):: self
REAL(8):: r(1:3)
END FUNCTION getCoord_interface
END INTERFACE
!Containers for nodes in the mesh
TYPE:: meshNodeCont
CLASS(meshNode), ALLOCATABLE:: obj
END TYPE meshNodeCont
!Parent of Edge element
TYPE, PUBLIC, ABSTRACT:: meshEdge
!Element index
INTEGER:: n = 0
!Connectivity to vols
CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL()
!Normal vector
REAL(8):: normal(1:3)
!Physical surface in mesh
INTEGER:: physicalSurface
!id for boundary condition
INTEGER:: bt = 0
CONTAINS
PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary
PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
PROCEDURE(randPos_interface), DEFERRED, PASS:: randPos
END TYPE meshEdge
ABSTRACT INTERFACE
SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
INTEGER, INTENT(in):: bt
INTEGER, INTENT(in):: physicalSurface
END SUBROUTINE initEdge_interface
SUBROUTINE boundary_interface(self, part)
USE moduleSpecies
IMPORT:: meshEdge
CLASS (meshEdge), INTENT(inout):: self
CLASS (particle), INTENT(inout):: part
END SUBROUTINE
PURE FUNCTION getNodesEdge_interface(self) RESULT(n)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
END FUNCTION
FUNCTION randPos_interface(self) RESULT(r)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self
REAL(8):: r(1:3)
END FUNCTION randPos_interface
END INTERFACE
!Containers for edges in the mesh
TYPE:: meshEdgeCont
CLASS(meshEdge), ALLOCATABLE:: obj
END TYPE meshEdgeCont
!Parent of Volume element
TYPE, PUBLIC, ABSTRACT:: meshVol
!Volume index
INTEGER:: n = 0
!Maximum collision rate
REAL(8):: sigmaVrelMax = 1.D-15
!Volume
REAL(8):: volume = 0.D0
!List of particles inside the volume
TYPE(listNode):: listPart_in
!Lock indicator for listPart_in
INTEGER(KIND=OMP_LOCK_KIND):: lock
!Number of collisions per volume
INTEGER:: nColl = 0
!Total weight of particles inside cell
REAL(8):: totalWeight = 0.D0
CONTAINS
PROCEDURE(initVol_interface), DEFERRED, PASS:: init
PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter
PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF
PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes
PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF
PROCEDURE, PASS:: findCell
PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log
PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside
PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement
PROCEDURE, PASS:: collision
PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput
END TYPE meshVol
ABSTRACT INTERFACE
SUBROUTINE initVol_interface(self, n, p)
IMPORT:: meshVol
CLASS(meshVol), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
END SUBROUTINE initVol_interface
SUBROUTINE scatter_interface(self, part)
USE moduleSpecies
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
CLASS(particle), INTENT(in):: part
END SUBROUTINE scatter_interface
PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: EF(1:3)
END FUNCTION gatherEF_interface
PURE FUNCTION getNodesVol_interface(self) RESULT(n)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
END FUNCTION getNodesVol_interface
PURE FUNCTION elemF_interface(self, source) RESULT(localF)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
END FUNCTION elemF_interface
SUBROUTINE nextElement_interface(self, xi, nextElement)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), POINTER, INTENT(out):: nextElement
END SUBROUTINE nextElement_interface
PURE FUNCTION phy2log_interface(self,r) RESULT(xN)
IMPORT:: meshVol
CLASS(meshVol), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xN(1:3)
END FUNCTION phy2log_interface
PURE FUNCTION inside_interface(xi) RESULT(ins)
IMPORT:: meshVol
REAL(8), INTENT(in):: xi(1:3)
LOGICAL:: ins
END FUNCTION inside_interface
SUBROUTINE collision_interface(self)
IMPORT:: meshVol
CLASS(meshVol), INTENT(inout):: self
END SUBROUTINE collision_interface
PURE SUBROUTINE resetOutput_interface(self)
IMPORT:: meshVol
CLASS(meshVol), INTENT(inout):: self
END SUBROUTINE resetOutput_interface
END INTERFACE
!Containers for volumes in the mesh
TYPE:: meshVolCont
CLASS(meshVol), ALLOCATABLE:: obj
END TYPE meshVolCont
!Abstract type of mesh
TYPE, PUBLIC, ABSTRACT:: meshGeneric
INTEGER:: numEdges, numNodes, numVols
!Array of nodes
TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
!Array of boundary elements
TYPE(meshEdgeCont), ALLOCATABLE:: edges(:)
!Array of volume elements
TYPE(meshVolCont), ALLOCATABLE:: vols(:)
!Global stiffness matrix
REAL(8), ALLOCATABLE, DIMENSION(:,:):: K
!Permutation matrix for P L U factorization
INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL()
PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL()
PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL()
CONTAINS
PROCEDURE(initMesh_interface), DEFERRED, PASS:: init
PROCEDURE(readMesh_interface), DEFERRED, PASS:: readMesh
END TYPE meshGeneric
ABSTRACT INTERFACE
!Inits the mesh
SUBROUTINE initMesh_interface(self, meshFormat)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
END SUBROUTINE initMesh_interface
!Reads the mesh from a file
SUBROUTINE readMesh_interface(self, filename)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
END SUBROUTINE readMesh_interface
!Prints Species data
SUBROUTINE printOutput_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printOutput_interface
!Prints number of collisions
SUBROUTINE printColl_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printColl_interface
!Prints EM info
SUBROUTINE printEM_interface(self, t)
IMPORT meshGeneric
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
END SUBROUTINE printEM_interface
END INTERFACE
!Generic mesh
CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh
CONTAINS
!Find next cell for particle
RECURSIVE SUBROUTINE findCell(self, part, oldCell)
USE moduleSpecies
USE OMP_LIB
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self
CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell
CLASS(particle), INTENT(inout), TARGET:: part
REAL(8):: xi(1:3)
CLASS(*), POINTER:: nextElement
xi = self%phy2log(part%r)
!Checks if particle is inside 'self' cell
IF (self%inside(xi)) THEN
part%vol = self%n
part%xi = xi
part%n_in = .TRUE.
!Assign particle to listPart_in
CALL OMP_SET_LOCK(self%lock)
CALL self%listPart_in%add(part)
self%totalWeight = self%totalWeight + part%weight
CALL OMP_UNSET_LOCK(self%lock)
ELSE
!If not, searches for a neighbour and repeats the process.
CALL self%nextElement(xi, nextElement)
!Defines the next step
SELECT TYPE(nextElement)
CLASS IS(meshVol)
!Particle moved to new cell, repeat find procedure
CALL nextElement%findCell(part, self)
CLASS IS (meshEdge)
!Particle encountered an edge, execute boundary
CALL nextElement%fBoundary(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(*,*) "ERROR, CHECK findCell"
END SELECT
END IF
END SUBROUTINE findCell
!Computes collisions in element
SUBROUTINE collision(self)
USE moduleCollisions
USE moduleSpecies
USE moduleList
use moduleRefParam
IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self
INTEGER:: nPart !Number of particles inside the cell
REAL(8):: pMax !Maximum probability of collision
INTEGER:: rnd !random index
TYPE(particle), POINTER:: part_i, part_j
INTEGER:: n !collision
INTEGER:: ij, k
REAL(8):: sigmaVrelMaxNew
TYPE(pointerArray), ALLOCATABLE:: partTemp(:)
self%nColl = 0
nPart = self%listPart_in%amount
IF (nPart > 1) THEN
pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume
self%nColl = INT(REAL(nPart)*pMax*0.5D0)
!Converts the list of particles to an array for easy access
IF (self%nColl > 0) THEN
partTemp = self%listPart_in%convert2Array()
END IF
DO n = 1, self%nColl
!Select random numbers
rnd = 1 + FLOOR(nPart*RAND())
part_i => partTemp(rnd)%part
rnd = 1 + FLOOR(nPart*RAND())
part_j => partTemp(rnd)%part
ij = interactionIndex(part_i%sp, part_j%sp)
sigmaVrelMaxNew = 0.D0
DO k = 1, interactionMatrix(ij)%amount
CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j)
END DO
!Update maximum cross section*v_rel per each collision
IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN
self%sigmaVrelMax = sigmaVrelMaxNew
END IF
END DO
END IF
self%totalWeight = 0.D0
!Reset output in nodes
CALL self%resetOutput()
!Erase the list of particles inside the cell
CALL self%listPart_in%erase()
END SUBROUTINE collision
SUBROUTINE printOutputGmsh(self, t)
USE moduleRefParam
USE moduleSpecies
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, i
TYPE(outputFormat):: output(1:self%numNodes)
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
time = DBLE(t)*tauMin*ti_ref
DO i = 1, nSpecies
WRITE(tstring, '(I6.6)') t
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)") '"Density (m^-3)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) 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)") '"Velocity (m/s)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 3
WRITE(60, *) 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)") '"Pressure (Pa)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) 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)") '"Temperature (K)"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numNodes
DO n=1, self%numNodes
WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
END DO
WRITE(60, "(A)") '$EndNodeData'
CLOSE (60)
END DO
END SUBROUTINE printOutputGmsh
SUBROUTINE printCollGmsh(self, t)
USE moduleRefParam
USE moduleCaseParam
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
IF (collOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tstring, '(I6.6)') t
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'
WRITE(60, "(A)") '$ElementData'
WRITE(60, "(A)") '1'
WRITE(60, "(A)") '"Collisions"'
WRITE(60, *) 1
WRITE(60, *) time
WRITE(60, *) 3
WRITE(60, *) t
WRITE(60, *) 1
WRITE(60, *) self%numVols
DO n=1, self%numVols
WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl
END DO
WRITE(60, "(A)") '$EndElementData'
CLOSE(60)
END IF
END SUBROUTINE printCollGmsh
SUBROUTINE printEMGmsh(self, t)
USE moduleRefParam
USE moduleCaseParam
USE moduleOutput
IMPLICIT NONE
CLASS(meshGeneric), INTENT(in):: self
INTEGER, INTENT(in):: t
INTEGER:: n, e
REAL(8):: time
CHARACTER(:), ALLOCATABLE:: fileName
CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
REAL(8):: xi(1:3)
xi = (/ 0.D0, 0.D0, 0.D0 /)
IF (emOutput) THEN
time = DBLE(t)*tauMin*ti_ref
WRITE(tstring, '(I6.6)') t
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
DO n=1, self%numNodes
WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
END DO
WRITE(20, "(A)") '$EndNodeData'
WRITE(20, "(A)") '$ElementData'
WRITE(20, "(A)") '1'
WRITE(20, "(A)") '"Electric Field (V/m)"'
WRITE(20, *) 1
WRITE(20, *) time
WRITE(20, *) 3
WRITE(20, *) t
WRITE(20, *) 3
WRITE(20, *) 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'
CLOSE(20)
END IF
END SUBROUTINE printEMGmsh
END MODULE moduleMesh

36
src/modules/moduleBoundary.f95
View file

@ -1,36 +0,0 @@
MODULE moduleBoundary
TYPE, PUBLIC:: boundaryGeneric
INTEGER:: id = 0
CHARACTER(:), ALLOCATABLE:: name
INTEGER:: physicalSurface = 0
CHARACTER(:), ALLOCATABLE:: boundaryType
END TYPE boundaryGeneric
TYPE:: boundaryCont
CLASS(boundaryGeneric), ALLOCATABLE:: obj
END TYPE boundaryCont
INTEGER:: nBoundary = 0
TYPE(boundaryCont), ALLOCATABLE:: 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)%obj%physicalSurface) id = boundary(i)%obj%id
END DO
END FUNCTION getBoundaryId
END MODULE moduleBoundary

9
src/modules/moduleCaseParam.f95
View file

@ -1,9 +0,0 @@
!Problems of the case
MODULE moduleCaseParam
!Maximum number of iterations and number of species
INTEGER:: tmax
REAL(8), ALLOCATABLE:: tau(:)
REAL(8):: tauMin
END MODULE moduleCaseParam

171
src/modules/moduleCollisions.f95
View file

@ -1,171 +0,0 @@
MODULE moduleCollisions
USE moduleTable
TYPE, ABSTRACT:: collisionBinary
REAL(8):: rMass !reduced mass
TYPE(table1D):: crossSec !cross section of collision
CONTAINS
PROCEDURE(initBinary_interface), PASS, DEFERRED:: init
PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide
END TYPE collisionBinary
ABSTRACT INTERFACE
SUBROUTINE initBinary_interface(self, crossSectionFilename, mass_i, mass_j)
IMPORT:: collisionBinary
CLASS(collisionBinary), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
REAL(8), INTENT(in):: mass_i, mass_j
END SUBROUTINE
SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j)
USE moduleSpecies
IMPORT:: collisionBinary
CLASS(collisionBinary), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew
TYPE(particle), INTENT(inout):: part_i, part_j
END SUBROUTINE
END INTERFACE
!Container for binary collisions
TYPE:: collisionCont
CLASS(collisionBinary), ALLOCATABLE:: obj
END TYPE collisionCont
TYPE, EXTENDS(collisionBinary):: collisionBinaryElastic
!Weight distribution for Maxwellian function
REAL(8):: w_i = (1.D0+DSQRT(3.D0))/2.D0
REAL(8):: w_j = (DSQRT(3.D0)-1.D0)/2.D0
CONTAINS
PROCEDURE, PASS:: init => initBinaryElastic
PROCEDURE, PASS:: collide => collideBinaryElastic
END TYPE collisionBinaryElastic
!Type for interaction matrix
TYPE:: interactionsBinary
INTEGER:: amount
TYPE(collisionCont), ALLOCATABLE:: collisions(:)
CONTAINS
PROCEDURE, PASS:: init => initInteractionBinary
END TYPE interactionsBinary
!Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory
TYPE(interactionsBinary), ALLOCATABLE:: interactionMatrix(:)
!Folder for collision cross section tables
CHARACTER(:), ALLOCATABLE:: pathCollisions
CONTAINS
SUBROUTINE initInteractionMatrix(interactionMatrix)
USE moduleSpecies
IMPLICIT NONE
TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:)
INTEGER:: nInteractions
nInteractions = (nSpecies*(nSpecies+1))/2
ALLOCATE(interactionMatrix(1:nInteractions))
END SUBROUTINE initInteractionMatrix
FUNCTION interactionIndex(i,j) RESULT(k)
INTEGER:: i, j
INTEGER:: p
INTEGER:: k
k = i + j
p = (k + ABS(i - j))/2
k = k + (p*(p-3))/2
END FUNCTION interactionIndex
SUBROUTINE initInteractionBinary(self, amount)
IMPLICIT NONE
CLASS(interactionsBinary), INTENT(inout):: self
INTEGER, INTENT(in):: amount
INTEGER:: k
self%amount = amount
ALLOCATE(self%collisions(1:self%amount))
DO k= 1, self%amount
!TODO: make type dependent
ALLOCATE(collisionBinaryElastic:: self%collisions(k)%obj)
END DO
END SUBROUTINE initInteractionBinary
SUBROUTINE initBinaryElastic(self, crossSectionFilename, mass_i, mass_j)
USE moduleTable
USE moduleRefParam
USE moduleConstParam
IMPLICIT NONE
CLASS(collisionBinaryElastic), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename
REAL(8), INTENT(in):: mass_i, mass_j
!Reads data from file
CALL self%crossSec%init(crossSectionFilename)
!Convert to no-dimensional units
CALL self%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
!Calculates reduced mass
self%rMass = (mass_i*mass_j)/(mass_i+mass_j)
END SUBROUTINE
SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, &
part_i, part_j)
USE moduleConstParam
USE moduleSpecies
USE moduleTable
IMPLICIT NONE
CLASS(collisionBinaryElastic), INTENT(in):: self
REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew
TYPE(particle), INTENT(inout):: part_i, part_j
REAL(8):: sigmaVrel
REAL(8):: vRel !relative velocity
REAL(8):: eRel !relative energy
REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities
REAL(8):: v_ij !sum of velocities modules
REAL(8):: alpha !random angle of scattering
!eRel (in units of [m][L]^2[s]^-2
vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1
eRel = self%rMass*vRel**2
sigmaVrel = self%crossSec%get(eRel)*vRel
sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel
IF (sigmaVrelMaxNew/sigmaVrelMax > RAND()) THEN
!Applies the collision
v_i = NORM2(part_i%v)
v_j = NORM2(part_j%v)
v_ij = v_i+v_j
vp_j = v_ij*self%w_i
vp_i = v_ij*self%w_j
alpha = PI*RAND()
part_i%v(1) = v_i*DCOS(alpha)
part_i%v(2) = v_i*DSIN(alpha)
alpha = PI*RAND()
part_j%v(1) = v_j*DCOS(alpha)
part_j%v(2) = v_j*DSIN(alpha)
END IF
END SUBROUTINE
END MODULE moduleCollisions

11
src/modules/moduleCompTime.f95
View file

@ -1,11 +0,0 @@
!Information to calculate computation time
MODULE moduleCompTime
REAL(8):: tStep=0.D0
REAL(8):: tPush=0.D0
REAL(8):: tReset=0.D0
REAL(8):: tColl=0.D0
REAL(8):: tWeight=0.D0
REAL(8):: tEMField=0.D0
END MODULE moduleCompTime

15
src/modules/moduleConstParam.f95
View file

@ -1,15 +0,0 @@
!Physical and mathematical constants
MODULE moduleConstParam
IMPLICIT NONE
PUBLIC
REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi
REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1
REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge
REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI
REAL(8), PARAMETER:: eV2J = qe !Electron volt to Joule conversion
REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0
END MODULE moduleConstParam

128
src/modules/moduleEM.f95
View file

@ -1,128 +0,0 @@
!Module to solve the electromagnetic field
MODULE moduleEM
IMPLICIT NONE
TYPE:: boundaryEM
CHARACTER(:), ALLOCATABLE:: typeEM
INTEGER:: physicalSurface
REAL(8):: potential
CONTAINS
PROCEDURE, PASS:: apply
END TYPE boundaryEM
INTEGER:: nBoundaryEM
TYPE(boundaryEM), ALLOCATABLE:: boundEM(:)
!Information of charge and reference parameters for rho vector
REAL(8), ALLOCATABLE:: qSpecies(:)
CONTAINS
SUBROUTINE apply(self, edge)
USE moduleMesh
IMPLICIT NONE
CLASS(boundaryEM), INTENT(in):: self
CLASS(meshEdge):: edge
INTEGER, ALLOCATABLE:: nodes(:)
INTEGER:: n
nodes = edge%getNodes()
DO n = 1, SIZE(nodes)
SELECT CASE(self%typeEM)
CASE ("dirichlet")
mesh%K(nodes(n), :) = 0.D0
mesh%K(nodes(n), nodes(n)) = 1.D0
mesh%nodes(nodes(n))%obj%emData%type = self%typeEM
mesh%nodes(nodes(n))%obj%emData%phi = self%potential
END SELECT
END DO
END SUBROUTINE
PURE FUNCTION gatherElecField(part) RESULT(elField)
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
TYPE(particle), INTENT(in):: part
REAl(8):: xi(1:3) !Logical coordinates of particle in element
REAL(8):: elField(1:3)
elField = 0.D0
xi = part%xi
elField = mesh%vols(part%vol)%obj%gatherEF(xi)
END FUNCTION gatherElecField
SUBROUTINE assembleSourceVector(vectorF)
USE moduleMesh
USE moduleRefParam
IMPLICIT NONE
REAL(8), INTENT(out):: vectorF(1:mesh%numNodes)
REAL(8), ALLOCATABLE:: localF(:)
INTEGER, ALLOCATABLE:: nodes(:)
REAL(8), ALLOCATABLE:: rho(:)
INTEGER:: nNodes
INTEGER:: e, i, ni
CLASS(meshNode), POINTER:: node
!$OMP SINGLE
vectorF = 0.D0
!$OMP END SINGLE
!$OMP DO REDUCTION(+:vectorF)
DO e = 1, mesh%numVols
nodes = mesh%vols(e)%obj%getNodes()
nNodes = SIZE(nodes)
!Calculates charge density (rho) in element nodes
ALLOCATE(rho(1:nNodes))
rho = 0.D0
DO i = 1, nNodes
ni = nodes(i)
node => mesh%nodes(ni)%obj
rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref))
END DO
!Calculates local F vector
localF = mesh%vols(e)%obj%elemF(rho)
!Assign local F to global F
DO i = 1, nNodes
ni = nodes(i)
vectorF(ni) = vectorF(ni) + localF(i)
END DO
DEALLOCATE(localF)
DEALLOCATE(nodes, rho)
END DO
!$OMP END DO
!Apply boundary conditions
!$OMP DO
DO i = 1, mesh%numNodes
node => mesh%nodes(i)%obj
SELECT CASE(node%emData%type)
CASE ("dirichlet")
vectorF(i) = node%emData%phi
END SELECT
END DO
!$OMP END DO
END SUBROUTINE assembleSourceVector
END MODULE moduleEM

43
src/modules/moduleErrors.f95
View file

@ -1,43 +0,0 @@
!moduleErrors: Manages the different type of errors the program can produce.
! By errors we understand critical errors (that stop the program),
! warnings (that only output a message with a WARNING tag),
! o verbose outputs that can be used for debugging porpouse.
MODULE moduleErrors
CONTAINS
SUBROUTINE criticalError(msg, pgr)
IMPLICIT NONE
CHARACTER(*), INTENT(in):: msg, pgr
CHARACTER(:), ALLOCATABLE:: errorMsg
errorMsg = "CRITICAL error in " // pgr // " with message:" // NEW_LINE('A') // msg
ERROR STOP errorMsg
END SUBROUTINE criticalError
SUBROUTINE warningError(msg)
IMPLICIT NONE
CHARACTER(*), INTENT(in):: msg
CHARACTER(:), ALLOCATABLE:: errorMsg
errorMsg = "WARNING: " // msg
WRITE (*, '(A)') errorMsg
END SUBROUTINE warningError
SUBROUTINE verboseError(msg)
IMPLICIT NONE
CHARACTER(*), INTENT(in):: msg
CHARACTER(:), ALLOCATABLE:: errorMsg
errorMsg = msg
WRITE (*, '(A)') errorMsg
END SUBROUTINE verboseError
END MODULE moduleErrors

110
src/modules/moduleInject.f90
View file

@ -1,21 +1,63 @@
!injection of particles
MODULE moduleInject
!Generic type for velocity distribution function
TYPE, ABSTRACT:: velDistGeneric
CONTAINS
!Returns random velocity from distribution function
PROCEDURE(randomVel_interface), DEFERRED, PASS:: randomVel
END TYPE velDistGeneric
ABSTRACT INTERFACE
FUNCTION randomVel_interface(self) RESULT(v)
IMPORT velDistGeneric
CLASS(velDistGeneric), INTENT(in):: self
REAL(8):: v
END FUNCTION randomVel_interface
END INTERFACE
!Container for velocity distributions
TYPE:: velDistCont
CLASS(velDistGeneric), ALLOCATABLE:: obj
END TYPE velDistCont
!Maxwellian distribution function
TYPE, EXTENDS(velDistGeneric):: velDistMaxwellian
REAL(8):: v !Velocity
REAL(8):: vTh !Thermal Velocity
CONTAINS
PROCEDURE, PASS:: randomVel => randomVelMaxwellian
END TYPE velDistMaxwellian
!Dirac's delta distribution function
TYPE, EXTENDS(velDistGeneric):: velDistDelta
REAL(8):: v !Velocity
CONTAINS
PROCEDURE, PASS:: randomVel => randomVelDelta
END TYPE velDistDelta
!Generic injection of particles
TYPE:: injectGeneric
INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: vMod !Velocity (module)
REAL(8):: T(1:3) !Temperature
REAL(8):: v(1:3) !Velocity(vector)
REAL(8):: vTh(1:3) !Thermal Velocity
REAL(8):: n(1:3) !Direction of injection
INTEGER:: nParticles !Number of particles to introduce each time step
INTEGER:: sp !Species of injection
INTEGER:: nEdges
INTEGER, ALLOCATABLE:: edges(:) !Array with edges
TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
CONTAINS
PROCEDURE, PASS:: init => initInject
PROCEDURE, PASS:: addParticles => addParticlesMaxwellian
PROCEDURE, PASS:: addParticles
END TYPE injectGeneric
@ -64,9 +106,6 @@ MODULE moduleInject
self%nParticles = self%nParticles * solver%pusher(sp)%every
self%sp = sp
self%v = self%vMod * self%n
self%vTh = DSQRT(self%T/species(self%sp)%obj%m)
!Gets the edge elements from which particles are injected
!TODO: Improve this A LOT
DO e = 1, mesh%numEdges
@ -88,7 +127,7 @@ MODULE moduleInject
END DO
! self%sumWeight = SUM(self%weight)
END SUBROUTINE
END SUBROUTINE initInject
!Injection of particles
SUBROUTINE doInjects()
@ -114,12 +153,35 @@ MODULE moduleInject
END SUBROUTINE doInjects
!Random velocity from maxwellian distribution
REAL(8) FUNCTION vBC(u, vth)
SUBROUTINE initVelDistMaxwellian(velDist, v, T, m)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: v, T, m
velDist = velDistMaxwellian(v = v, vTh = DSQRT(T/m))
END SUBROUTINE initVelDistMaxwellian
SUBROUTINE initVelDistDelta(velDist, v)
IMPLICIT NONE
CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
REAL(8), INTENT(in):: v
velDist = velDistDelta(v = v)
END SUBROUTINE initVelDistDelta
!Random velocity from Maxwellian distribution
FUNCTION randomVelMaxwellian(self) RESULT (v)
USE moduleConstParam, ONLY: PI
REAL(8), INTENT (in):: u, vth
IMPLICIT NONE
CLASS(velDistMaxwellian), INTENT(in):: self
REAL(8):: v
REAL(8):: x, y
vBC=0.D0
v = 0.D0
x = 0.D0
DO WHILE (x == 0.D0)
@ -127,11 +189,23 @@ MODULE moduleInject
END DO
CALL RANDOM_NUMBER(y)
vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
v = self%v + self%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
END FUNCTION vBC
END FUNCTION randomVelMaxwellian
SUBROUTINE addParticlesMaxwellian(self)
!Random velocity from Dirac's delta distribution
PURE FUNCTION randomVelDelta(self) RESULT(v)
IMPLICIT NONE
CLASS(velDistDelta), INTENT(in):: self
REAL(8):: v
v = self%v
END FUNCTION randomVelDelta
!Add particles for the injection
SUBROUTINE addParticles(self)
USE moduleSpecies
USE moduleSolver
USE moduleMesh
@ -185,9 +259,9 @@ MODULE moduleInject
END IF
partInj(n)%v = (/ vBC(self%v(1), self%vTh(1)), &
vBC(self%v(2), self%vTh(2)), &
vBC(self%v(3), self%vTh(3)) /)
partInj(n)%v = (/ self%v(1)%obj%randomVel(), &
self%v(2)%obj%randomVel(), &
self%v(3)%obj%randomVel() /)
!Push new particle
CALL solver%pusher(self%sp)%pushParticle(partInj(n))
@ -197,6 +271,6 @@ MODULE moduleInject
END DO
!$OMP END DO
END SUBROUTINE addParticlesMaxwellian
END SUBROUTINE addParticles
END MODULE moduleInject

202
src/modules/moduleInject.f95
View file

@ -1,202 +0,0 @@
!injection of particles
MODULE moduleInject
TYPE:: injectGeneric
INTEGER:: id
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: vMod !Velocity (module)
REAL(8):: T(1:3) !Temperature
REAL(8):: v(1:3) !Velocity(vector)
REAL(8):: vTh(1:3) !Thermal Velocity
REAL(8):: n(1:3) !Direction of injection
INTEGER:: nParticles !Number of particles to introduce each time step
INTEGER:: sp !Species of injection
INTEGER:: nEdges
INTEGER, ALLOCATABLE:: edges(:) !Array with edges
! REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection
! REAL(8):: sumWeight
CONTAINS
PROCEDURE, PASS:: init => initInject
PROCEDURE, PASS:: addParticles => addParticlesMaxwellian
END TYPE injectGeneric
INTEGER:: nInject
TYPE(injectGeneric), ALLOCATABLE:: inject(:)
CONTAINS
!Initialize an injection of particles
SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
USE moduleMesh
USE moduleRefParam
USE moduleConstParam
USE moduleSpecies
USE moduleSolver
USE moduleErrors
IMPLICIT NONE
CLASS(injectGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: i
REAL(8), INTENT(in):: v, n(1:3), T(1:3)
INTEGER, INTENT(in):: sp, physicalSurface
REAL(8), INTENT(in):: flow
CHARACTER(:), ALLOCATABLE, INTENT(in):: units
INTEGER:: e, et
INTEGER:: phSurface(1:mesh%numEdges)
self%id = i
self%vMod = v/v_ref
self%n = n
self%T = T/T_ref
SELECT CASE(units)
CASE ("sccm")
!Standard cubic centimeter per minute
self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight)
CASE ("A")
!Input current in Ampers
self%nParticles = INT(flow*tauMin*ti_ref/(qe*species(sp)%obj%weight))
CASE DEFAULT
CALL criticalError("No support for units: " // units, 'initInject')
END SELECT
IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
self%nParticles = self%nParticles * solver%pusher(sp)%every
self%sp = sp
self%v = self%vMod * self%n
self%vTh = DSQRT(self%T/species(self%sp)%obj%m)
!Gets the edge elements from which particles are injected
!TODO: Improve this A LOT
DO e = 1, mesh%numEdges
phSurface(e) = mesh%edges(e)%obj%physicalSurface
END DO
self%nEdges = COUNT(phSurface == physicalSurface)
ALLOCATE(inject(i)%edges(1:self%nEdges))
! ALLOCATE(inject(i)%weight(1:self%nEdges))
et = 0
DO e=1, mesh%numEdges
IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
et = et + 1
self%edges(et) = mesh%edges(e)%obj%n
END IF
END DO
! self%sumWeight = SUM(self%weight)
END SUBROUTINE
!Injection of particles
SUBROUTINE doInjects()
USE moduleSpecies
USE moduleSolver
IMPLICIT NONE
INTEGER:: i
!$OMP SINGLE
nPartInj = 0
DO i = 1, nInject
IF (solver%pusher(inject(i)%sp)%pushSpecies) nPartInj = nPartInj + inject(i)%nParticles
END DO
IF (ALLOCATED(partInj)) DEALLOCATE(partInj)
ALLOCATE(partInj(1:nPartInj))
!$OMP END SINGLE
DO i=1, nInject
IF (solver%pusher(inject(i)%sp)%pushSpecies) CALL inject(i)%addParticles()
END DO
END SUBROUTINE doInjects
!Random velocity from maxwellian distribution
REAL(8) FUNCTION vBC(u, vth)
USE moduleConstParam, ONLY: PI
REAL(8), INTENT (in):: u, vth
REAL(8):: x, y
vBC=0.D0
x = 0.D0
DO WHILE (x == 0.D0)
x=RAND()
END DO
y=RAND()
vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y)
END FUNCTION vBC
SUBROUTINE addParticlesMaxwellian(self)
USE moduleSpecies
USE moduleSolver
USE moduleMesh
IMPLICIT NONE
CLASS(injectGeneric), INTENT(in):: self
INTEGER:: randomX
INTEGER:: i!, j
INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array
INTEGER:: n
CLASS(meshEdge), POINTER:: randomEdge
!Insert particles
!$OMP SINGLE
nMin = 0
DO i = 1, self%id - 1
IF (solver%pusher(inject(i)%sp)%pushSpecies) nMin = nMin + inject(i)%nParticles
END DO
nMin = nMin + 1
nMax = nMin + self%nParticles - 1
!Assign particle type
partInj(nMin:nMax)%sp = self%sp
!Assign weight to particle.
partInj(nMin:nMax)%weight = species(self%sp)%obj%weight
!Particle is considered to be outside the domain
partInj(nMin:nMax)%n_in = .FALSE.
!Assign charge/mass to charged particle.
SELECT TYPE(sp => species(self%sp)%obj)
TYPE IS(speciesCharged)
partInj(nMin:nMax)%qm = sp%qm
END SELECT
!$OMP END SINGLE
!$OMP DO
DO n = nMin, nMax
randomX = INT(DBLE(self%nEdges-1)*RAND()) + 1
randomEdge => mesh%edges(self%edges(randomX))%obj
!Random position in edge
partInj(n)%r = randomEdge%randPos()
!Volume associated to the edge:
IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN
partInj(n)%vol = randomEdge%e1%n
ELSE
partInj(n)%vol = randomEdge%e2%n
END IF
partInj(n)%v = (/ vBC(self%v(1), self%vTh(1)), &
vBC(self%v(2), self%vTh(2)), &
vBC(self%v(3), self%vTh(3)) /)
!Push new particle
CALL solver%pusher(self%sp)%pushParticle(partInj(n))
!Assign cell to new particle
CALL solver%updateParticleCell(partInj(n))
END DO
!$OMP END DO
END SUBROUTINE addParticlesMaxwellian
END MODULE moduleInject

48
src/modules/moduleInput.f90
View file

@ -520,6 +520,8 @@ MODULE moduleInput
CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
CALL readVelDistr(config, inject(i), object)
END DO
!Allocate array for injected particles
@ -530,6 +532,52 @@ MODULE moduleInput
END SUBROUTINE readInject
!Reads the velocity distribution functions for each inject
SUBROUTINE readVelDistr(config, inj, object)
USE moduleErrors
USE moduleInject
USE moduleSpecies
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
TYPE(injectGeneric), INTENT(inout):: inj
CHARACTER(:), ALLOCATABLE, INTENT(in):: object
INTEGER:: i
CHARACTER(2):: istring
CHARACTER(:), ALLOCATABLE:: fvType
LOGICAL:: found
REAL(8):: v, T, m
!Reads species mass
m = species(inj%sp)%obj%m
!Reads distribution functions for velocity
DO i = 1, 3
WRITE(istring, '(i2)') i
CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found)
IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // &
" found for " // object, 'readVelDistr')
SELECT CASE(fvType)
CASE ("Maxwellian")
v = inj%vMod*inj%n(i)
T = inj%T(i)
CALL initVelDistMaxwellian(inj%v(i)%obj, v, t, m)
CASE ("Delta")
v = inj%vMod*inj%n(i)
CALL initVelDistDelta(inj%v(i)%obj, v)
CASE DEFAULT
CALL criticalError("No velocity distribution type " // fvType // " defined", 'readVelDistr')
END SELECT
END DO
END SUBROUTINE readVelDistr
!Reads configuration for parallel run
SUBROUTINE readParallel(config)
USE moduleParallel
USE moduleErrors

559
src/modules/moduleInput.f95
View file

@ -1,559 +0,0 @@
! moduleInput: Reads JSON configuration file
MODULE moduleInput
USE json_module
IMPLICIT NONE
CONTAINS
!Main routine to read the input JSON file
SUBROUTINE readConfig(inputFile)
USE json_module
USE moduleErrors
USE moduleBoundary
USE moduleInject
IMPLICIT NONE
CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile
TYPE(json_file):: config
!Initialize the json file variable
CALL config%initialize()
!Loads the config file
CALL verboseError('Loading input file...')
CALL config%load(filename = inputFile)
!Reads reference parameters
CALL readReference(config)
!Reads output parameters
CALL readOutput(config)
!Read species
CALL readSpecies(config)
!Read interactions between species
CALL readInteractions(config)
!Read boundaries
CALL readBoundary(config)
!Read Geometry
CALL verboseError('Reading Geometry...')
CALL readGeometry(config)
!Reads case parameters
CALL verboseError('Reading Case Parameters...')
CALL readCase(config)
!Read injection of particles
CALL verboseError('Reading Interactions between species...')
CALL readInject(config)
!Read parallel parameters
CALL verboseError('Reading Parallel configuration...')
CALL readParallel(config)
END SUBROUTINE readConfig
!Reads the reference parameters
SUBROUTINE readReference(config)
USE moduleRefParam
USE moduleConstParam
USE moduleErrors
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
LOGICAL:: found, found_r
CHARACTER(:), ALLOCATABLE:: object
object = 'reference'
!Mandatory parameters that define the case and computes derived parameters
CALL config%get(object // '.density', n_ref, found)
IF (.NOT. found) CALL criticalError('Reference density not found','readReference')
CALL config%get(object // '.mass', m_ref, found)
IF (.NOT. found) CALL criticalError('Reference mass not found','readReference')
CALL config%get(object // '.temperature', T_ref, found)
IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference')
CALL config%get(object // '.radius', r_ref, found_r)
!Derived parameters
v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity
!TODO: Make this solver dependent
IF (found_r) THEN
sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section
L_ref = 1.D0/(sigma_ref*n_ref) !mean free path
ELSE
L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !Debye length
!TODO: Obtain right sigma_ref for PIC case
sigma_ref = PI*(4.D-10)**2 !reference cross section
END IF
ti_ref = L_ref/v_ref !reference time
Vol_ref = L_ref**3 !reference volume
EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field
Volt_ref = EF_ref*L_ref !reference voltage
END SUBROUTINE readReference
!Reads the specific case parameters
SUBROUTINE readCase(config)
USE moduleRefParam
USE moduleErrors
USE moduleCaseParam
USE moduleSolver
USE moduleSpecies
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: object
REAL(8):: time !simulation time in [t]
CHARACTER(:), ALLOCATABLE:: pusherType, EMType, NAType
INTEGER:: nTau, nSolver
INTEGER:: i
CHARACTER(2):: iString
object = 'case'
!Time parameters
CALL config%info(object // '.tau', found, n_children = nTau)
IF (.NOT. found .OR. nTau == 0) CALL criticalError('Required parameter tau not found','readCase')
ALLOCATE(tau(1:nSpecies))
DO i = 1, nTau
WRITE(iString, '(I2)') i
CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(i), found)
END DO
IF (nTau < nSpecies) THEN
CALL warningError('Using minimum time step for some species')
tau(nTau+1:nSpecies) = MINVAL(tau(1:nTau))
END IF
tauMin = MINVAL(tau)
!Gets the simulation time
CALL config%get(object // '.time', time, found)
IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase')
!Convert simulation time to number of iterations
tmax = INT(time/tauMin)
!Gest the pusher for each species
CALL config%info(object // '.pusher', found, n_children = nSolver)
IF (.NOT. found .OR. nSolver /= nSpecies) CALL criticalError('Required parameter pusher not found','readCase')
!Allocates all the pushers for particles
ALLOCATE(solver%pusher(1:nSpecies))
!Initialize pushers
DO i = 1, nSolver
WRITE(iString, '(I2)') i
CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found)
!Associate the type of solver
CALL solver%pusher(i)%init(pusherType, tauMin, tau(i))
END DO
!Gets the solver for the electromagnetic field
CALL config%get(object // '.EMSolver', EMType, found)
CALL solver%initEM(EMType)
SELECT CASE(EMType)
CASE("Electrostatic")
CALL readEMBoundary(config)
END SELECT
!Gest the non-analogue scheme
CALL config%get(object // '.NAScheme', NAType, found)
CALL solver%initNA(NAType)
!Makes tau non-dimensional
tau = tau / ti_ref
tauMin = tauMin / ti_ref
END SUBROUTINE readCase
!Reads configuration for the output files
SUBROUTINE readOutput(config)
USE moduleErrors
USE moduleOutput
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: object
CHARACTER(8) :: date_now=''
CHARACTER(10) :: time_now=''
object = 'output'
CALL config%get(object // '.path', path, found)
CALL config%get(object // '.triggerOutput', triggerOutput, found)
IF (.NOT. found) THEN
triggerOutput = 100
CALL warningError('Using default trigger for output file of 100 iterations')
END IF
!Creates output folder
!TODO: Add option for custon name output_folder
CALL DATE_AND_TIME(date_now, time_now)
folder = date_now(1:4) // '-' // date_now(5:6) // '-' // date_now(7:8) // '_' &
// time_now(1:2) // '.' // time_now(3:4) // '.' // time_now(5:6)
CALL SYSTEM('mkdir ' // path // folder )
CALL config%get(object // '.cpuTime', timeOutput, found)
CALL config%get(object // '.numColl', collOutput, found)
CALL config%get(object // '.EMField', emOutput, found)
CALL config%get(object // '.triggerCPUTime', triggerCPUTime, found)
IF (.NOT. found) THEN
triggerCPUTime = triggerOutput
IF (timeOutput) CALL warningError('No triggerCPUTime found, using same vale as triggerOutput')
END IF
END SUBROUTINE readOutput
!Reads information about the case species
SUBROUTINE readSpecies(config)
USE moduleSpecies
USE moduleErrors
USE moduleRefParam
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(2):: iString
CHARACTER(:), ALLOCATABLE:: object
CHARACTER(:), ALLOCATABLE:: speciesType
REAL(8):: mass, charge
LOGICAL:: found
INTEGER:: i
!Gets the number of species
CALL config%info('species', found, n_children = nSpecies)
!Zero species means critical error
IF (nSpecies == 0) CALL criticalError("No species found", "configRead")
ALLOCATE(species(1:nSpecies))
!Reads information of individual species
DO i = 1, nSpecies
WRITE(iString, '(I2)') i
object = 'species(' // TRIM(iString) // ')'
CALL config%get(object // '.type', speciesType, found)
CALL config%get(object // '.mass', mass, found)
mass = mass/m_ref
!Allocate species depending on type and assign specific parameters
SELECT CASE(speciesType)
CASE ("neutral")
ALLOCATE(species(i)%obj, source=speciesNeutral())
CASE ("charged")
CALL config%get(object // '.charge', charge, found)
ALLOCATE(species(i)%obj, source=speciesCharged(q = charge, &
qm = charge/mass))
CASE DEFAULT
CALL warningError("Species " // speciesType // " not supported yet")
END SELECT
!Assign shared parameters for all species
CALL config%get(object // '.name', species(i)%obj%name, found)
CALL config%get(object // '.weight', species(i)%obj%weight, found)
species(i)%obj%sp = i
species(i)%obj%m = mass
END DO
!Set number of particles to 0 for init state
!TODO: In a future, this should include the particles from init states
nPartOld = 0
!Initialize the lock for the non-analogue (NA) list of particles
CALL OMP_INIT_LOCK(lockNAScheme)
END SUBROUTINE readSpecies
!Reads information about interactions between species
SUBROUTINE readInteractions(config)
USE moduleSpecies
USE moduleCollisions
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(2):: iString, kString
CHARACTER(:), ALLOCATABLE:: object
CHARACTER(:), ALLOCATABLE:: species_i, species_j
CHARACTER(:), ALLOCATABLE:: crossSecFile
CHARACTER(:), ALLOCATABLE:: crossSecFilePath
LOGICAL:: found
INTEGER:: nInteractions, nCollisions
INTEGER:: i, k, ij
INTEGER:: pt_i, pt_j
CALL initInteractionMatrix(interactionMatrix)
CALL config%get('interactions.folderCollisions', pathCollisions, found)
CALL config%info('interactions.collisions', found, n_children = nInteractions)
DO i = 1, nInteractions
WRITE(iString, '(I2)') i
object = 'interactions.collisions(' // TRIM(iString) // ')'
CALL config%get(object // '.species_i', species_i, found)
pt_i = speciesName2Index(species_i)
CALL config%get(object // '.species_j', species_j, found)
pt_j = speciesName2Index(species_j)
CALL config%info(object // '.crossSections', found, n_children = nCollisions)
ij = interactionIndex(pt_i,pt_j)
CALL interactionMatrix(ij)%init(nCollisions)
DO k = 1, nCollisions
WRITE (kString, '(I2)') k
CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile, found)
crossSecFilePath = pathCollisions // crossSecFile
CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m)
END DO
END DO
END SUBROUTINE readInteractions
!Reads boundary conditions for the mesh
SUBROUTINE readBoundary(config)
USE moduleBoundary
USE moduleErrors
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(2):: istring
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
INTEGER:: i
CALL config%info('boundary', found, n_children = nBoundary)
ALLOCATE(boundary(1:nBoundary))
DO i = 1, nBoundary
WRITE(istring, '(i2)') i
object = 'boundary(' // trim(istring) // ')'
ALLOCATE(boundaryGeneric:: boundary(i)%obj)
CALL config%get(object // '.type', boundary(i)%obj%boundaryType, found)
CALL config%get(object // '.physicalSurface', boundary(i)%obj%physicalSurface, found)
boundary(i)%obj%id = i
END DO
END SUBROUTINE readBoundary
!Read the geometry (mesh) for the case
SUBROUTINE readGeometry(config)
USE moduleMesh
USE moduleMeshCylRead, ONLY: meshCylGeneric
USE moduleMesh1DRead, ONLY: mesh1DGeneric
USE moduleErrors
USE moduleOutput
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: geometryType, meshFormat, meshFile
CHARACTER(:), ALLOCATABLE:: fullPath
!Selects the type of geometry.
CALL config%get('geometry.type', geometryType, found)
SELECT CASE(geometryType)
CASE ("2DCyl")
!Creates a 2D cylindrical mesh
ALLOCATE(meshCylGeneric:: mesh)
CASE ("1DCart")
!Creates a 1D cartesian mesh
ALLOCATE(mesh1DGeneric:: mesh)
CASE DEFAULT
CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry")
END SELECT
!Gets the type of mesh
CALL config%get('geometry.meshType', meshFormat, found)
CALL mesh%init(meshFormat)
!Reads the mesh
CALL config%get('geometry.meshFile', meshFile, found)
fullpath = path // meshFile
CALL mesh%readMesh(fullPath)
END SUBROUTINE readGeometry
SUBROUTINE readEMBoundary(config)
USE moduleMesh
USE moduleOutput
USE moduleErrors
USE moduleEM
USE moduleRefParam
USE moduleSpecies
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(2):: istring
INTEGER:: i, e, s
INTEGER:: info
EXTERNAL:: dgetrf
CALL config%info('boundaryEM', found, n_children = nBoundaryEM)
IF (found) ALLOCATE(boundEM(1:nBoundaryEM))
DO i = 1, nBoundaryEM
WRITE(istring, '(i2)') i
object = 'boundaryEM(' // trim(istring) // ')'
CALL config%get(object // '.type', boundEM(i)%typeEM, found)
SELECT CASE(boundEM(i)%typeEM)
CASE ("dirichlet")
CALL config%get(object // '.potential', boundEM(i)%potential, found)
IF (.NOT. found) &
CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary')
boundEM(i)%potential = boundEM(i)%potential/Volt_ref
CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found)
IF (.NOT. found) &
CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary')
CASE DEFAULT
CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary')
END SELECT
END DO
ALLOCATE(qSpecies(1:nSpecies))
DO s = 1, nSpecies
SELECT TYPE(sp => species(s)%obj)
TYPE IS (speciesCharged)
qSpecies(s) = sp%q
CLASS DEFAULT
qSpecies(s) = 0.D0
END SELECT
END DO
!TODO: Improve this
IF (ALLOCATED(boundEM)) THEN
DO e = 1, mesh%numEdges
IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN
DO i = 1, nBoundaryEM
IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN
CALL boundEM(i)%apply(mesh%edges(e)%obj)
END IF
END DO
END IF
END DO
END IF
!Compute the PLU factorization of K once boundary conditions have been read
CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info)
IF (info /= 0) CALL criticalError('Factorization of K matrix failed', 'readEMBoundary')
END SUBROUTINE readEMBoundary
!Reads the injection of particles from the boundaries
SUBROUTINE readInject(config)
USE moduleSpecies
USE moduleErrors
USE moduleInject
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
INTEGER:: i
CHARACTER(2):: istring
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
CHARACTER(:), ALLOCATABLE:: speciesName
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: v
REAL(8), ALLOCATABLE:: T(:), normal(:)
REAL(8):: flow
CHARACTER(:), ALLOCATABLE:: units
INTEGER:: physicalSurface
INTEGER:: sp
CALL config%info('inject', found, n_children = nInject)
ALLOCATE(inject(1:nInject))
nPartInj = 0
DO i = 1, nInject
WRITE(istring, '(i2)') i
object = 'inject(' // trim(istring) // ')'
!Find species
CALL config%get(object // '.species', speciesName, found)
sp = speciesName2Index(speciesName)
CALL config%get(object // '.name', name, found)
CALL config%get(object // '.v', v, found)
CALL config%get(object // '.T', T, found)
CALL config%get(object // '.n', normal, found)
CALL config%get(object // '.flow', flow, found)
CALL config%get(object // '.units', units, found)
CALL config%get(object // '.physicalSurface', physicalSurface, found)
CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
END DO
!Allocate array for injected particles
IF (nPartInj > 0) THEN
ALLOCATE(partInj(1:nPartInj))
END IF
END SUBROUTINE readInject
SUBROUTINE readParallel(config)
USE moduleParallel
USE moduleErrors
USE json_module
IMPLICIT NONE
TYPE(json_file), INTENT(inout):: config
CHARACTER(:), ALLOCATABLE:: object
LOGICAL:: found
!Reads OpenMP parameters
object = 'parallel.OpenMP'
CALL config%get(object // '.nThreads', openMP%nThreads, found)
IF (.NOT. found) THEN
openMP%nthreads = 8
CALL warningError('No OpenMP configuration detected, using 8 threads as default.')
END IF
CALL initParallel()
END SUBROUTINE readParallel
END MODULE moduleInput

91
src/modules/moduleList.f95
View file

@ -1,91 +0,0 @@
!Linked list of particles
MODULE moduleList
USE moduleSpecies
IMPLICIT NONE
TYPE lNode
TYPE(particle), POINTER:: part => NULL()
TYPE(lNode), POINTER:: next => NULL()
END TYPE lNode
TYPE listNode
INTEGER:: amount = 0!TODO: Make private
TYPE(lNode),POINTER:: head => NULL()
TYPE(lNode),POINTER:: tail => NULL()
CONTAINS
PROCEDURE,PASS:: add => addToList
PROCEDURE,PASS:: convert2Array
PROCEDURE,PASS:: erase => eraseList
END TYPE listNode
TYPE(listNode):: partNAScheme !Particles comming from the nonAnalogue scheme
TYPE pointerArray
TYPE(particle), POINTER:: part
END TYPE
CONTAINS
!Adds element to list
SUBROUTINE addToList(self,part)
USE moduleSpecies
CLASS(listNode), INTENT(inout):: self
TYPE(particle),INTENT(in), TARGET:: part
TYPE(lNode),POINTER:: temp
ALLOCATE(temp)
temp%part => part
temp%next => NULL()
self%amount = self%amount + 1
IF (.NOT. ASSOCIATED(self%head)) THEN
!First element
self%head => temp
self%tail => temp
ELSE
!Append element
self%tail%next => temp
self%tail => temp
END IF
END SUBROUTINE addToList
!converts list to array
FUNCTION convert2Array(self) RESULT(partArray)
IMPLICIT NONE
CLASS(listNode), INTENT(in):: self
TYPE(pointerArray), ALLOCATABLE:: partArray(:)
TYPE(lNode), POINTER:: tempNode
INTEGER:: n
ALLOCATE(partArray(1:self%amount))
tempNode => self%head
DO n=1, self%amount
!Point element in array to element in list
partArray(n)%part => tempNode%part
!Go to next element
tempNode => tempNode%next
END DO
END FUNCTION convert2Array
!Erase list
SUBROUTINE eraseList(self)
CLASS(listNode):: self
TYPE(lNode),POINTER:: current, next
current => self%head
DO WHILE (ASSOCIATED(current))
next => current%next
DEALLOCATE(current)
current => next
END DO
IF (ASSOCIATED(self%head)) NULLIFY(self%head)
IF (ASSOCIATED(self%tail)) NULLIFY(self%tail)
self%amount = 0
END SUBROUTINE eraseList
END MODULE moduleList

131
src/modules/moduleOutput.f95
View file

@ -1,131 +0,0 @@
!Contains information about output
MODULE moduleOutput
IMPLICIT NONE
!Output for each node
TYPE outputNode
REAL(8):: den, mom(1:3), tensorS(1:3,1:3)
END TYPE
!Type for EM data in node
TYPE emNode
CHARACTER(:), ALLOCATABLE:: type
REAL(8):: phi
END TYPE emNode
!Output in dimensional units to print
TYPE outputFormat
REAL(8):: density, velocity(1:3), pressure, temperature
END TYPE
CHARACTER(:), ALLOCATABLE:: path
CHARACTER(:), ALLOCATABLE:: folder
INTEGER:: triggerOutput, counterOutput = 0
INTEGER:: triggerCPUTime, counterCPUTime = 0
LOGICAL:: timeOutput = .FALSE.
LOGICAL:: collOutput = .FALSE.
LOGICAL:: emOutput = .FALSE.
CONTAINS
FUNCTION outerProduct(a,b) RESULT(s)
IMPLICIT NONE
REAL(8), DIMENSION(1:3):: a, b
REAL(8):: s(1:3,1:3)
s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3)
END FUNCTION outerProduct
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
SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn)
USE moduleConstParam
USE moduleRefParam
USE moduleSpecies
IMPLICIT NONE
TYPE(outputNode), INTENT(in):: rawValues
TYPE(outputFormat), INTENT(out):: formatValues
REAL(8), INTENT(in):: nodeVol
CLASS(speciesGeneric), INTENT(in):: speciesIn
REAL(8), DIMENSION(1:3,1:3):: tensorTemp
REAL(8), DIMENSION(1:3):: tempVel
REAL(8):: tempVol
!Resets the node outputs
formatValues%density = 0.D0
formatValues%velocity = 0.D0
formatValues%pressure = 0.D0
formatValues%temperature = 0.D0
tempVol = 1.D0/(nodeVol*Vol_ref)
IF (rawValues%den > 0.D0) THEN
tempVel = rawValues%mom(:)/rawValues%den
IF ((tempVel(1) - 1.D0) .EQ. tempVel(1)) THEN
PRINT *, rawValues%mom
END IF
tensorTemp = (rawValues%tensorS(:,:) - rawValues%den*outerProduct(tempVel,tempVel))
formatValues%density = rawValues%den*tempVol
formatValues%velocity(:) = tempVel
IF (tensorTrace(tensorTemp) > 0.D0) THEN
formatValues%pressure = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0
formatValues%temperature = formatValues%pressure/(formatValues%density*kb)
END IF
END IF
formatValues%velocity = formatValues%velocity*v_ref
formatValues%pressure = formatValues%pressure*m_ref*v_ref**2
formatValues%temperature = formatValues%temperature*m_ref*v_ref**2
END SUBROUTINE calculateOutput
SUBROUTINE printTime(t, first)
USE moduleSpecies
USE moduleCompTime
IMPLICIT NONE
INTEGER, INTENT(in):: t
LOGICAL, INTENT(in), OPTIONAL:: first
CHARACTER(:), ALLOCATABLE:: fileName
fileName = 'cpuTime.dat'
IF (timeOutput) THEN
IF (PRESENT(first)) THEN
IF (first) THEN
OPEN(20, file = path // folder // '/' // fileName, action = 'write')
WRITE(20, "(A1, 8X, A1, 9X, A1, 5(A20))") "#","t","n","total","push","reset","collision","weighting"
WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
ELSE
END IF
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
ELSE
OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write')
END IF
WRITE (20, "(I10, I10, 5(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tWeight
CLOSE(20)
END IF
END SUBROUTINE printTime
END MODULE moduleOutput

20
src/modules/moduleParallel.f95
View file

@ -1,20 +0,0 @@
MODULE moduleParallel
IMPLICIT NONE
TYPE openMP_type
INTEGER:: nThreads
END TYPE
TYPE(openMP_type):: openMP
CONTAINS
SUBROUTINE initParallel()
IMPLICIT NONE
!Starts threads for OpenMP parallelization
CALL OMP_SET_NUM_THREADS(openMP%nThreads)
END SUBROUTINE initParallel
END MODULE moduleParallel

9
src/modules/moduleRefParam.f95
View file

@ -1,9 +0,0 @@
!Reference parameters
MODULE moduleRefParam
!Parameters that define the problem (inputs)
REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref
!Reference parameters for non-dimensional problem
REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref
END MODULE moduleRefParam

584
src/modules/moduleSolver.f95
View file

@ -1,584 +0,0 @@
MODULE moduleSolver
!Generic type for pusher of particles
TYPE, PUBLIC:: pusherGeneric
PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL()
!Boolean to indicate if the species is moved in the iteration
LOGICAL:: pushSpecies
!How many interations between advancing the species
INTEGER:: every
CONTAINS
PROCEDURE, PASS:: init => initPusher
END TYPE pusherGeneric
!Generic type for solver
TYPE, PUBLIC:: solverGeneric
TYPE(pusherGeneric), ALLOCATABLE:: pusher(:)
PROCEDURE(solveEM_interface), POINTER, NOPASS:: solveEM => NULL()
PROCEDURE(nonAnalogue_interface), POINTER, NOPASS:: nonAnalogue => NULL()
CONTAINS
PROCEDURE, PASS:: initEM
PROCEDURE, PASS:: initNA
PROCEDURE, PASS:: updateParticleCell
PROCEDURE, PASS:: updatePushSpecies
END TYPE solverGeneric
INTERFACE
!Push a particle
PURE SUBROUTINE push_interafece(part)
USE moduleSpecies
TYPE(particle), INTENT(inout):: part
END SUBROUTINE push_interafece
!Solve the electromagnetic field
SUBROUTINE solveEM_interface()
END SUBROUTINE solveEM_interface
!Apply nonAnalogue scheme to a particle
SUBROUTINE nonAnalogue_interface(part, volOld, volNew)
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
CLASS(meshVol), POINTER, INTENT(in):: volOld
CLASS(meshVol), POINTER, INTENT(inout):: volNew
END SUBROUTINE nonAnalogue_interface
END INTERFACE
TYPE(solverGeneric):: solver
CONTAINS
!Init Pusher
SUBROUTINE initPusher(self, pusherType, tau, tauSp)
USE moduleErrors
IMPLICIT NONE
CLASS(pusherGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE:: pusherType
REAL(8):: tau, tauSp
SELECT CASE(pusherType)
CASE('2DCylNeutral')
self%pushParticle => pushCylNeutral
CASE('2DCylCharged')
self%pushParticle => pushCylCharged
CASE('1DCartCharged')
self%pushParticle => push1DCharged
CASE DEFAULT
CALL criticalError('Solver ' // pusherType // ' not found','readCase')
END SELECT
self%pushSpecies = .FALSE.
self%every = INT(tauSp/tau)
END SUBROUTINE initPusher
SUBROUTINE initEM(self, EMType)
IMPLICIT NONE
CLASS(solverGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE:: EMType
SELECT CASE(EMType)
CASE('Electrostatic')
self%solveEM => solveElecField
END SELECT
END SUBROUTINE initEM
!Initialize the non-analogue scheme
SUBROUTINE initNA(self, NAType)
USE moduleMesh
IMPLICIT NONE
CLASS(solverGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE:: NAType
SELECT CASE(NAType)
CASE ('Volume')
self%nonAnalogue => volumeNAScheme
END SELECT
END SUBROUTINE initNA
!Do all pushes for particles
SUBROUTINE doPushes()
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
INTEGER:: n
INTEGER:: sp
!$OMP DO
DO n=1, nPartOld
!Select species type
sp = partOld(n)%sp
!Checks if the species sp is update this iteration
IF (solver%pusher(sp)%pushSpecies) THEN
!Push particle
CALL solver%pusher(sp)%pushParticle(partOld(n))
!Find cell in wich particle reside
CALL solver%updateParticleCell(partOld(n))
END IF
END DO
!$OMP END DO
END SUBROUTINE doPushes
!Push one particle. Boris pusher for 2D Cyl Neutral particle
PURE SUBROUTINE pushCylNeutral(part)
USE moduleSpecies
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
TYPE(particle):: part_temp
REAL(8):: tauSp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: v_p_oh_star(2:3)
part_temp = part
!Time step for the species
tauSp = tau(part_temp%sp)
!z
part_temp%v(1) = part%v(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp
!r,theta
v_p_oh_star(2) = part%v(2)
x_new = part%r(2) + v_p_oh_star(2)*tauSp
v_p_oh_star(3) = part%v(3)
y_new = v_p_oh_star(3)*tauSp
r = DSQRT(x_new**2+y_new**2)
part_temp%r(2) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE pushCylNeutral
!Push one particle. Boris pusher for 2D Cyl Charged particle
PURE SUBROUTINE pushCylCharged(part)
USE moduleSpecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
REAL(8):: v_p_oh_star(2:3)
TYPE(particle):: part_temp
REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha
REAL(8):: tauSp
REAL(8):: qmEFt(1:3)!charge*tauSp*EF/mass
part_temp = part
!Time step for the species
tauSp = tau(part_temp%sp)
!Get electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp
!z
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp
!r,theta
v_p_oh_star(2) = part%v(2) + qmEFt(2)
x_new = part%r(2) + v_p_oh_star(2)*tauSp
v_p_oh_star(3) = part%v(3) + qmEFt(3)
y_new = v_p_oh_star(3)*tauSp
r = DSQRT(x_new**2+y_new**2)
part_temp%r(2) = r
IF (r > 0.D0) THEN
sin_alpha = y_new/r
cos_alpha = x_new/r
ELSE
sin_alpha = 0.D0
cos_alpha = 1.D0
END IF
part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3)
part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3)
part_temp%n_in = .FALSE. !Assume particle is outside until cell is found
!Copy temporal particle to particle
part=part_temp
END SUBROUTINE pushCylCharged
!Push charged particles in 1D cartesian coordinates
PURE SUBROUTINE push1DCharged(part)
USE moduleSPecies
USE moduleEM
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
TYPE(particle):: part_temp
REAL(8):: tauSp
REAL(8):: qmEFt(1:3)
part_temp = part
!Time step for particle species
tauSp = tau(part_temp%sp)
!Get the electric field at particle position
qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp
!x
part_temp%v(1) = part%v(1) + qmEFt(1)
part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp
part_temp%n_in = .FALSE.
part = part_temp
END SUBROUTINE push1DCharged
!Do the collisions in all the cells
SUBROUTINE doCollisions()
USE moduleMesh
IMPLICIT NONE
INTEGER:: e
!$OMP DO SCHEDULE(DYNAMIC)
DO e=1, mesh%numVols
CALL mesh%vols(e)%obj%collision()
END DO
!$OMP END DO
END SUBROUTINE doCollisions
SUBROUTINE doReset()
USE moduleSpecies
USE moduleList
IMPLICIT NONE
INTEGER:: nn, n
INTEGER, SAVE:: nPartNew
INTEGER, SAVE:: nInjIn, nOldIn, nNAScheme
TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:)
TYPE(lNode), POINTER:: partCurr, partNext
!$OMP SECTIONS
!$OMP SECTION
nInjIn = 0
IF (ALLOCATED(partInj)) THEN
nInjIn = COUNT(partInj%n_in)
END IF
!$OMP SECTION
nOldIn = 0
IF (ALLOCATED(partOld)) THEN
nOldIn = COUNT(partOld%n_in)
END IF
!$OMP SECTION
nNAScheme = partNAScheme%amount
!$OMP END SECTIONS
!$OMP BARRIER
!$OMP SINGLE
CALL MOVE_ALLOC(partOld, partTemp)
nPartNew = nInjIn + nOldIn + nNAScheme
ALLOCATE(partOld(1:nPartNew))
!$OMP END SINGLE
!$OMP SECTIONS
!$OMP SECTION
nn = 0
DO n = 1, nPartInj
IF (partInj(n)%n_in) THEN
nn = nn + 1
partOld(nn) = partInj(n)
END IF
END DO
!$OMP SECTION
nn = nInjIn
DO n = 1, nPartOld
IF (partTemp(n)%n_in) THEN
nn = nn + 1
partOld(nn) = partTemp(n)
END IF
END DO
!$OMP SECTION
nn = nInjIn + nOldIn
partCurr => partNAScheme%head
DO n = 1, nNAScheme
partNext => partCurr%next
partOld(nn+n) = partCurr%part
DEALLOCATE(partCurr)
partCurr => partNext
END DO
IF (ASSOCIATED(partNAScheme%head)) NULLIFY(partNAScheme%head)
IF (ASSOCIATED(partNAScheme%tail)) NULLIFY(partNAScheme%tail)
partNAScheme%amount = 0
!$OMP END SECTIONS
!$OMP SINGLE
nPartOld = nPartNew
!$OMP END SINGLE
END SUBROUTINE doReset
!Scatter particles in the grid
SUBROUTINE doScatter
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
INTEGER:: n
!Loops over the particles to scatter them
!$OMP DO
DO n=1, nPartOld
CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n))
END DO
!$OMP END DO
END SUBROUTINE doScatter
SUBROUTINE doEMField()
IMPLICIT NONE
IF (ASSOCIATED(solver%solveEM)) THEN
CALL solver%solveEM()
END IF
END SUBROUTINE doEMField
!Solving the Poisson equation for electrostatic potential
SUBROUTINE solveElecField()
USE moduleEM
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
INTEGER, SAVE:: INFO
INTEGER:: n
REAL(8), ALLOCATABLE, SAVE:: tempF(:)
EXTERNAL:: dgetrs
!$OMP SINGLE
ALLOCATE(tempF(1:mesh%numNodes))
!$OMP END SINGLE
CALL assembleSourceVector(tempF)
!$OMP SINGLE
CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, &
mesh%IPIV, tempF, mesh%numNodes, info)
!$OMP END SINGLE
IF (info == 0) THEN
!Suscessful resolution of Poission equation
!$OMP DO
DO n = 1, mesh%numNodes
mesh%nodes(n)%obj%emData%phi = tempF(n)
END DO
!$OMP END DO
ELSE
!$OMP SINGLE
CALL criticalError('Poisson equation failed', 'solveElecField')
!$OMP END SINGLE
END IF
!$OMP SINGLE
DEALLOCATE(tempF)
!$OMP END SINGLE
END SUBROUTINE solveElecField
!Modify particle weight as a function of cell volume and splits particle
SUBROUTINE volumeNAScheme(part, volOld, volNew)
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
CLASS(meshVol), POINTER, INTENT(in):: volOld
CLASS(meshVol), POINTER, INTENT(inout):: volNew
REAL(8):: fractionVolume, fractionWeight
INTEGER:: nSplit
!If particle has change cell, call nonAnalogue scheme
IF (volOld%n /= volNew%n) THEN
fractionVolume = volOld%volume/volNew%volume
part%weight = part%weight * fractionVolume
fractionWeight = part%weight / species(part%sp)%obj%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.
END IF
END IF
END SUBROUTINE volumeNAScheme
!Subroutine to split the particle 'part' into a number 'nSplit' of particles.
!'nSplit-1' particles are added to the partNAScheme list
SUBROUTINE splitParticle(part, nSplit, vol)
USE moduleSpecies
USE moduleList
USE moduleMesh
USE OMP_LIB
IMPLICIT NONE
TYPE(particle), INTENT(inout):: part
INTEGER, INTENT(in):: nSplit
CLASS(meshVol), INTENT(inout):: vol
REAL(8):: newWeight
TYPE(particle), POINTER:: newPart
INTEGER:: p
newWeight = part%weight / nSplit
!Assign new weight to original particle
part%weight = newWeight
!Add new particles to list of NA particles
DO p = 2, nSplit
!Allocate the pointer for the new particles
ALLOCATE(newPart)
!Copy data from original particle
newPart = part
CALL OMP_SET_LOCK(lockNAScheme)
CALL partNAScheme%add(newPart)
CALL OMP_UNSET_LOCK(lockNASCheme)
!Add particle to cell list
CALL OMP_SET_lock(vol%lock)
CALL vol%listPart_in%add(newPart)
CALL OMP_UNSET_lock(vol%lock)
END DO
END SUBROUTINE splitParticle
SUBROUTINE updateParticleCell(self, part)
USE moduleSpecies
USE moduleMesh
IMPLICIT NONE
CLASS(solverGeneric), INTENT(in):: self
TYPE(particle), INTENT(inout):: part
CLASS(meshVol), POINTER:: volOld, volNew
volOld => mesh%vols(part%vol)%obj
CALL volOld%findCell(part)
volNew => mesh%vols(part%vol)%obj
!Call the NA shcme
IF (ASSOCIATED(self%nonAnalogue)) THEN
CALL self%nonAnalogue(part, volOld, volNew)
END IF
END SUBROUTINE updateParticleCell
!Update the information about if a species needs to be moved this iteration
SUBROUTINE updatePushSpecies(self, t)
USE moduleSpecies
IMPLICIT NONE
CLASS(solverGeneric), INTENT(inout):: self
INTEGER, INTENT(in):: t
INTEGER:: s
DO s=1, nSpecies
self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0
END DO
END SUBROUTINE updatePushSpecies
!Output the different data and information
SUBROUTINE doOutput(t)
USE moduleMesh
USE moduleOutput
USE moduleSpecies
USE moduleCompTime
IMPLICIT NONE
INTEGER, INTENT(in):: t
counterOutput = counterOutput + 1
IF (counterOutput >= triggerOutput .OR. &
t == tmax .OR. t == 0) THEN
!Resets output counter
counterOutput=0
CALL mesh%printOutput(t)
CALL mesh%printColl(t)
CALL mesh%printEM(t)
WRITE(*, "(5X,A21,I10,A1,I10)") "t/tmax: ", t, "/", tmax
WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld
IF (t == 0) THEN
WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms"
ELSE
WRITE(*, "(5X,A21,F8.1,A2)") "iteration time: ", 1.D3*tStep, "ms"
END IF
IF (nPartOld > 0) THEN
WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns"
END IF
WRITE(*,*)
END IF
counterCPUTime = counterCPUTime + 1
IF (counterCPUTime >= triggerCPUTime .OR. &
t == tmax .OR. t == 0) THEN
!Reset CPU Time counter
counterCPUTime = 0
CALL printTime(t, t == 0)
END IF
END SUBROUTINE doOutput
END MODULE moduleSolver

71
src/modules/moduleSpecies.f95
View file

@ -1,71 +0,0 @@
!Contains the information about species (particles)
MODULE moduleSpecies
USE moduleCaseParam
USE OMP_LIB
IMPLICIT NONE
TYPE, ABSTRACT:: speciesGeneric
CHARACTER(:), ALLOCATABLE:: name
REAL(8):: m=0.D0, weight=0.D0
INTEGER:: sp=0
END TYPE speciesGeneric
TYPE, EXTENDS(speciesGeneric):: speciesNeutral
END TYPE speciesNeutral
TYPE, EXTENDS(speciesGeneric):: speciesCharged
REAL(8):: q=0.D0, qm=0.D0
END TYPE speciesCharged
TYPE:: speciesCont
CLASS(speciesGeneric), ALLOCATABLE:: obj
END TYPE
INTEGER:: nSpecies
TYPE(speciesCont), ALLOCATABLE:: species(:)
TYPE particle
REAL(8):: r(1:3) !Position
REAL(8):: v(1:3) !Velocity
INTEGER:: sp !Particle species id
INTEGER:: vol !Index of element in which the particle is located
REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p.
LOGICAL:: n_in !Flag that indicates if a particle is in the domain
REAL(8):: weight=0.D0 !weight of particle
REAL(8):: qm = 0.D0 !charge over mass
END TYPE particle
!Number of old particles
INTEGER:: nPartOld
INTEGER:: nPartInj
!Arrays that contain the particles
TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration
TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles
INTEGER(KIND=OMP_LOCK_KIND):: lockNAScheme !Lock for the NA list of particles
CONTAINS
FUNCTION speciesName2Index(speciesName) RESULT(sp)
USE moduleErrors
IMPLICIT NONE
CHARACTER(:), ALLOCATABLE:: speciesName
INTEGER:: sp
INTEGER:: n
sp = 0
DO n = 1, nSpecies
IF (speciesName == species(n)%obj%name) THEN
sp = species(n)%obj%sp
EXIT
END IF
END DO
!If no species is found, call a critical error
IF (sp == 0) CALL criticalError('Species ' // speciesName // ' not found.', 'speciesName2Index')
END FUNCTION speciesName2Index
END MODULE moduleSpecies

121
src/modules/moduleTable.f95
View file

@ -1,121 +0,0 @@
MODULE moduleTable
TYPE:: table1D
REAL(8):: xMin, xMax
REAL(8):: fMin, fMax
REAL(8), ALLOCATABLE, DIMENSION(:):: x, f, k
CONTAINS
PROCEDURE, PASS:: init => initTable1D
PROCEDURE, PASS:: get => getValueTable1D
PROCEDURE, PASS:: convert => convertUnits
END TYPE table1D
CONTAINS
SUBROUTINE initTable1D(self, tableFile)
USE moduleErrors
IMPLICIT NONE
CLASS(table1D), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(IN):: tableFile
CHARACTER(100):: dummy
INTEGER:: amount
INTEGER:: i
INTEGER:: stat
INTEGER:: id = 20
OPEN(id, file = tableFile)
amount = 0
DO
READ(id, '(A)', iostat = stat) dummy
!Skip comment
IF (INDEX(dummy,'#') /= 0) CYCLE
!If EOF or error, exit file
IF (stat /= 0) EXIT
!Add row
amount = amount + 1
END DO
IF (amount == 0) CALL criticalError('Table ' // tableFile // ' is empty', 'initTable1D')
IF (amount == 1) CALL criticalError('Table ' // tableFile // ' has only one row', 'initTable1D')
!Go bback to initial point
REWIND(id)
!Allocate table arrays
ALLOCATE(self%x(1:amount), self%f(1:amount), self%k(1:amount))
self%x = 0.D0
self%f = 0.D0
self%k = 0.D0
i = 0
DO
READ(id, '(A)', iostat = stat) dummy
IF (INDEX(dummy,'#') /= 0) CYCLE
IF (stat /= 0) EXIT
!Add data
!TODO: substitute with extracting information from dummy
BACKSPACE(id)
i = i + 1
READ(id, *) self%x(i), self%f(i)
END DO
CLOSE(10)
self%xMin = self%x(1)
self%xMax = self%x(amount)
self%fMin = self%f(1)
self%fMax = self%f(amount)
DO i = 1, amount - 1
self%k(i) = (self%f(i+1) - self%f(i))/(self%x(i+1) - self%x(i))
END DO
END SUBROUTINE initTable1D
FUNCTION getValueTable1D(self, x) RESULT(f)
IMPLICIT NONE
CLASS(table1D), INTENT(in):: self
REAL(8):: x
REAL(8):: f
REAL(8):: deltaX
INTEGER:: i
IF (x <= self%xMin) THEN
f = self%fMin
ELSEIF (x >= self%xMax) THEN
f = self%fMax
ELSE
i = MINLOC(x - self%x, 1)
deltaX = x - self%x(i)
IF (deltaX < 0 ) THEN
i = i - 1
deltaX = x - self%x(i)
END IF
f = self%f(i) + self%k(i)*deltaX
END IF
END FUNCTION getValueTable1D
SUBROUTINE convertUnits(self, data_x, data_f)
IMPLICIT NONE
CLASS(table1D), INTENT(inout):: self
REAL(8):: data_x, data_f
self%x = self%x * data_x
self%f = self%f * data_f
self%k = self%k * data_f / data_x
END SUBROUTINE convertUnits
END MODULE moduleTable