diff --git a/doc/user-manual/.gitignore b/doc/user-manual/.gitignore index b2ba4d5..a8475d7 100644 --- a/doc/user-manual/.gitignore +++ b/doc/user-manual/.gitignore @@ -6,6 +6,7 @@ *.aux *.ps bibliography.bib.bak +bibliography.bib.sav *.bbl *.blg *.out diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index 2453bb7..ef1494a 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index daf42c3..13be4cd 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -611,7 +611,7 @@ make \begin{itemize} \item \textbf{A}: Ampere. \item \textbf{Am2}: Ampere per square meter. - This value will be multiplied by the square of the reference length to convert to A, it will not consider the surface area of injection. + This value will be multiplied by the surface of injection. \item \textbf{sccm}: Standard cubic centimetre. \item \textbf{part/s}: Particles (real) per second. \end{itemize} @@ -638,6 +638,11 @@ make Temperature in each direction. \item \textbf{physicalSurface}: Integer. Identification of the edge in the mesh file. + \item \textbf{particlesPerEdge}: Integer. + Optional. + Number of particles to be injected by each edge in the numerical domain. + The weight of the particles for each edge will modified by the surface of the edge to ensure the right flux is injected. + If no value is provided, the number of particles to inject per edge will be calculated with the species weight and the surface of the edge respect to the total one. \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{reference} @@ -708,11 +713,15 @@ make \begin{itemize} \item \textbf{species}: Character. Name of species as defined in the object \textbf{species}. - \item \textbf{file}: Character. - Output file from previous run used as an initial state for the species. - The file format must be the same as in \textbf{geometry.meshType} - Initial particles are assumed to have a Maxwellian distribution. - File must be located at \textbf{output.path}. + \item \textbf{file}: Character. + Output file from previous run used as an initial state for the species. + The file format must be the same as in \textbf{geometry.meshType} + Initial particles are assumed to have a Maxwellian distribution. + File must be located at \textbf{output.path}. + \item \textbf{particlesPerCell}: Integer. + Optional. + Initial number of particles per cell. + If not, the number of particles per cell will be assigned based on the species weight and the cell volume. \end{itemize} \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index 0237de1..ae5c548 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -40,10 +40,10 @@ MODULE moduleRandom INTEGER:: rnd REAL(8):: rnd01 - rnd = 0.D0 + rnd = 0 CALL RANDOM_NUMBER(rnd01) - rnd = INT(REAL(b - a) * rnd01) + 1 + rnd = a + FLOOR((b+1-a)*rnd01) END FUNCTION randomIntAB diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 7a1b206..6f1d5bb 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -308,7 +308,7 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: object INTEGER:: nInitial - INTEGER:: i, j, p, e + INTEGER:: i, p, e CHARACTER(LEN=2):: iString CHARACTER(:), ALLOCATABLE:: spName INTEGER:: sp @@ -324,7 +324,8 @@ MODULE moduleInput REAL(8):: densityCen !Mean velocity and temperature at particle position REAL(8):: velocityXi(1:3), temperatureXi - INTEGER:: nNewPart = 0.D0 + INTEGER:: nNewPart = 0 + REAL(8):: weight = 0.D0 CLASS(meshCell), POINTER:: cell TYPE(particle), POINTER:: partNew REAL(8):: vTh @@ -343,6 +344,9 @@ MODULE moduleInput !Reads node values at the nodes filename = path // spFile CALL mesh%readInitial(filename, density, velocity, temperature) + !Check if initial number of particles is given + CALL config%get(object // '.particlesPerCell', nNewPart, found) + !For each volume in the node, create corresponding particles DO e = 1, mesh%numCells !Scale variables @@ -355,7 +359,11 @@ MODULE moduleInput densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar) !Calculate number of particles - nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) + IF (.NOT. found) THEN + nNewPart = FLOOR(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) + + END IF + weight = densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / REAL(nNewPart) !Allocate new particles DO p = 1, nNewPart @@ -392,7 +400,7 @@ MODULE moduleInput partNew%n_in = .TRUE. - partNew%weight = species(sp)%obj%weight + partNew%weight = weight !Assign particle to temporal list of particles CALL partInitial%add(partNew) @@ -915,7 +923,6 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile REAL(8):: volume - CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT object = 'geometry' @@ -1228,6 +1235,7 @@ MODULE moduleInput REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface + INTEGER:: particlesPerEdge INTEGER:: sp CALL config%info('inject', found, n_children = nInject) @@ -1252,8 +1260,10 @@ MODULE moduleInput CALL config%get(object // '.flow', flow, found) CALL config%get(object // '.units', units, found) CALL config%get(object // '.physicalSurface', physicalSurface, found) + particlesPerEdge = 0 + CALL config%get(object // '.particlesPerEdge', particlesPerEdge, found) - CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface) + CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface, particlesPerEdge) CALL readVelDistr(config, inject(i), object) diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 269f157..f400ab0 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -104,6 +104,7 @@ MODULE moduleMesh1DCart USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(out):: self @@ -122,6 +123,8 @@ MODULE moduleMesh1DCart self%x = r1(1) + self%surface = 1.D0 / L_ref**2 + self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index d998267..fd617bd 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -104,6 +104,7 @@ MODULE moduleMesh1DRad USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge1DRad), INTENT(out):: self @@ -122,6 +123,8 @@ MODULE moduleMesh1DRad self%r = r1(1) + self%surface = 1.D0 / L_ref**2 + self%normal = (/ 1.D0, 0.D0, 0.D0 /) !Boundary index diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 7bc94ce..dbc8b25 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -144,6 +144,7 @@ MODULE moduleMesh2DCart USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(out):: self @@ -163,7 +164,7 @@ MODULE moduleMesh2DCart r2 = self%n2%getCoordinates() self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) - self%weight = 1.D0 + self%surface = SQRT((self%x(2) - self%x(1))**2 + (self%y(2) - self%y(1))**2) / L_ref !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -556,6 +557,7 @@ MODULE moduleMesh2DCart !Compute element volume PURE SUBROUTINE volumeQuad(self) + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(inout):: self @@ -573,7 +575,7 @@ MODULE moduleMesh2DCart fPsi = self%fPsi(Xi, 4) !Compute total volume of the cell - self%volume = detJ*4.D0 + self%volume = detJ*4.D0/L_ref !Compute volume per node self%n1%v = self%n1%v + fPsi(1)*self%volume self%n2%v = self%n2%v + fPsi(2)*self%volume diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 7bae6b3..ae1eb92 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -144,6 +144,7 @@ MODULE moduleMesh2DCyl USE moduleSpecies USE moduleBoundary USE moduleErrors + USE moduleConstParam, ONLY: PI IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(out):: self @@ -163,7 +164,15 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - self%weight = DABS(self%r(2)**2 - self%r(1)**2) + !Edge surface + IF (self%z(2) /= self%z(1)) THEN + self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1)) + + ELSE + self%surface = ABS(self%r(2)**2 - self%r(1)**2) + + END IF + self%surface = self%surface * PI !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & @@ -226,11 +235,6 @@ MODULE moduleMesh2DCyl REAL(8):: p1(1:2), p2(1:2) rnd = random() - IF (self%r(1) == 0.D0 .OR. & - self%r(2) == 0.D0) THEN - rnd = rnd**(1.D0/3.D0) - - END IF p1 = (/self%z(1), self%r(1) /) p2 = (/self%z(2), self%r(2) /) @@ -243,7 +247,6 @@ MODULE moduleMesh2DCyl !QUAD FUNCTIONS !Init element SUBROUTINE initCellQuad2DCyl(self, n, p, nodes) - USE moduleRefParam IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(out):: self @@ -589,12 +592,20 @@ MODULE moduleMesh2DCyl fPsi = self%fPsi(Xi, 4) r = DOT_PRODUCT(fPsi,self%r) !Computes total volume of the cell - self%volume = r*detJ*PI8 !4*2*pi - !Computes volume per node - self%n1%v = self%n1%v + fPsi(1)*self%volume - self%n2%v = self%n2%v + fPsi(2)*self%volume - self%n3%v = self%n3%v + fPsi(3)*self%volume - self%n4%v = self%n4%v + fPsi(4)*self%volume + self%volume = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral) + !Computes volume per node. Change the radius point to calculate the area to improve accuracy near the axis. + Xi = (/-5.D-1, -0.25D0, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8 + Xi = (/ 5.D-1, -0.25D0, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8 + Xi = (/ 5.D-1, 0.75D0, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8 + Xi = (/-5.D-1, 0.75D0, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8 END SUBROUTINE volumeQuad diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index c451689..8170df7 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -109,6 +109,7 @@ MODULE moduleMesh3DCart USE moduleBoundary USE moduleErrors USE moduleMath + USE moduleRefParam, ONLY: L_ref IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(out):: self @@ -142,6 +143,8 @@ MODULE moduleMesh3DCart self%normal = crossProduct(vec1, vec2) self%normal = normalize(self%normal) + self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D + !Boundary index self%boundary => boundary(bt) ALLOCATE(self%fBoundary(1:nSpecies)) diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 index 0df1289..1ad6a36 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 @@ -108,6 +108,7 @@ MODULE moduleMeshInputGmsh2 READ(10, *) totalNumElem !Count edges and volume elements + numEdges = 0 SELECT TYPE(self) TYPE IS(meshParticles) self%numEdges = 0 @@ -328,7 +329,7 @@ MODULE moduleMeshInputGmsh2 DO i = 1, numNodes !Reads the density - READ(10, *), e, density(i) + READ(10, *) e, density(i) END DO @@ -339,7 +340,7 @@ MODULE moduleMeshInputGmsh2 DO i = 1, numNodes !Reads the velocity - READ(10, *), e, velocity(i, 1:3) + READ(10, *) e, velocity(i, 1:3) END DO diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index 763517f..184d645 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 @@ -275,6 +275,7 @@ MODULE moduleMeshInputVTU END DO !Count the number of edges + numEdges = 0 SELECT CASE(self%dimen) CASE(3) !Edges are triangles, type 5 in VTK @@ -495,7 +496,7 @@ MODULE moduleMeshInputVTU END SELECT END DO - + !Call mesh connectivity CALL self%connectMesh diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 index 6286cfc..da90b6b 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 +++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 @@ -209,7 +209,7 @@ MODULE moduleMeshOutputVTU WRITE(fileID,"(8X,A)") '' !Electric field WRITE(fileID,"(10X,A, A, A)") '' - WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) + WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) WRITE(fileID,"(10X,A)") '' WRITE(fileID,"(8X,A)") '' @@ -315,9 +315,8 @@ MODULE moduleMeshOutputVTU CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t - INTEGER:: n, i, fileID + INTEGER:: i, fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection - TYPE(outputFormat):: output(1:self%numNodes) fileID = 60 @@ -352,10 +351,9 @@ MODULE moduleMeshOutputVTU CLASS(meshGeneric), INTENT(in):: self INTEGER, INTENT(in):: t - INTEGER:: n, i, fileID + INTEGER:: fileID CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection CHARACTER (LEN=iterationDigits):: tstring - TYPE(outputFormat):: output(1:self%numNodes) fileID = 62 @@ -424,9 +422,8 @@ MODULE moduleMeshOutputVTU IMPLICIT NONE CLASS(meshParticles), INTENT(in):: self - INTEGER:: n, i, fileIDMean, fileIDDeviation + INTEGER:: i, fileIDMean, fileIDDeviation CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation - TYPE(outputFormat):: output(1:self%numNodes) fileIDMean = 66 fileIDDeviation = 67 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index e96ff2a..ae53aa0 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -76,8 +76,8 @@ MODULE moduleMesh CLASS(meshCell), POINTER:: eColl => NULL() !Normal vector REAL(8):: normal(1:3) - !Weight for random injection of particles - REAL(8):: weight = 1.D0 + ! Surface of edge + REAL(8):: surface = 0.D0 !Pointer to boundary type TYPE(boundaryCont), POINTER:: boundary !Array of functions for boundary conditions @@ -613,6 +613,7 @@ MODULE moduleMesh INTEGER:: sp INTEGER:: i CLASS(meshNode), POINTER:: node + REAL(8):: pFraction !Particle fraction cellNodes = self%getNodes(nNodes) fPsi = self%fPsi(part%Xi, nNodes) @@ -623,10 +624,11 @@ MODULE moduleMesh DO i = 1, nNodes node => mesh%nodes(cellNodes(i))%obj + pFraction = fPsi(i)*part%weight CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS + node%output(sp)%den = node%output(sp)%den + pFraction + node%output(sp)%mom(:) = node%output(sp)%mom(:) + pFraction*part%v(:) + node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS CALL OMP_UNSET_LOCK(node%lock) END DO @@ -1023,6 +1025,9 @@ MODULE moduleMesh ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3)) ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3)) ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount)) + deltaV_ij = 0.D0 + p_ij = 0.D0 + mass_ij = 0.D0 !Loop over particles of species_i partTemp => cell%listPart_in(i)%head p = 1 @@ -1107,6 +1112,9 @@ MODULE moduleMesh ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3)) ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3)) ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount)) + deltaV_ji = 0.D0 + p_ji = 0.D0 + mass_ji = 0.D0 !Loop over particles of species_j partTemp => cell%listPart_in(j)%head p = 1 diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index fb01f75..18c7fbb 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -61,8 +61,9 @@ MODULE moduleInject CLASS(speciesGeneric), POINTER:: species !Species of injection INTEGER:: nEdges INTEGER, ALLOCATABLE:: edges(:) !Array with edges - REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability - REAL(8):: sumWeight + INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge + REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge + REAL(8):: surface ! Total surface of injection TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction CONTAINS PROCEDURE, PASS:: init => initInject @@ -75,7 +76,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface) + SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface, particlesPerEdge) USE moduleMesh USE moduleRefParam USE moduleConstParam @@ -87,52 +88,28 @@ MODULE moduleInject 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 + INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge REAL(8):: tauInject REAL(8), INTENT(in):: flow CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et INTEGER:: phSurface(1:mesh%numEdges) 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%species => species(sp)%obj - tauInject = tau(self%species%n) - SELECT CASE(units) - CASE ("sccm") - !Standard cubic centimeter per minute - self%nParticles = INT(flow*sccm2atomPerS*tauInject*ti_ref/species(sp)%obj%weight) - - CASE ("A") - !Input current in Ampers - self%nParticles = INT(flow*tauInject*ti_ref/(qe*species(sp)%obj%weight)) - - CASE ("Am2") - !Input current in Ampers per square meter - self%nParticles = INT(flow*tauInject*ti_ref*L_ref**2/(qe*species(sp)%obj%weight)) - - CASE ("part/s") - !Input current in Ampers - self%nParticles = INT(flow*tauInject*ti_ref/species(sp)%obj%weight) - - CASE DEFAULT - CALL criticalError("No support for units: " // units, 'initInject') - - END SELECT - !Scale particles for different species steps - IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') - + self%id = i + self%vMod = v / v_ref + self%n = n / NORM2(n) + self%T = T / T_ref !Gets the edge elements from which particles are injected 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(self%edges(1:self%nEdges)) + ALLOCATE(self%particlesPerEdge(1:self%nEdges)) + ALLOCATE(self%weightPerEdge(1:self%nEdges)) et = 0 DO e=1, mesh%numEdges IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN @@ -164,15 +141,63 @@ MODULE moduleInject END DO - !Calculates cumulative probability - ALLOCATE(self%cumWeight(1:self%nEdges)) - et = 1 - self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight - DO et = 2, self%nEdges - self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1) + !Calculates total area + self%surface = 0.D0 + DO et = 1, self%nEdges + self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface END DO - self%sumWeight = self%cumWeight(self%nEdges) + + ! Information about species and flux + self%species => species(sp)%obj + tauInject = tau(self%species%n) + ! Convert units + SELECT CASE(units) + CASE ("sccm") + !Standard cubic centimeter per minute + fluxPerStep = flow*sccm2atomPerS + + CASE ("A") + !Current in Ampers + fluxPerStep = flow/qe + + CASE ("Am2") + !Input current in Ampers per square meter + fluxPerStep = flow*self%surface*L_ref**2/qe + + CASE ("part/s") + !Input current in Ampers + fluxPerStep = flow + + CASE DEFAULT + CALL criticalError("No support for units: " // units, 'initInject') + + END SELECT + fluxPerStep = fluxPerStep * tauInject * ti_ref / self%surface + + !Assign particles per edge + IF (particlesPerEdge > 0) THEN + ! Particles per edge defined by the user + self%particlesPerEdge = particlesPerEdge + DO et = 1, self%nEdges + self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge) + + END DO + + ELSE + ! No particles assigned per edge, use the species weight + self%weightPerEdge = self%species%weight + DO et = 1, self%nEdges + self%particlesPerEdge(et) = FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface /self%species%weight) + + END DO + + END IF + + self%nParticles = SUM(self%particlesPerEdge) + + !Scale particles for different species steps + IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') END SUBROUTINE initInject @@ -286,9 +311,8 @@ MODULE moduleInject IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self - INTEGER:: randomX - INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array - INTEGER:: i + INTEGER, SAVE:: nMin + INTEGER:: i, e INTEGER:: n, sp CLASS(meshEdge), POINTER:: randomEdge REAL(8):: direction(1:3) @@ -303,58 +327,62 @@ MODULE moduleInject END IF END DO - nMin = nMin + 1 - nMax = nMin + self%nParticles - 1 - !Particle is considered to be outside the domain - partInj(nMin:nMax)%n_in = .FALSE. !$OMP END SINGLE !$OMP DO - DO n = nMin, nMax - randomX = randomWeighted(self%cumWeight, self%sumWeight) - randomEdge => mesh%edges(self%edges(randomX))%obj - !Random position in edge - partInj(n)%r = randomEdge%randPos() - !Assign weight to particle. - partInj(n)%weight = self%species%weight - !Volume associated to the edge: - IF (ASSOCIATED(randomEdge%e1)) THEN - partInj(n)%cell = randomEdge%e1%n + DO e = 1, self%nEdges + ! Select edge for injection + randomEdge => mesh%edges(self%edges(e))%obj + ! Inject particles in edge + DO i = 1, self%particlesPerEdge(e) + ! Index in the global partInj array + n = nMin - 1 + SUM(self%particlesPerEdge(1:e-1)) + i + !Particle is considered to be outside the domain + partInj(n)%n_in = .FALSE. + !Random position in edge + partInj(n)%r = randomEdge%randPos() + !Assign weight to particle. + partInj(n)%weight = self%weightPerEdge(e) + !Volume associated to the edge: + IF (ASSOCIATED(randomEdge%e1)) THEN + partInj(n)%cell = randomEdge%e1%n - ELSEIF (ASSOCIATED(randomEdge%e2)) THEN - partInj(n)%cell = randomEdge%e2%n + ELSEIF (ASSOCIATED(randomEdge%e2)) THEN + partInj(n)%cell = randomEdge%e2%n - ELSE - CALL criticalError("No Volume associated to edge", 'addParticles') + ELSE + CALL criticalError("No Volume associated to edge", 'addParticles') - END IF - partInj(n)%cellColl = randomEdge%eColl%n - sp = self%species%n + END IF + partInj(n)%cellColl = randomEdge%eColl%n + sp = self%species%n - !Assign particle type - partInj(n)%species => self%species + !Assign particle type + partInj(n)%species => self%species - direction = self%n + direction = self%n - partInj(n)%v = 0.D0 + partInj(n)%v = 0.D0 - partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & - self%v(2)%obj%randomVel(), & - self%v(3)%obj%randomVel() /) + partInj(n)%v = self%vMod*direction + (/ self%v(1)%obj%randomVel(), & + self%v(2)%obj%randomVel(), & + self%v(3)%obj%randomVel() /) - !If velocity is not in the right direction, invert it - IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN - partInj(n)%v = - partInj(n)%v + !If velocity is not in the right direction, invert it + IF (DOT_PRODUCT(direction, partInj(n)%v) < 0.D0) THEN + partInj(n)%v = - partInj(n)%v - END IF + END IF - !Obtain natural coordinates of particle in cell - partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) - !Push new particle with the minimum time step - CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) - !Assign cell to new particle - CALL solver%updateParticleCell(partInj(n)) + !Obtain natural coordinates of particle in cell + partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r) + !Push new particle with the minimum time step + CALL solver%pusher(sp)%pushParticle(partInj(n), tau(sp)) + !Assign cell to new particle + CALL solver%updateParticleCell(partInj(n)) + + END DO END DO !$OMP END DO diff --git a/src/modules/moduleProbe.f90 b/src/modules/moduleProbe.f90 index c7d3cf5..c29eeda 100644 --- a/src/modules/moduleProbe.f90 +++ b/src/modules/moduleProbe.f90 @@ -101,7 +101,7 @@ MODULE moduleProbe !Maximum radius !TODO: Make this an input parameter - self%maxR = 1.D0 + self%maxR = 1.D-2/L_ref !Init the probe lock CALL OMP_INIT_LOCK(self%lock) @@ -148,7 +148,7 @@ MODULE moduleProbe deltaR = NORM2(self%r - part%r) !Only include particle if it is inside the maximum radius - IF (deltaR < self%maxR) THEN + ! IF (deltaR < self%maxR) THEN !find lower index for all dimensions CALL self%findLowerIndex(part%v, i, j, k, inside) @@ -162,28 +162,28 @@ MODULE moduleProbe fk = self%vk(k+1) - part%v(3) fk1 = part%v(3) - self%vk(k) - ! weight = part%weight * DEXP(deltaR/self%maxR) - weight = part%weight + weight = part%weight * DEXP(-deltaR/self%maxR) + ! weight = part%weight !Lock the probe CALL OMP_SET_LOCK(self%lock) !Assign particle weight to distribution function - self%f(i , j , k ) = fi * fj * fk * weight - self%f(i+1, j , k ) = fi1 * fj * fk * weight - self%f(i , j+1, k ) = fi * fj1 * fk * weight - self%f(i+1, j+1, k ) = fi1 * fj1 * fk * weight - self%f(i , j , k+1) = fi * fj * fk1 * weight - self%f(i+1, j , k+1) = fi1 * fj * fk1 * weight - self%f(i , j+1, k+1) = fi * fj1 * fk1 * weight - self%f(i+1, j+1, k+1) = fi1 * fj1 * fk1 * weight + self%f(i , j , k ) = self%f(i , j , k ) + fi * fj * fk * weight + self%f(i+1, j , k ) = self%f(i+1, j , k ) + fi1 * fj * fk * weight + self%f(i , j+1, k ) = self%f(i , j+1, k ) + fi * fj1 * fk * weight + self%f(i+1, j+1, k ) = self%f(i+1, j+1, k ) + fi1 * fj1 * fk * weight + self%f(i , j , k+1) = self%f(i , j , k+1) + fi * fj * fk1 * weight + self%f(i+1, j , k+1) = self%f(i+1, j , k+1) + fi1 * fj * fk1 * weight + self%f(i , j+1, k+1) = self%f(i , j+1, k+1) + fi * fj1 * fk1 * weight + self%f(i+1, j+1, k+1) = self%f(i+1, j+1, k+1) + fi1 * fj1 * fk1 * weight !Unlock the probe CALL OMP_UNSET_LOCK(self%lock) END IF - END IF + ! END IF END IF diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index bdf6b03..d5d0793 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -30,8 +30,9 @@ MODULE moduleEM INTEGER, ALLOCATABLE:: nodes(:) INTEGER:: n + nNodes = 1 nNodes = edge%nNodes - nodes = edge%getNodes(nNodes) + nodes = edge%getNodes(nNodes) DO n = 1, nNodes SELECT CASE(self%typeEM)