diff --git a/src/fpakc.f90 b/src/fpakc.f90 new file mode 100644 index 0000000..a24b492 --- /dev/null +++ b/src/fpakc.f90 @@ -0,0 +1,103 @@ +! 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 + diff --git a/src/modules/mesh/1DCart/moduleMesh1D.f90 b/src/modules/mesh/1DCart/moduleMesh1D.f90 new file mode 100644 index 0000000..63acc27 --- /dev/null +++ b/src/modules/mesh/1DCart/moduleMesh1D.f90 @@ -0,0 +1,489 @@ +!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:2) + + 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 + diff --git a/src/modules/mesh/1DCart/moduleMesh1DBoundary.f90 b/src/modules/mesh/1DCart/moduleMesh1DBoundary.f90 new file mode 100644 index 0000000..bd4e6d6 --- /dev/null +++ b/src/modules/mesh/1DCart/moduleMesh1DBoundary.f90 @@ -0,0 +1,40 @@ +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 diff --git a/src/modules/mesh/1DCart/moduleMesh1DRead.f90 b/src/modules/mesh/1DCart/moduleMesh1DRead.f90 new file mode 100644 index 0000000..1c52795 --- /dev/null +++ b/src/modules/mesh/1DCart/moduleMesh1DRead.f90 @@ -0,0 +1,268 @@ +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 diff --git a/src/modules/mesh/2DCyl/makefile b/src/modules/mesh/2DCyl/makefile new file mode 100644 index 0000000..73cd926 --- /dev/null +++ b/src/modules/mesh/2DCyl/makefile @@ -0,0 +1,11 @@ +all : moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.o + +moduleMeshCyl.o: moduleMeshCyl.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +moduleMeshCylRead.o: moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + diff --git a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 new file mode 100644 index 0000000..dd3a283 --- /dev/null +++ b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 @@ -0,0 +1,1022 @@ +!moduleMeshCyl: 2D axial symmetric extension of generic mesh from GMSH format. +! x == z +! y == r +! z == theta (unused) +MODULE moduleMeshCyl + USE moduleMesh + IMPLICIT NONE + + !Values for Gauss integral + REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) + REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) + + REAL(8), PARAMETER:: xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /) + REAL(8), PARAMETER:: xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /) + REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /) + + TYPE, PUBLIC, EXTENDS(meshNode):: meshNodeCyl + !Element coordinates + REAL(8):: r = 0.D0, z = 0.D0 + CONTAINS + PROCEDURE, PASS:: init => initNodeCyl + PROCEDURE, PASS:: getCoordinates => getCoordCyl + + END TYPE meshNodeCyl + + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdgeCyl + !Element coordinates + REAL(8):: r(1:2) = 0.D0, z(1:2) = 0.D0 + !Connectivity to nodes + CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() + CONTAINS + PROCEDURE, PASS:: init => initEdgeCyl + PROCEDURE, PASS:: getNodes => getNodesCyl + PROCEDURE, PASS:: randPos => randPosCyl + + END TYPE meshEdgeCyl + + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl + CONTAINS + PROCEDURE, PASS:: detJac => detJCyl + PROCEDURE, PASS:: invJac => invJCyl + PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi + PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi + PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer + + END TYPE meshVolCyl + + 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, dz, dr) + IMPORT meshVolCyl + CLASS(meshVolCyl), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + + END SUBROUTINE partialDer_interface + + END INTERFACE + + TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylQuad + !Element coordinates + REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 + !Connectivity to nodes + CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL() + !Connectivity to adjacent elements + CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() + REAL(8):: arNodes(1:4) = 0.D0 + + CONTAINS + PROCEDURE, PASS:: init => initVolQuadCyl + PROCEDURE, PASS:: area => areaQuad + PROCEDURE, NOPASS:: fPsi => fPsiQuad + PROCEDURE, NOPASS:: dPsi => dPsiQuad + PROCEDURE, NOPASS:: dPsiXi1 => dPsiQuadXi1 + PROCEDURE, NOPASS:: dPsiXi2 => dPsiQuadXi2 + PROCEDURE, PASS:: partialDer => partialDerQuad + PROCEDURE, PASS:: elemK => elemKQuad + PROCEDURE, PASS:: elemF => elemFQuad + PROCEDURE, NOPASS:: weight => weightQuad + PROCEDURE, NOPASS:: inside => insideQuad + PROCEDURE, PASS:: scatter => scatterQuad + PROCEDURE, PASS:: gatherEF => gatherEFQuad + PROCEDURE, PASS:: getNodes => getNodesQuad + PROCEDURE, PASS:: phy2log => phy2logQuad + PROCEDURE, PASS:: nextElement => nextElementQuad + PROCEDURE, PASS:: resetOutput => resetOutputQuad + + END TYPE meshVolCylQuad + + TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylTria + !Element coordinates + REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0 + !Connectivity to nodes + CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() + !Connectivity to adjacent elements + CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL() + REAL(8):: arNodes(1:3) = 0.D0 + !Derivatives in z,r real space + REAL(8), DIMENSION(1:3):: dPsiZ, dPsiR + + CONTAINS + PROCEDURE, PASS:: init => initVolTriaCyl + PROCEDURE, PASS:: area => areaTria + PROCEDURE, NOPASS:: fPsi => fPsiTria + PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1 + PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2 + PROCEDURE, PASS:: partialDer => partialDerTria + PROCEDURE, PASS:: elemK => elemKTria + PROCEDURE, PASS:: elemF => elemFTria + PROCEDURE, NOPASS:: weight => weightTria + PROCEDURE, NOPASS:: inside => insideTria + PROCEDURE, PASS:: scatter => scatterTria + PROCEDURE, PASS:: gatherEF => gatherEFTria + PROCEDURE, PASS:: getNodes => getNodesTria + PROCEDURE, PASS:: phy2log => phy2logTria + PROCEDURE, PASS:: nextElement => nextElementTria + PROCEDURE, PASS:: resetOutput => resetOutputTria + + END TYPE meshVolCylTria + + CONTAINS + + !NODE FUNCTIONS + !Inits node element + SUBROUTINE initNodeCyl(self, n, r) + USE moduleSpecies + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshNodeCyl), INTENT(out):: self + INTEGER, INTENT(in):: n + REAL(8), INTENT(in):: r(1:3) + + self%n = n + self%z = r(2)/L_ref + self%r = r(1)/L_ref + !Node volume, to be determined in mesh + self%v = 0.D0 + + !Allocates output: + ALLOCATE(self%output(1:nSpecies)) + + END SUBROUTINE initNodeCyl + + PURE FUNCTION getCoordCyl(self) RESULT(r) + IMPLICIT NONE + + CLASS(meshNodeCyl), INTENT(in):: self + REAL(8):: r(1:3) + + r = (/self%z, self%r, 0.D0/) + + END FUNCTION getCoordCyl + + !EDGE FUNCTIONS + !Inits edge element + SUBROUTINE initEdgeCyl(self, n, p, bt, physicalSurface) + IMPLICIT NONE + + CLASS(meshEdgeCyl), 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, 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%z = (/r1(1), r2(1)/) + self%r = (/r1(2), r2(2)/) + !Normal vector + self%normal = (/ self%r(2)-self%r(1), & + self%z(2)-self%z(1), & + 0.D0 /) + !Boundary index + self%bt = bt + !Phyiscal Surface + self%physicalSurface = physicalSurface + + END SUBROUTINE initEdgeCyl + + !Get nodes from edge + PURE FUNCTION getNodesCyl(self) RESULT(n) + IMPLICIT NONE + + CLASS(meshEdgeCyl), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + ALLOCATE(n(1:2)) + n = (/self%n1%n, self%n2%n /) + + END FUNCTION getNodesCyl + + !Calculates a random position in edge + FUNCTION randPosCyl(self) RESULT(r) + CLASS(meshEdgeCyl), INTENT(in):: self + REAL(8):: rnd + REAL(8):: r(1:3) + REAL(8):: p1(1:2), p2(1:2) + + CALL RANDOM_NUMBER(rnd) + + p1 = (/self%z(1), self%r(1) /) + p2 = (/self%z(2), self%r(2) /) + r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 + r(3) = 0.D0 + + END FUNCTION randPosCyl + + !VOLUME FUNCTIONS + !QUAD FUNCTIONS + !Inits quadrilateral element + SUBROUTINE initVolQuadCyl(self, n, p) + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(out):: self + INTEGER, INTENT(in):: n + INTEGER, INTENT(in):: p(:) + REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 + + self%n = n + self%n1 => mesh%nodes(p(1))%obj + self%n2 => mesh%nodes(p(2))%obj + self%n3 => mesh%nodes(p(3))%obj + self%n4 => mesh%nodes(p(4))%obj + !Get element coordinates + r1 = self%n1%getCoordinates() + r2 = self%n2%getCoordinates() + r3 = self%n3%getCoordinates() + r4 = self%n4%getCoordinates() + self%z = (/r1(1), r2(1), r3(1), r4(1)/) + self%r = (/r1(2), r2(2), r3(2), r4(2)/) + + !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%n3%v = self%n3%v + self%arNodes(3) + self%n4%v = self%n4%v + self%arNodes(4) + + self%sigmaVrelMax = sigma_ref/L_ref**2 + + CALL OMP_INIT_LOCK(self%lock) + + END SUBROUTINE initVolQuadCyl + + !Computes element area + PURE SUBROUTINE areaQuad(self) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(inout):: self + REAL(8):: r, xi(1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:4) + + self%volume = 0.D0 + self%arNodes = 0.D0 + !2D 1 point Gauss Quad Integral + xi = 0.D0 + detJ = self%detJac(xi)*8.D0*PI !4*2*pi + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + self%volume = r*detJ + self%arNodes = fPsi*r*detJ + + END SUBROUTINE areaQuad + + !Computes element functions in point xi + PURE FUNCTION fPsiQuad(xi) RESULT(fPsi) + IMPLICIT NONE + + REAL(8),INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + ALLOCATE(fPsi(1:4)) + + fPsi(1) = (1.D0-xi(1))*(1.D0-xi(2)) + fPsi(2) = (1.D0+xi(1))*(1.D0-xi(2)) + fPsi(3) = (1.D0+xi(1))*(1.D0+xi(2)) + fPsi(4) = (1.D0-xi(1))*(1.D0+xi(2)) + fPsi = fPsi*0.25D0 + + END FUNCTION fPsiQuad + + !Derivative element function at coordinates xi + PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: dPsi(:,:) + + ALLOCATE(dPsi(1:2,1:4)) + + dPsi(1,:) = dPsiQuadXi1(xi(2)) + dPsi(2,:) = dPsiQuadXi2(xi(1)) + + END FUNCTION dPsiQuad + + !Derivative element function (xi1) + PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1) + IMPLICIT NONE + + REAL(8),INTENT(in):: xi2 + REAL(8):: dPsiXi1(1:4) + + dPsiXi1(1) = -(1.D0-xi2) + dPsiXi1(2) = (1.D0-xi2) + dPsiXi1(3) = (1.D0+xi2) + dPsiXi1(4) = -(1.D0+xi2) + dPsiXi1 = dPsiXi1*0.25D0 + + END FUNCTION dPsiQuadXi1 + + !Derivative element function (xi2) + PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2) + IMPLICIT NONE + + REAL(8),INTENT(in):: xi1 + REAL(8):: dPsiXi2(1:4) + + dPsiXi2(1) = -(1.D0-xi1) + dPsiXi2(2) = -(1.D0+xi1) + dPsiXi2(3) = (1.D0+xi1) + dPsiXi2(4) = (1.D0-xi1) + dPsiXi2 = dPsiXi2*0.25D0 + + END FUNCTION dPsiQuadXi2 + + PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + + dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) + + END SUBROUTINE partialDerQuad + + !Computes element local stiffness matrix + PURE FUNCTION elemKQuad(self) RESULT(ke) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8):: r, xi(1:3) + REAL(8):: fPsi(1:4), dPsi(1:2,1:4) + REAL(8):: ke(1:4,1:4) + REAL(8):: invJ(1:2,1:2), detJ + INTEGER:: l, m + + ke=0.D0 + xi=0.D0 + !Start 2D Gauss Quad Integral + DO l=1, 3 + xi(2) = corQuad(l) + dPsi(1,:) = self%dPsiXi1(xi(2)) + DO m = 1, 3 + xi(1) = corQuad(m) + dPsi(2,:) = self%dPsiXi2(xi(1)) + fPsi = self%fPsi(xi) + detJ = self%detJac(xi,dPsi) + invJ = self%invJac(xi,dPsi) + r = DOT_PRODUCT(fPsi,self%r) + ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wQuad(l)*wQuad(m)/detJ + + END DO + END DO + ke = ke*2.D0*PI + + END FUNCTION elemKQuad + + !Computes the local source vector for a force f + PURE FUNCTION elemFQuad(self, source) RESULT(localF) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: source(1:) + REAL(8), ALLOCATABLE:: localF(:) + REAL(8):: r, xi(1:3) + REAL(8):: fPsi(1:4) + REAL(8):: detJ, f + INTEGER:: l, m + + ALLOCATE(localF(1:4)) + localF = 0.D0 + xi = 0.D0 + DO l=1, 3 + xi(1) = corQuad(l) + DO m = 1, 3 + xi(2) = corQuad(m) + detJ = self%detJac(xi) + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + f = DOT_PRODUCT(fPsi,source) + localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ + + END DO + END DO + localF = localF*2.D0*PI + + END FUNCTION elemFQuad + + !Computes weights in the element nodes + PURE FUNCTION weightQuad(xi) RESULT(w) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: w(1:4) + + w = fPsiQuad(xi) + + END FUNCTION weightQuad + + !Checks if a particle is inside a quad element + PURE FUNCTION insideQuad(xi) RESULT(ins) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + LOGICAL:: ins + + ins = (xi(1) >= -1.D0 .AND. xi(1) <= 1.D0) .AND. & + (xi(2) >= -1.D0 .AND. xi(2) <= 1.D0) + + END FUNCTION insideQuad + + !Scatter properties of particle into element nodes + SUBROUTINE scatterQuad(self, part) + USE moduleOutput + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(particle), INTENT(in):: part + TYPE(outputNode), POINTER:: vertex + REAL(8):: w_p(1:4) + 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 + + vertex => self%n3%output(part%sp) + vertex%den = vertex%den + part%weight*w_p(3) + vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS + + vertex => self%n4%output(part%sp) + vertex%den = vertex%den + part%weight*w_p(4) + vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS + + END SUBROUTINE scatterQuad + + !Gathers the electric field at position xi + PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: dPsi(1:2,1:4) + REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates + REAL(8):: invJ(1:2,1:2), detJ + REAL(8):: phi(1:4) + REAL(8):: EF(1:3) + + phi = (/self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi, & + self%n4%emData%phi /) + + dPsi = self%dPsi(xi) + detJ = self%detJac(xi,dPsi) + invJ = self%invJac(xi,dPsi) + dPsiR = MATMUL(invJ, dPsi)/detJ + EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) + EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) + EF(3) = 0.D0 + + END FUNCTION gatherEFQuad + + !Gets nodes from quadrilateral element + PURE FUNCTION getNodesQuad(self) RESULT(n) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + ALLOCATE(n(1:4)) + n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) + + END FUNCTION getNodesQuad + + !Transforms physical coordinates to element coordinates + PURE FUNCTION phy2logQuad(self,r) RESULT(xN) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: r(1:3) + REAL(8):: xN(1:3) + REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2) + REAL(8):: dPsi(1:2,1:4), fPsi(1:4) + REAL(8):: conv + + !Iterative newton method to transform coordinates + conv=1.D0 + xO=0.D0 + + DO WHILE(conv>1.D-4) + dPsi = self%dPsi(xO) + invJ = self%invJac(xO, dPsi) + fPsi = self%fPsi(xO) + f(1) = DOT_PRODUCT(fPsi,self%z)-r(1) + f(2) = DOT_PRODUCT(fPsi,self%r)-r(2) + detJ = self%detJac(xO,dPsi) + xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ + conv=MAXVAL(DABS(xN-xO),1) + xO=xN + + END DO + + END FUNCTION phy2logQuad + + !Gets the next element for a logical position xi + SUBROUTINE nextElementQuad(self, xi, nextElement) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + CLASS(*), POINTER, INTENT(out):: nextElement + REAL(8):: xiArray(1:4) + INTEGER:: nextInt + + xiArray = (/ -xi(2), xi(1), xi(2), -xi(1) /) + nextInt = MAXLOC(xiArray,1) + !Selects the higher value of directions and searches in that direction + NULLIFY(nextElement) + SELECT CASE (nextInt) + CASE (1) + nextElement => self%e1 + CASE (2) + nextElement => self%e2 + CASE (3) + nextElement => self%e3 + CASE (4) + nextElement => self%e4 + END SELECT + + END SUBROUTINE nextElementQuad + + !Reset the output of nodes in quad element + PURE SUBROUTINE resetOutputQuad(self) + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(meshVolCylQuad), 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 + + self%n3%output(k)%den = 0.D0 + self%n3%output(k)%mom = 0.D0 + self%n3%output(k)%tensorS = 0.D0 + + self%n4%output(k)%den = 0.D0 + self%n4%output(k)%mom = 0.D0 + self%n4%output(k)%tensorS = 0.D0 + + END DO + + END SUBROUTINE resetOutputQuad + + !TRIA ELEMENT + !Init tria element + SUBROUTINE initVolTriaCyl(self, n, p) + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(out):: self + INTEGER, INTENT(in):: n + INTEGER, INTENT(in):: p(:) + REAL(8), DIMENSION(1:3):: r1, r2, r3 + REAL(8):: A + + !Assign node index + self%n = n + + !Assign nodes to element + self%n1 => mesh%nodes(p(1))%obj + self%n2 => mesh%nodes(p(2))%obj + self%n3 => mesh%nodes(p(3))%obj + !Get element coordinates + r1 = self%n1%getCoordinates() + r2 = self%n2%getCoordinates() + r3 = self%n3%getCoordinates() + self%z = (/r1(1), r2(1), r3(1)/) + self%r = (/r1(2), r2(2), r3(2)/) + !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%n3%v = self%n3%v + self%arNodes(3) + + !Derivatives in z/r for shape functions (node independent) + !TODO: This is used because invJ.dPsi does not produce the right Electric field + A = self%z(2)*self%r(3) - self%z(3)*self%r(2) + & + self%z(3)*self%r(1) - self%z(1)*self%r(3) + & + self%z(1)*self%r(2) - self%z(2)*self%r(1) + + self%dPsiZ = (/ self%r(2)-self%r(3), & + self%r(3)-self%r(1), & + self%r(1)-self%r(2) /)/A + + self%dPsiR = (/ self%z(3)-self%z(2), & + self%z(1)-self%z(3), & + self%z(2)-self%z(1) /)/A + + self%sigmaVrelMax = sigma_ref/L_ref**2 + + CALL OMP_INIT_LOCK(self%lock) + + END SUBROUTINE initVolTriaCyl + + !Calculates area for triangular element + PURE SUBROUTINE areaTria(self) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(inout):: self + REAL(8):: r, xi(1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:3) + + self%volume = 0.D0 + self%arNodes = 0.D0 + !2D 1 point Gauss Quad Integral + xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) + detJ = self%detJac(xi)*PI !2PI*1/2 + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + self%volume = r*detJ + self%arNodes = fPsi*r*detJ + + END SUBROUTINE areaTria + + !Shape functions for triangular element + PURE FUNCTION fPsiTria(xi) RESULT(fPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + ALLOCATE(fPsi(1:3)) + + fPsi(1) = 1.D0 - xi(1) - xi(2) + fPsi(2) = xi(1) + fPsi(3) = xi(2) + + END FUNCTION fPsiTria + + !Derivative element function at coordinates xi + PURE FUNCTION dPsiTria(xi) RESULT(dPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: dPsi(:,:) + + ALLOCATE(dPsi(1:2,1:3)) + + dPsi(1,:) = dPsiTriaXi1(xi(2)) + dPsi(2,:) = dPsiTriaXi2(xi(1)) + + END FUNCTION dPsiTria + + !Derivative element function (xi1) + PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi2 + REAL(8):: dPsiXi1(1:3) + + dPsiXi1(1) = -1.D0 + dPsiXi1(2) = 1.D0 + dPsiXi1(3) = 0.D0 + + END FUNCTION dPsiTriaXi1 + + !Derivative element function (xi1) + PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi1 + REAL(8):: dPsiXi2(1:3) + + dPsiXi2(1) = -1.D0 + dPsiXi2(2) = 0.D0 + dPsiXi2(3) = 1.D0 + + END FUNCTION dPsiTriaXi2 + + PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + + dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) + + END SUBROUTINE partialDerTria + + !Computes element local stiffness matrix + PURE FUNCTION elemKTria(self) RESULT(ke) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8):: r, xi(1:3) + REAL(8):: fPsi(1:3), dPsi(1:2,1:3) + REAL(8):: ke(1:3,1:3) + REAL(8):: invJ(1:2,1:2), detJ + INTEGER:: l + + ke=0.D0 + xi=0.D0 + !Start 2D Gauss Quad Integral + DO l=1, 4 + xi(1) = xi1Tria(l) + xi(2) = xi2Tria(l) + dPsi = self%dPsi(xi) + detJ = self%detJac(xi,dPsi) + invJ = self%invJac(xi,dPsi) + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ + + END DO + ke = ke*2.D0*PI + + END FUNCTION elemKTria + + !Computes element local source vector + PURE FUNCTION elemFTria(self, source) RESULT(localF) + USE moduleConstParam + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: source(1:) + REAL(8), ALLOCATABLE:: localF(:) + REAL(8):: fPsi(1:3) + REAL(8):: r, xi(1:3) + REAL(8):: detJ, f + INTEGER:: l + + ALLOCATE(localF(1:3)) + localF = 0.D0 + xi = 0.D0 + !Start 2D Gauss Quad Integral + DO l=1, 4 + xi(1) = xi1Tria(l) + xi(2) = xi2Tria(l) + detJ = self%detJac(xi) + fPsi = self%fPsi(xi) + r = DOT_PRODUCT(fPsi,self%r) + f = DOT_PRODUCT(fPsi,source) + localF = localF + r*f*fPsi*wTria(l)*detJ + + END DO + localF = localF*2.D0*PI + + END FUNCTION elemFTria + + !Computes weights in the element nodes + PURE FUNCTION weightTria(xi) RESULT(w) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + REAL(8), ALLOCATABLE:: w(:) + + w = fPsiTria(xi) + + END FUNCTION weightTria + + PURE FUNCTION insideTria(xi) RESULT(ins) + IMPLICIT NONE + + REAL(8), INTENT(in):: xi(1:3) + LOGICAL:: ins + + ins = xi(1) >= 0.D0 .AND. & + xi(2) >= 0.D0 .AND. & + 1.D0 - xi(1) - xi(2) >= 0.D0 + + END FUNCTION insideTria + + !Scatter properties of particles into element + SUBROUTINE scatterTria(self, part) + USE moduleOutput + USE moduleSpecies + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + CLASS(particle), INTENT(in):: part + TYPE(outputNode), POINTER:: vertex + REAL(8):: w_p(1:3) + 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 + + vertex => self%n3%output(part%sp) + vertex%den = vertex%den + part%weight*w_p(3) + vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) + vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS + + END SUBROUTINE scatterTria + + !Gathers the electric field at position xi + PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: phi(1:3) + REAL(8):: EF(1:3) + + phi = (/self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi/) + + EF(1) = -DOT_PRODUCT(self%dPsiZ(:), phi) + EF(2) = -DOT_PRODUCT(self%dPsiR(:), phi) + EF(3) = 0.D0 + + END FUNCTION gatherEFTria + + !Gets node indexes from triangular element + PURE FUNCTION getNodesTria(self) RESULT(n) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + ALLOCATE(n(1:3)) + n = (/self%n1%n, self%n2%n, self%n3%n /) + + END FUNCTION getNodesTria + + !Transforms physical coordinates to element coordinates + PURE FUNCTION phy2logTria(self,r) RESULT(xi) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: r(1:3) + REAL(8):: xi(1:3) + REAL(8):: invJ(1:2,1:2), detJ + REAL(8):: deltaR(1:2) + REAL(8):: dPsi(1:2,1:3) + + !Direct method to convert coordinates + xi = 0.D0 !Irrelevant, required for input + deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /) + dPsi = self%dPsi(xi) + invJ = self%invJac(xi, dPsi) + detJ = self%detJac(xi, dPsi) + xi(1:2) = MATMUL(invJ,deltaR)/detJ + + END FUNCTION phy2logTria + + SUBROUTINE nextElementTria(self, xi, nextElement) + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + CLASS(*), POINTER, INTENT(out):: nextElement + REAL(8):: xiArray(1:3) + INTEGER:: nextInt + + xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /) + nextInt = MINLOC(xiArray,1) + NULLIFY(nextElement) + SELECT CASE (nextInt) + CASE (1) + nextElement => self%e1 + CASE (2) + nextElement => self%e2 + CASE (3) + nextElement => self%e3 + END SELECT + + END SUBROUTINE nextElementTria + + !Reset the output of nodes in tria element + PURE SUBROUTINE resetOutputTria(self) + USE moduleSpecies + USE moduleOutput + IMPLICIT NONE + + CLASS(meshVolCylTria), 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 + + self%n3%output(k)%den = 0.D0 + self%n3%output(k)%mom = 0.D0 + self%n3%output(k)%tensorS = 0.D0 + + END DO + + END SUBROUTINE resetOutputTria + + + !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS + !Computes element Jacobian determinant + PURE FUNCTION detJCyl(self, xi, dPsi_in) RESULT(dJ) + IMPLICIT NONE + + CLASS(meshVolCyl), 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):: dz(1:2), dr(1:2) + + IF(PRESENT(dPsi_in)) THEN + dPsi = dPsi_in + + ELSE + dPsi = self%dPsi(xi) + + END IF + + CALL self%partialDer(dPsi, dz, dr) + dJ = dz(1)*dr(2)-dz(2)*dr(1) + + END FUNCTION detJCyl + + !Computes element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJCyl(self,xi,dPsi_in) RESULT(invJ) + IMPLICIT NONE + + CLASS(meshVolCyl), 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):: dz(1:2), dr(1:2) + REAL(8):: invJ(1:2,1:2) + + IF(PRESENT(dPsi_in)) THEN + dPsi=dPsi_in + + ELSE + dPsi = self%dPsi(xi) + + END IF + + CALL self%partialDer(dPsi, dz, dr) + invJ(1,:) = (/ dr(2), -dz(2) /) + invJ(2,:) = (/ -dr(1), dz(1) /) + + END FUNCTION invJCyl + +END MODULE moduleMeshCyl diff --git a/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 b/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 new file mode 100644 index 0000000..7c29115 --- /dev/null +++ b/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 @@ -0,0 +1,80 @@ +!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 diff --git a/src/modules/mesh/2DCyl/moduleMeshCylRead.f90 b/src/modules/mesh/2DCyl/moduleMeshCylRead.f90 new file mode 100644 index 0000000..2da8c31 --- /dev/null +++ b/src/modules/mesh/2DCyl/moduleMeshCylRead.f90 @@ -0,0 +1,610 @@ +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 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 new file mode 100644 index 0000000..285c0a1 --- /dev/null +++ b/src/modules/mesh/moduleMesh.f90 @@ -0,0 +1,609 @@ +!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 + REAL(8):: rndReal + 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 + CALL RANDOM_NUMBER(rndReal) + rnd = 1 + FLOOR(nPart*rndReal) + part_i => partTemp(rnd)%part + CALL RANDOM_NUMBER(rndReal) + rnd = 1 + FLOOR(nPart*rndReal) + 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 diff --git a/src/modules/moduleBoundary.f90 b/src/modules/moduleBoundary.f90 new file mode 100644 index 0000000..71caa16 --- /dev/null +++ b/src/modules/moduleBoundary.f90 @@ -0,0 +1,36 @@ +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 diff --git a/src/modules/moduleCaseParam.f90 b/src/modules/moduleCaseParam.f90 new file mode 100644 index 0000000..5f1023a --- /dev/null +++ b/src/modules/moduleCaseParam.f90 @@ -0,0 +1,9 @@ +!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 + diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 new file mode 100644 index 0000000..c95c90e --- /dev/null +++ b/src/modules/moduleCollisions.f90 @@ -0,0 +1,175 @@ +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 + REAL(8):: rnd + + !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 + CALL RANDOM_NUMBER(rnd) + IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) 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 + CALL RANDOM_NUMBER(rnd) + alpha = PI*rnd + part_i%v(1) = v_i*DCOS(alpha) + part_i%v(2) = v_i*DSIN(alpha) + CALL RANDOM_NUMBER(rnd) + alpha = PI*rnd + part_j%v(1) = v_j*DCOS(alpha) + part_j%v(2) = v_j*DSIN(alpha) + + END IF + + + END SUBROUTINE + +END MODULE moduleCollisions diff --git a/src/modules/moduleCompTime.f90 b/src/modules/moduleCompTime.f90 new file mode 100644 index 0000000..28273ce --- /dev/null +++ b/src/modules/moduleCompTime.f90 @@ -0,0 +1,11 @@ +!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 + diff --git a/src/modules/moduleConstParam.f90 b/src/modules/moduleConstParam.f90 new file mode 100644 index 0000000..634f845 --- /dev/null +++ b/src/modules/moduleConstParam.f90 @@ -0,0 +1,15 @@ +!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 + diff --git a/src/modules/moduleEM.f90 b/src/modules/moduleEM.f90 new file mode 100644 index 0000000..876093e --- /dev/null +++ b/src/modules/moduleEM.f90 @@ -0,0 +1,128 @@ +!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 diff --git a/src/modules/moduleErrors.f90 b/src/modules/moduleErrors.f90 new file mode 100644 index 0000000..0a95db8 --- /dev/null +++ b/src/modules/moduleErrors.f90 @@ -0,0 +1,43 @@ +!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 diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 new file mode 100644 index 0000000..eb3c737 --- /dev/null +++ b/src/modules/moduleInject.f90 @@ -0,0 +1,202 @@ +!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 + 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) + CALL RANDOM_NUMBER(x) + END DO + CALL RANDOM_NUMBER(y) + + 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 + REAL(8):: rnd + 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 + CALL RANDOM_NUMBER(rnd) + randomX = INT(DBLE(self%nEdges-1)*rnd) + 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 diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 new file mode 100644 index 0000000..d31b5cf --- /dev/null +++ b/src/modules/moduleInput.f90 @@ -0,0 +1,559 @@ +! 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 EXECUTE_COMMAND_LINE('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 diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 new file mode 100644 index 0000000..bea9ba0 --- /dev/null +++ b/src/modules/moduleList.f90 @@ -0,0 +1,91 @@ +!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 diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90 new file mode 100644 index 0000000..33d9a94 --- /dev/null +++ b/src/modules/moduleOutput.f90 @@ -0,0 +1,131 @@ +!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 + diff --git a/src/modules/moduleParallel.f90 b/src/modules/moduleParallel.f90 new file mode 100644 index 0000000..e945e6d --- /dev/null +++ b/src/modules/moduleParallel.f90 @@ -0,0 +1,22 @@ +MODULE moduleParallel + IMPLICIT NONE + + TYPE openMP_type + INTEGER:: nThreads + + END TYPE + + TYPE(openMP_type):: openMP + + CONTAINS + SUBROUTINE initParallel() + IMPLICIT NONE + + EXTERNAL:: OMP_SET_NUM_THREADS + + !Starts threads for OpenMP parallelization + CALL OMP_SET_NUM_THREADS(openMP%nThreads) + + END SUBROUTINE initParallel + +END MODULE moduleParallel diff --git a/src/modules/moduleRefParam.f90 b/src/modules/moduleRefParam.f90 new file mode 100644 index 0000000..9d2b8d9 --- /dev/null +++ b/src/modules/moduleRefParam.f90 @@ -0,0 +1,9 @@ +!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 + diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 new file mode 100644 index 0000000..60c1bd3 --- /dev/null +++ b/src/modules/moduleSolver.f90 @@ -0,0 +1,584 @@ +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 + diff --git a/src/modules/moduleSpecies.f90 b/src/modules/moduleSpecies.f90 new file mode 100644 index 0000000..da5128c --- /dev/null +++ b/src/modules/moduleSpecies.f90 @@ -0,0 +1,71 @@ +!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 diff --git a/src/modules/moduleTable.f90 b/src/modules/moduleTable.f90 new file mode 100644 index 0000000..1ea89e5 --- /dev/null +++ b/src/modules/moduleTable.f90 @@ -0,0 +1,121 @@ +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