From ac277259408bd6670ea17085527e3f91881fb04e Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 12 Jul 2024 23:08:19 +0200 Subject: [PATCH] Big one... I should've commited before, but I wanted to make things compile. The big change is that I've added a global time step so the parameter does not need to be passed in each function. This is useful as we are moving towards using time profiles for boundary conditions and injection of particles (not in this branch, but in the future and the procedure will be quite similar) --- src/fpakc.f90 | 25 ++- src/modules/common/moduleCaseParam.f90 | 5 + src/modules/init/moduleInput.f90 | 134 +++++++------ .../mesh/inout/0D/moduleMeshOutput0D.f90 | 21 +- .../inout/gmsh2/moduleMeshOutputGmsh2.f90 | 38 ++-- .../mesh/inout/moduleMeshInoutCommon.f90 | 8 +- .../mesh/inout/vtu/moduleMeshInputVTU.f90 | 2 +- .../mesh/inout/vtu/moduleMeshOutputVTU.f90 | 35 ++-- src/modules/mesh/moduleMesh.f90 | 15 +- src/modules/moduleInject.f90 | 26 +-- src/modules/moduleProbe.f90 | 28 +-- src/modules/output/moduleOutput.f90 | 6 +- .../solver/electromagnetic/moduleEM.f90 | 180 +++++++++++++----- src/modules/solver/moduleSolver.f90 | 37 ++-- 14 files changed, 340 insertions(+), 220 deletions(-) diff --git a/src/fpakc.f90 b/src/fpakc.f90 index e90a5e0..ae6e7bb 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -11,12 +11,12 @@ PROGRAM fpakc USE OMP_LIB IMPLICIT NONE - ! t = time step - INTEGER:: t ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 ! inputFile = path+name of input file CHARACTER(:), ALLOCATABLE:: inputFile + ! generic integer for time step + INTEGER:: t tStep = omp_get_wtime() !Gets the input file @@ -32,10 +32,13 @@ PROGRAM fpakc CALL initOutput(inputFile) !Do '0' iteration - t = tInitial + timeStep = tInitial !$OMP PARALLEL DEFAULT(SHARED) !$OMP SINGLE + ! Initial reset of probes + CALL resetProbes() + CALL verboseError("Initial scatter of particles...") !$OMP END SINGLE CALL doScatter() @@ -49,19 +52,21 @@ PROGRAM fpakc tStep = omp_get_wtime() - tStep !Output initial state - CALL doOutput(t) + CALL doOutput() CALL verboseError('Starting main loop...') !$OMP PARALLEL DEFAULT(SHARED) DO t = tInitial + 1, tFinal - !Insert new particles and push them !$OMP SINGLE tStep = omp_get_wtime() + ! Update global time step index + timeStep = t + !Checks if a species needs to me moved in this iteration - CALL solver%updatePushSpecies(t) + CALL solver%updatePushSpecies() !Checks if probes need to be calculated this iteration - CALL resetProbes(t) + CALL resetProbes() tPush = omp_get_wtime() !$OMP END SINGLE @@ -79,7 +84,7 @@ PROGRAM fpakc !$OMP END SINGLE IF (doMCCollisions) THEN - CALL meshForMCC%doCollisions(t) + CALL meshForMCC%doCollisions() END IF @@ -124,12 +129,12 @@ PROGRAM fpakc !$OMP SINGLE tEMField = omp_get_wtime() - tEMField - CALL doAverage(t) + CALL doAverage() tStep = omp_get_wtime() - tStep !Output data - CALL doOutput(t) + CALL doOutput() !$OMP END SINGLE END DO diff --git a/src/modules/common/moduleCaseParam.f90 b/src/modules/common/moduleCaseParam.f90 index 8df3210..551d867 100644 --- a/src/modules/common/moduleCaseParam.f90 +++ b/src/modules/common/moduleCaseParam.f90 @@ -2,8 +2,13 @@ MODULE moduleCaseParam !Final and initial iterations INTEGER:: tFinal, tInitial = 0 + ! Global index of current iteration + INTEGER:: timeStep + ! Time step for all species REAL(8), ALLOCATABLE:: tau(:) + ! Minimum time step REAL(8):: tauMin + ! Time step for Monte-Carlo Collisions REAL(8):: tauColl END MODULE moduleCaseParam diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 4216c73..6800ac1 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -264,8 +264,8 @@ MODULE moduleInput CALL readEMBoundary(config) !Read constant magnetic field DO i = 1, 3 - WRITE(istring, '(i2)') i - CALL config%get(object // '.B(' // istring // ')', B(i), found) + WRITE(iString, '(i2)') i + CALL config%get(object // '.B(' // iString // ')', B(i), found) IF (.NOT. found) THEN CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver') @@ -799,7 +799,7 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config INTEGER:: i, s - CHARACTER(2):: istring, sString + CHARACTER(2):: iString, sString CHARACTER(:), ALLOCATABLE:: object, bType REAL(8):: Tw, cw !Wall temperature and specific heat !Neutral Properties @@ -815,8 +815,8 @@ MODULE moduleInput CALL config%info('boundary', found, n_children = nBoundary) ALLOCATE(boundary(1:nBoundary)) DO i = 1, nBoundary - WRITE(istring, '(i2)') i - object = 'boundary(' // TRIM(istring) // ')' + WRITE(iString, '(i2)') i + object = 'boundary(' // TRIM(iString) // ')' boundary(i)%n = i CALL config%get(object // '.name', boundary(i)%name, found) @@ -1100,13 +1100,13 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(2):: istring + CHARACTER(2):: iString INTEGER:: i CHARACTER(:), ALLOCATABLE:: speciesName REAL(8), ALLOCATABLE, DIMENSION(:):: r REAL(8), ALLOCATABLE, DIMENSION(:):: v1, v2, v3 INTEGER, ALLOCATABLE, DIMENSION(:):: points - REAL(8):: timeStep + REAL(8):: everyTimeStep CALL config%info('output.probes', found, n_children = nProbes) @@ -1114,7 +1114,7 @@ MODULE moduleInput DO i = 1, nProbes WRITE(iString, '(I2)') i - object = 'output.probes(' // trim(istring) // ')' + object = 'output.probes(' // trim(iString) // ')' CALL config%get(object // '.species', speciesName, found) CALL config%get(object // '.position', r, found) @@ -1122,16 +1122,14 @@ MODULE moduleInput CALL config%get(object // '.velocity_2', v2, found) CALL config%get(object // '.velocity_3', v3, found) CALL config%get(object // '.points', points, found) - CALL config%get(object // '.timeStep', timeStep, found) + CALL config%get(object // '.timeStep', everyTimeStep, found) IF (ANY(points < 2)) CALL criticalError("Number of points in probe " // iString // " incorrect", 'readProbes') - CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, timeStep) + CALL probe(i)%init(i, speciesName, r, v1, v2, v3, points, everyTimeStep) END DO - CALL resetProbes(tInitial) - END SUBROUTINE readProbes SUBROUTINE readEMBoundary(config) @@ -1147,40 +1145,56 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found - CHARACTER(2):: istring - INTEGER:: i, e, s + CHARACTER(:), ALLOCATABLE:: typeEM + REAL(8):: potential + INTEGER:: physicalSurface + CHARACTER(:), ALLOCATABLE:: timeProfile + INTEGER:: b, e, s, n, ni + CHARACTER(2):: bString INTEGER:: info EXTERNAL:: dgetrf CALL config%info('boundaryEM', found, n_children = nBoundaryEM) - IF (found) ALLOCATE(boundEM(1:nBoundaryEM)) + IF (found) ALLOCATE(boundaryEM(1:nBoundaryEM)) - DO i = 1, nBoundaryEM - WRITE(istring, '(I2)') i - object = 'boundaryEM(' // trim(istring) // ')' + DO b = 1, nBoundaryEM + WRITE(bString, '(I2)') b + object = 'boundaryEM(' // trim(bString) // ')' - CALL config%get(object // '.type', boundEM(i)%typeEM, found) + CALL config%get(object // '.type', typeEM, found) - SELECT CASE(boundEM(i)%typeEM) + SELECT CASE(typeEM) CASE ("dirichlet") - CALL config%get(object // '.potential', boundEM(i)%potential, found) - IF (.NOT. found) & + CALL config%get(object // '.potential', potential, found) + IF (.NOT. found) THEN 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) & + END IF + + CALL config%get(object // '.physicalSurface', physicalSurface, found) + IF (.NOT. found) THEN CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') + END IF + + CALL initDirichlet(boundaryEM(b)%obj, physicalSurface, potential) + CASE ("dirichletTime") - CALL config%get(object // '.potential', boundEM(i)%potential, found) - IF (.NOT. found) & + CALL config%get(object // '.potential', potential, found) + IF (.NOT. found) THEN CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') - boundEM(i)%potential = boundEM(i)%potential/Volt_ref + + END IF + + CALL config%get(object // '.physicalSurface', physicalSurface, found) + IF (.NOT. found) THEN + CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') + + END IF CASE DEFAULT - CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary') + CALL criticalError('Boundary type ' // typeEM // ' not yet supported', 'readEMBoundary') END SELECT @@ -1199,18 +1213,28 @@ MODULE moduleInput END DO - 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) + ! Modify K matrix due to boundary conditions + DO b = 1, nBoundaryEM + SELECT TYPE(boundary => boundaryEM(b)%obj) + TYPE IS(boundaryEMDirichlet) + DO n = 1, boundary%nNodes + ni = boundary%nodes(n)%obj%n + mesh%K(ni, :) = 0.D0 + mesh%K(ni, ni) = 1.D0 - END IF - END DO - END IF - END DO - END IF + END DO + + TYPE IS(boundaryEMDirichletTime) + DO n = 1, boundary%nNodes + ni = boundary%nodes(n)%obj%n + mesh%K(ni, :) = 0.D0 + mesh%K(ni, ni) = 1.D0 + + END DO + + END SELECT + + END DO !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) @@ -1231,13 +1255,13 @@ MODULE moduleInput TYPE(json_file), INTENT(inout):: config INTEGER:: i - CHARACTER(2):: istring + CHARACTER(2):: iString CHARACTER(:), ALLOCATABLE:: object LOGICAL:: found CHARACTER(:), ALLOCATABLE:: speciesName CHARACTER(:), ALLOCATABLE:: name REAL(8):: v - REAL(8), ALLOCATABLE:: T(:), normal(:) + REAL(8), ALLOCATABLE:: temperature(:), normal(:) REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface @@ -1248,8 +1272,8 @@ MODULE moduleInput ALLOCATE(inject(1:nInject)) nPartInj = 0 DO i = 1, nInject - WRITE(istring, '(i2)') i - object = 'inject(' // trim(istring) // ')' + WRITE(iString, '(i2)') i + object = 'inject(' // trim(iString) // ')' !Find species CALL config%get(object // '.species', speciesName, found) @@ -1257,7 +1281,7 @@ MODULE moduleInput 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 // '.T', temperature, found) CALL config%get(object // '.n', normal, found) IF (.NOT. found) THEN ALLOCATE(normal(1:3)) @@ -1269,7 +1293,7 @@ MODULE moduleInput particlesPerEdge = 0 CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) - CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge) + CALL inject(i)%init(i, v, normal, temperature, flow, units, sp, physicalSurface, particlesPerEdge) CALL readVelDistr(config, inject(i), object) @@ -1329,28 +1353,28 @@ MODULE moduleInput TYPE(injectGeneric), INTENT(inout):: inj CHARACTER(:), ALLOCATABLE, INTENT(in):: object INTEGER:: i - CHARACTER(2):: istring + CHARACTER(2):: iString CHARACTER(:), ALLOCATABLE:: fvType LOGICAL:: found - REAL(8):: v, T, m + REAL(8):: v, temperature, m !Reads species mass m = inj%species%m !Reads distribution functions for velocity DO i = 1, 3 - WRITE(istring, '(i2)') i - CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found) - IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // & + WRITE(iString, '(i2)') i + CALL config%get(object // '.velDist('// TRIM(iString) //')', fvType, found) + IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // iString // & " found for " // object, 'readVelDistr') SELECT CASE(fvType) CASE ("Maxwellian") - T = inj%T(i) - CALL initVelDistMaxwellian(inj%v(i)%obj, t, m) + temperature = inj%temperature(i) + CALL initVelDistMaxwellian(inj%v(i)%obj, temperature, m) CASE ("Half-Maxwellian") - T = inj%T(i) - CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m) + temperature = inj%temperature(i) + CALL initVelDistHalfMaxwellian(inj%v(i)%obj, temperature, m) CASE ("Delta") v = inj%vMod*inj%n(i) diff --git a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 index 97ec729..c0dcfbb 100644 --- a/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 +++ b/src/modules/mesh/inout/0D/moduleMeshOutput0D.f90 @@ -1,22 +1,22 @@ MODULE moduleMeshOutput0D CONTAINS - SUBROUTINE printOutput0D(self, t) + SUBROUTINE printOutput0D(self) USE moduleMesh USE moduleRefParam USE moduleSpecies USE moduleOutput + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: i TYPE(outputFormat):: output CHARACTER(:), ALLOCATABLE:: fileName DO i = 1, nSpecies fileName='OUTPUT_' // species(i)%obj%name // '.dat' - IF (t == 0) THEN + IF (timeStep == 0) THEN OPEN(20, file = path // folder // '/' // fileName, action = 'write') WRITE(20, "(A1, 14X, A5, A20, 40X, A20, 2(A20))") "#","t (s)","density (m^-3)", "velocity (m/s)", & "pressure (Pa)", "temperature (K)" @@ -27,14 +27,17 @@ MODULE moduleMeshOutput0D OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') CALL calculateOutput(self%nodes(1)%obj%output(i), output, self%nodes(1)%obj%v, species(i)%obj) - WRITE(20, "(7(ES20.6E3))") REAL(t)*tauMin*ti_ref, output%density, output%velocity, output%pressure, output%temperature + WRITE(20, "(7(ES20.6E3))") REAL(timeStep)*tauMin*ti_ref, output%density, & + output%velocity, & + output%pressure, & + output%temperature CLOSE(20) END DO END SUBROUTINE printOutput0D - SUBROUTINE printColl0D(self, t) + SUBROUTINE printColl0D(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -43,12 +46,11 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t CHARACTER(:), ALLOCATABLE:: fileName INTEGER:: k fileName='OUTPUT_Collisions.dat' - IF (t == tInitial) THEN + IF (timeStep == tInitial) THEN OPEN(20, file = path // folder // '/' // fileName, action = 'write') WRITE(20, "(A1, 14X, A5, A20)") "#","t (s)","collisions" WRITE(*, "(6X,A15,A)") "Creating file: ", fileName @@ -57,12 +59,12 @@ MODULE moduleMeshOutput0D END IF OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE(20, "(ES20.6E3, 10I20)") REAL(t)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) + WRITE(20, "(ES20.6E3, 10I20)") REAL(timeStep)*tauMin*ti_ref, (self%cells(1)%obj%tallyColl(k)%tally, k=1,nCollPairs) CLOSE(20) END SUBROUTINE printColl0D - SUBROUTINE printEM0D(self, t) + SUBROUTINE printEM0D(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -70,7 +72,6 @@ MODULE moduleMeshOutput0D IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printEM0D diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index ccdcf3d..8176bb5 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -80,50 +80,50 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE writeGmsh2FooterElementData !Prints the scattered properties of particles into the nodes - SUBROUTINE printOutputGmsh2(self, t) + SUBROUTINE printOutputGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleSpecies USE moduleOutput USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: n, i TYPE(outputFormat):: output(1:self%numNodes) REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'msh', t) + fileName = formatFileName(prefix, species(i)%obj%name, 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) CALL writeGmsh2HeaderMesh(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' density (m^-3)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) WRITE(60, "(I6,ES20.6E3)") n, output(n)%density END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', t, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' velocity (m s^-1)', timeStep, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Pressure (Pa)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure END DO CALL writeGmsh2FooterNodeData(60) - CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(60, species(i)%obj%name // ' Temperature (K)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature END DO @@ -135,7 +135,7 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printOutputGmsh2 !Prints the number of collisions into the volumes - SUBROUTINE printCollGmsh2(self, t) + SUBROUTINE printCollGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -145,7 +145,6 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: numEdges INTEGER:: k, c INTEGER:: n @@ -167,9 +166,9 @@ MODULE moduleMeshOutputGmsh2 END SELECT IF (collOutput) THEN - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref - fileName = formatFileName(prefix, 'Collisions', 'msh', t) + fileName = formatFileName(prefix, 'Collisions', 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) @@ -179,7 +178,7 @@ MODULE moduleMeshOutputGmsh2 DO c = 1, interactionMatrix(k)%amount WRITE(cString, "(I2)") c title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString - CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numCells) + CALL writeGmsh2HeaderElementData(60, title, timeStep, time, 1, self%numCells) DO n=1, self%numCells WRITE(60, "(I6,I10)") n + numEdges, self%cells(n)%obj%tallyColl(k)%tally(c) END DO @@ -196,7 +195,7 @@ MODULE moduleMeshOutputGmsh2 END SUBROUTINE printCollGmsh2 !Prints the electrostatic EM properties into the nodes and volumes - SUBROUTINE printEMGmsh2(self, t) + SUBROUTINE printEMGmsh2(self) USE moduleMesh USE moduleRefParam USE moduleCaseParam @@ -205,7 +204,6 @@ MODULE moduleMeshOutputGmsh2 IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: n, e REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName @@ -214,27 +212,27 @@ MODULE moduleMeshOutputGmsh2 Xi = (/ 0.D0, 0.D0, 0.D0 /) IF (emOutput) THEN - time = DBLE(t)*tauMin*ti_ref + time = DBLE(timeStep)*tauMin*ti_ref - fileName = formatFileName(prefix, 'EMField', 'msh', t) + fileName = formatFileName(prefix, 'EMField', 'msh', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (20, file = path // folder // '/' // fileName) CALL writeGmsh2HeaderMesh(20) - CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', t, time, 1, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Potential (V)', timeStep, time, 1, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref END DO CALL writeGmsh2FooterNodeData(20) - CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', t, time, 3, self%numCells) + CALL writeGmsh2HeaderElementData(20, 'Electric Field (V m^-1)', timeStep, time, 3, self%numCells) DO e=1, self%numCells WRITE(20, *) e+self%numEdges, self%cells(e)%obj%gatherElectricField(Xi)*EF_ref END DO CALL writeGmsh2FooterElementData(20) - CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', t, time, 3, self%numNodes) + CALL writeGmsh2HeaderNodeData(20, 'Magnetic Field (T)', timeStep, time, 3, self%numNodes) DO n=1, self%numNodes WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref END DO diff --git a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 index e4a6c72..7dc2e84 100644 --- a/src/modules/mesh/inout/moduleMeshInoutCommon.f90 +++ b/src/modules/mesh/inout/moduleMeshInoutCommon.f90 @@ -3,17 +3,17 @@ MODULE moduleMeshInoutCommon CHARACTER(LEN=4):: prefix = 'Step' CONTAINS - PURE FUNCTION formatFileName(prefix, suffix, extension, t) RESULT(fileName) + PURE FUNCTION formatFileName(prefix, suffix, extension, timeStep) RESULT(fileName) USE moduleOutput IMPLICIT NONE CHARACTER(*), INTENT(in):: prefix, suffix, extension - INTEGER, INTENT(in), OPTIONAL:: t + INTEGER, INTENT(in), OPTIONAL:: timeStep CHARACTER (LEN=iterationDigits):: tString CHARACTER(:), ALLOCATABLE:: fileName - IF (PRESENT(t)) THEN - WRITE(tString, iterationFormat) t + IF (PRESENT(timeStep)) THEN + WRITE(tString, iterationFormat) timeStep fileName = prefix // '_' // tString // '_' // suffix // '.' // extension ELSE diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index 184d645..e07db01 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -167,7 +167,7 @@ MODULE moduleMeshInputVTU CLASS(meshGeneric), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8):: r(1:3) !3 generic coordinates - INTEGER:: fileID, error, found + INTEGER:: fileID CHARACTER(LEN=256):: line INTEGER:: numNodes, numElements, numEdges INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 index da90b6b..81e4bbf 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 @@ -215,17 +215,16 @@ MODULE moduleMeshOutputVTU END SUBROUTINE writeEM - SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection) + SUBROUTINE writeCollection(fileID, fileNameStep, fileNameCollection) USE moduleCaseParam USE moduleOutput USE moduleRefParam IMPLICIT NONE INTEGER:: fileID - INTEGER, INTENT(in):: t CHARACTER(*):: fileNameStep, fileNameCollection - IF (t == tInitial) THEN + IF (timeStep == tInitial) THEN !Create collection file WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection) @@ -237,10 +236,11 @@ MODULE moduleMeshOutputVTU !Write iteration file in collection OPEN (fileID + 1, file = path // folder // '/' // fileNameCollection, ACCESS='APPEND') - WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '' + WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") & + '' !Close collection file - IF (t == tFinal) THEN + IF (timeStep == tFinal) THEN WRITE (fileID + 1, "(2X, A)") '' WRITE (fileID + 1, "(A)") '' @@ -307,21 +307,21 @@ MODULE moduleMeshOutputVTU END SUBROUTINE writeAverage - SUBROUTINE printOutputVTU(self,t) + SUBROUTINE printOutputVTU(self) USE moduleMesh USE moduleSpecies USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: i, fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection fileID = 60 DO i = 1, nSpecies - fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t) + fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -337,28 +337,27 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END DO END SUBROUTINE printOutputVTU - SUBROUTINE printCollVTU(self,t) + SUBROUTINE printCollVTU(self) USE moduleMesh USE moduleOutput USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection - CHARACTER (LEN=iterationDigits):: tstring fileID = 62 IF (collOutput) THEN - fileName = formatFileName(prefix, 'Collisions', 'vtu', t) + fileName = formatFileName(prefix, 'Collisions', 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -374,26 +373,26 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END IF END SUBROUTINE printCollVTU - SUBROUTINE printEMVTU(self, t) + SUBROUTINE printEMVTU(self) USE moduleMesh USE moduleMeshInoutCommon + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t INTEGER:: fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection fileID = 64 IF (emOutput) THEN - fileName = formatFileName(prefix, 'EMField', 'vtu', t) + fileName = formatFileName(prefix, 'EMField', 'vtu', timeStep) WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (fileID, file = path // folder // '/' // fileName) @@ -409,7 +408,7 @@ MODULE moduleMeshOutputVTU !Write collection file for time plotting fileNameCollection = formatFileName('Collection', 'EMField', 'pvd') - CALL writeCollection(fileID, t, fileName, filenameCollection) + CALL writeCollection(fileID, fileName, filenameCollection) END IF diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index ae53aa0..f933dd1 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -372,10 +372,9 @@ MODULE moduleMesh END SUBROUTINE connectMesh_interface !Prints number of collisions in each cell - SUBROUTINE printColl_interface(self, t) + SUBROUTINE printColl_interface(self) IMPORT meshGeneric CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printColl_interface @@ -403,18 +402,16 @@ MODULE moduleMesh ABSTRACT INTERFACE !Prints Species data - SUBROUTINE printOutput_interface(self, t) + SUBROUTINE printOutput_interface(self) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printOutput_interface !Prints EM info - SUBROUTINE printEM_interface(self, t) + SUBROUTINE printEM_interface(self) IMPORT meshParticles CLASS(meshParticles), INTENT(in):: self - INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface @@ -789,7 +786,7 @@ MODULE moduleMesh END FUNCTION findCellBrute !Computes collisions in element - SUBROUTINE doCollisions(self, t) + SUBROUTINE doCollisions(self) USE moduleCollisions USE moduleSpecies USE moduleList @@ -797,10 +794,10 @@ MODULE moduleMesh USE moduleRandom USE moduleOutput USE moduleMath + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(meshGeneric), INTENT(inout), TARGET:: self - INTEGER, INTENT(in):: t INTEGER:: e CLASS(meshCell), POINTER:: cell INTEGER:: k, i, j @@ -816,7 +813,7 @@ MODULE moduleMesh REAL(8):: rnd_real !Random number for collision INTEGER:: rnd_int !Random number for collision - IF (MOD(t, everyColl) == 0) THEN + IF (MOD(timeStep, everyColl) == 0) THEN !Collisions need to be performed in this iteration !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j) DO e=1, self%numCells diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 18c7fbb..5b96be2 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -54,7 +54,7 @@ MODULE moduleInject INTEGER:: id CHARACTER(:), ALLOCATABLE:: name REAL(8):: vMod !Velocity (module) - REAL(8):: T(1:3) !Temperature + REAL(8):: temperature(1:3) !Temperature REAL(8):: n(1:3) !Direction of injection LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n INTEGER:: nParticles !Number of particles to introduce each time step @@ -76,7 +76,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge) + SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge) USE moduleMesh USE moduleRefParam USE moduleConstParam @@ -87,7 +87,7 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i - REAL(8), INTENT(in):: v, n(1:3), T(1:3) + REAL(8), INTENT(in):: v, n(1:3), temperature(1:3) INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge REAL(8):: tauInject REAL(8), INTENT(in):: flow @@ -97,10 +97,10 @@ MODULE moduleInject INTEGER:: nVolColl REAL(8):: fluxPerStep = 0.D0 - self%id = i - self%vMod = v / v_ref - self%n = n / NORM2(n) - self%T = T / T_ref + self%id = i + self%vMod = v / v_ref + self%n = n / NORM2(n) + self%temperature = temperature / T_ref !Gets the edge elements from which particles are injected DO e = 1, mesh%numEdges phSurface(e) = mesh%edges(e)%obj%physicalSurface @@ -232,23 +232,23 @@ MODULE moduleInject END SUBROUTINE doInjects - SUBROUTINE initVelDistMaxwellian(velDist, T, m) + SUBROUTINE initVelDistMaxwellian(velDist, temperature, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: T, m + REAL(8), INTENT(in):: temperature, m - velDist = velDistMaxwellian(vTh = DSQRT(T/m)) + velDist = velDistMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistMaxwellian - SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m) + SUBROUTINE initVelDistHalfMaxwellian(velDist, temperature, m) IMPLICIT NONE CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist - REAL(8), INTENT(in):: T, m + REAL(8), INTENT(in):: temperature, m - velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m)) + velDist = velDistHalfMaxwellian(vTh = DSQRT(temperature/m)) END SUBROUTINE initVelDistHalfMaxwellian diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90 index c29eeda..754d56a 100644 --- a/src/modules/moduleProbe.f90 +++ b/src/modules/moduleProbe.f90 @@ -27,7 +27,7 @@ MODULE moduleProbe CONTAINS !Functions for probeDistFunc type - SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, timeStep) + SUBROUTINE init(self, id, speciesName, r, v1, v2, v3, points, everyTimeStep) USE moduleCaseParam USE moduleRefParam USE moduleSpecies @@ -41,7 +41,7 @@ MODULE moduleProbe REAL(8), INTENT(in):: r(1:3) REAL(8), INTENT(in):: v1(1:2), v2(1:2), v3(1:2) INTEGER, INTENT(in):: points(1:3) - REAL(8), INTENT(in):: timeStep + REAL(8), INTENT(in):: everyTimeStep INTEGER:: sp, i REAL(8):: dv(1:3) @@ -91,11 +91,11 @@ MODULE moduleProbe 1:self%nv(3))) !Number of iterations between output - IF (timeStep == 0.D0) THEN + IF (everyTimeStep == 0.D0) THEN self%every = 1 ELSE - self%every = NINT(timeStep/ tauMin / ti_ref) + self%every = NINT(everyTimeStep/ tauMin / ti_ref) END IF @@ -189,13 +189,13 @@ MODULE moduleProbe END SUBROUTINE calculate - SUBROUTINE output(self, t) + SUBROUTINE output(self) USE moduleOutput USE moduleRefParam + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE CLASS(probeDistFunc), INTENT(inout):: self - INTEGER, INTENT(in):: t CHARACTER (LEN=iterationDigits):: tstring CHARACTER (LEN=3):: pstring CHARACTER(:), ALLOCATABLE:: filename @@ -204,14 +204,14 @@ MODULE moduleProbe !Divide by the velocity cube volume self%f = self%f * self%dvInv - WRITE(tstring, iterationFormat) t + WRITE(tstring, iterationFormat) timeStep WRITE(pstring, "(I3.3)") self%id fileName='Probe_' // tstring// '_f_' // pstring // '.dat' WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (10, file = path // folder // '/' // fileName) WRITE(10, "(A1, 1X, A)") "# ", self%species%name WRITE(10, "(A6, 3(ES15.6E3), A2)") "# r = ", self%r(:)*L_ref, " m" - WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(t)*tauMin*ti_ref, " s" + WRITE(10, "(A6, ES15.6E3, A2)") "# t = ", REAL(timeStep)*tauMin*ti_ref, " s" WRITE(10, "(A1, A19, 3(A20))") "#", "v1 (m s^-1)", "v2 (m s^-1)", "v3 (m s^-1)", "f" DO i = 1, self%nv(1) DO j = 1, self%nv(2) @@ -252,15 +252,15 @@ MODULE moduleProbe END SUBROUTINE doProbes - SUBROUTINE outputProbes(t) + SUBROUTINE outputProbes() + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes IF (probe(i)%update) THEN - CALL probe(i)%output(t) + CALL probe(i)%output() END IF @@ -268,15 +268,15 @@ MODULE moduleProbe END SUBROUTINE outputProbes - SUBROUTINE resetProbes(t) + SUBROUTINE resetProbes() + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: i DO i = 1, nProbes probe(i)%f = 0.D0 - probe(i)%update = t == tFinal .OR. t == tInitial .OR. MOD(t, probe(i)%every) == 0 + probe(i)%update = timeStep == tFinal .OR. timeStep == tInitial .OR. MOD(timeStep, probe(i)%every) == 0 END DO diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index 18fbb7f..e6dc91f 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -160,12 +160,12 @@ MODULE moduleOutput END SUBROUTINE calculateOutput - SUBROUTINE printTime(t, first) + SUBROUTINE printTime(first) USE moduleSpecies USE moduleCompTime + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t LOGICAL, INTENT(in), OPTIONAL:: first CHARACTER(:), ALLOCATABLE:: fileName @@ -187,7 +187,7 @@ MODULE moduleOutput OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - WRITE (20, "(I10, I10, 7(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField + WRITE (20, "(I10, I10, 7(ES20.6E3))") timeStep, nPartOld, tStep, tPush, tReset, tColl, tCoul, tWeight, tEMField CLOSE(20) diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index a141eac..30a23ad 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -1,6 +1,7 @@ !Module to solve the electromagnetic field MODULE moduleEM USE moduleMesh + USE moduleTable IMPLICIT NONE ! Array of pointers to nodes. @@ -11,57 +12,155 @@ MODULE moduleEM END TYPE meshNodePointer - - ! TODO: Make this a derived type. - TYPE:: boundaryEM - CHARACTER(:), ALLOCATABLE:: typeEM - INTEGER:: physicalSurface - REAL(8):: potential + ! Generic type for electromagnetic boundary conditions + TYPE, PUBLIC, ABSTRACT:: boundaryEMGeneric + INTEGER:: nNodes TYPE(meshNodePointer), ALLOCATABLE:: nodes(:) CONTAINS - PROCEDURE, PASS:: apply + PROCEDURE(applyEM_interface), DEFERRED, PASS:: apply !PROCEDURE, PASS:: update !only for time dependent boundary conditions or maybe change apply????? That might be better. - END TYPE boundaryEM + END TYPE boundaryEMGeneric + + ABSTRACT INTERFACE + ! Apply boundary condition to the load vector for the Poission equation + SUBROUTINE applyEM_interface(self, vectorF) + IMPORT boundaryEMGeneric + CLASS(boundaryEMGeneric), INTENT(in):: self + REAL(8), INTENT(inout):: vectorF(:) + + END SUBROUTINE applyEM_interface + + END INTERFACE + + TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichlet + REAL(8):: potential + + CONTAINS + ! boundaryEMGeneric DEFERRED PROCEDURES + PROCEDURE, PASS:: apply => applyDirichlet + + END TYPE boundaryEMDirichlet + + TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime + REAL(8):: potential + TYPE(table1D):: temporalProfile + + CONTAINS + ! boundaryEMGeneric DEFERRED PROCEDURES + PROCEDURE, PASS:: apply => applyDirichletTime + + END TYPE boundaryEMDirichletTime + + ! Container for boundary conditions + TYPE:: boundaryEMCont + CLASS(boundaryEMGeneric), ALLOCATABLE:: obj + + END TYPE boundaryEMCont INTEGER:: nBoundaryEM - TYPE(boundaryEM), ALLOCATABLE:: boundEM(:) + TYPE(boundaryEMCont), ALLOCATABLE:: boundaryEM(:) !Information of charge and reference parameters for rho vector REAL(8), ALLOCATABLE:: qSpecies(:) CONTAINS - !Apply boundary conditions to the K matrix for Poisson's equation - SUBROUTINE apply(self, edge) + ! Initialize Dirichlet boundary condition + SUBROUTINE initDirichlet(self, physicalSurface, potential) + USE moduleMesh + USE moduleRefParam, ONLY: Volt_ref + IMPLICIT NONE + + CLASS(boundaryEMGeneric), ALLOCATABLE, INTENT(out):: self + INTEGER:: physicalSurface + REAL(8), INTENT(in):: potential + CLASS(meshEdge), POINTER:: edge + INTEGER, ALLOCATABLE:: nodes(:), nodesEdge(:) + INTEGER:: nNodes, nodesNew + INTEGER:: e, n + + ! Allocate boundary edge + ALLOCATE(boundaryEMDirichlet:: self) + + SELECT TYPE(self) + TYPE IS(boundaryEMDirichlet) + self%potential = potential / Volt_ref + + !TODO: This is going into a function + !Temporal array to hold nodes + ALLOCATE(nodes(0)) + + ! Loop thorugh the edges and identify those that are part of the boundary + DO e = 1, mesh%numEdges + edge => mesh%edges(e)%obj + IF (edge%physicalSurface == physicalSurface) THEN + ! Edge is of the right boundary index + ! Get nodes in the edge + nNodes = edge%nNodes + nodesEdge = edge%getNodes(nNodes) + ! Collect all nodes that are not already in the temporal array + DO n = 1, nNodes + IF (ANY(nodes == nodesEdge(n))) THEN + ! Node already in array, skip + CYCLE + + ELSE + ! If not, add element to array of nodes + nodes = [nodes, nodesEdge(n)] + + END IF + + END DO + + END IF + + END DO + + ! Point boundary to nodes + nNodes = SIZE(nodes) + ALLOCATE(self%nodes(nNodes)) + DO n = 1, nNodes + self%nodes(n)%obj => mesh%nodes(nodes(n))%obj + + END DO + + END SELECT + + END SUBROUTINE initDirichlet + + !Apply Dirichlet boundary condition to the poisson equation + SUBROUTINE applyDirichlet(self, vectorF) USE moduleMesh IMPLICIT NONE - CLASS(boundaryEM), INTENT(in):: self - CLASS(meshEdge):: edge - INTEGER:: nNodes - INTEGER, ALLOCATABLE:: nodes(:) - INTEGER:: n + CLASS(boundaryEMDirichlet), INTENT(in):: self + REAL(8), INTENT(inout):: vectorF(:) + INTEGER:: n, ni - nNodes = 1 - nNodes = edge%nNodes - nodes = edge%getNodes(nNodes) - - DO n = 1, nNodes - SELECT CASE(self%typeEM) - CASE ("dirichlet") - mesh%K(nodes(n), :) = 0.D0 - mesh%K(nodes(n), nodes(n)) = 1.D0 - - ! TODO: Change this to pointer - mesh%nodes(nodes(n))%obj%emData%type = self%typeEM - mesh%nodes(nodes(n))%obj%emData%phi = self%potential - - END SELECT + DO n = 1, self%nNodes + self%nodes(n)%obj%emData%phi = self%potential END DO - END SUBROUTINE apply + END SUBROUTINE applyDirichlet + + !Apply Dirichlet boundary condition with time temporal profile + SUBROUTINE applyDirichletTime(self, vectorF) + USE moduleMesh + USE moduleCaseParam, ONLY: timeStep, tauMin + IMPLICIT NONE + + CLASS(boundaryEMDirichletTime), INTENT(in):: self + REAL(8), INTENT(inout):: vectorF(:) + INTEGER:: n, ni + + DO n = 1, self%nNodes + self%nodes(n)%obj%emData%phi = self%potential + + END DO + + END SUBROUTINE applyDirichletTime !Assemble the source vector based on the charge density to solve Poisson's equation SUBROUTINE assembleSourceVector(vectorF) @@ -74,13 +173,14 @@ MODULE moduleEM INTEGER, ALLOCATABLE:: nodes(:) REAL(8), ALLOCATABLE:: rho(:) INTEGER:: nNodes - INTEGER:: e, i, ni + INTEGER:: e, i, ni, b CLASS(meshNode), POINTER:: node !$OMP SINGLE vectorF = 0.D0 !$OMP END SINGLE + ! Calculate charge density in each node !$OMP DO REDUCTION(+:vectorF) DO e = 1, mesh%numCells nNodes = mesh%cells(e)%obj%nNodes @@ -112,18 +212,12 @@ MODULE moduleEM !$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 + !$OMP SINGLE + DO b = 1, nBoundaryEM + CALL boundaryEM(b)%obj%apply(vectorF) END DO - !$OMP END DO + !$OMP END SINGLE END SUBROUTINE assembleSourceVector diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 0b2837f..5928508 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -491,47 +491,46 @@ MODULE moduleSolver END SUBROUTINE updateParticleCell !Update the information about if a species needs to be moved this iteration - SUBROUTINE updatePushSpecies(self, t) + SUBROUTINE updatePushSpecies(self) USE moduleSpecies + USE moduleCaseparam, ONLY: timeStep 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 + self%pusher(s)%pushSpecies = MOD(timeStep, self%pusher(s)%every) == 0 END DO END SUBROUTINE updatePushSpecies !Output the different data and information - SUBROUTINE doOutput(t) + SUBROUTINE doOutput() USE moduleMesh USE moduleOutput USE moduleSpecies USE moduleCompTime USE moduleProbe + USE moduleCaseParam, ONLY: timeStep IMPLICIT NONE - INTEGER, INTENT(in):: t - - CALL outputProbes(t) + CALL outputProbes() counterOutput = counterOutput + 1 IF (counterOutput >= triggerOutput .OR. & - t == tFinal .OR. t == tInitial) THEN + timeStep == tFinal .OR. timeStep == tInitial) THEN !Resets output counter counterOutput=0 - CALL mesh%printOutput(t) - IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t) - CALL mesh%printEM(t) - WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal + CALL mesh%printOutput() + IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl() + CALL mesh%printEM() + WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", timeStep, "/", tFinal WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld - IF (t == 0) THEN + IF (timeStep == 0) THEN WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms" ELSE @@ -549,34 +548,32 @@ MODULE moduleSolver counterCPUTime = counterCPUTime + 1 IF (counterCPUTime >= triggerCPUTime .OR. & - t == tFinal .OR. t == tInitial) THEN + timeStep == tFinal .OR. timeStep == tInitial) THEN !Reset CPU Time counter counterCPUTime = 0 - CALL printTime(t, t == 0) + CALL printTime(timeStep == 0) END IF !Output average values - IF (useAverage .AND. t == tFinal) THEN + IF (useAverage .AND. timeStep == tFinal) THEN CALL mesh%printAverage() END IF END SUBROUTINE doOutput - SUBROUTINE doAverage(t) + SUBROUTINE doAverage() USE moduleAverage USE moduleMesh IMPLICIT NONE - INTEGER, INTENT(in):: t INTEGER:: tAverage, n - IF (useAverage) THEN - tAverage = t - tAverageStart + tAverage = timeStep - tAverageStart IF (tAverage == 1) THEN !First iteration in which average scheme is used