diff --git a/doc/user-manual/bibliography.bib b/doc/user-manual/bibliography.bib index 878b7ee..9dc17fc 100644 --- a/doc/user-manual/bibliography.bib +++ b/doc/user-manual/bibliography.bib @@ -62,4 +62,27 @@ publisher = {Taylor \& Francis}, } +@Article{sherlock2008monte, + author = {Sherlock, Mark}, + journal = {Journal of Computational Physics}, + title = {A Monte-Carlo method for Coulomb collisions in hybrid plasma models}, + year = {2008}, + number = {4}, + pages = {2286--2292}, + volume = {227}, + groups = {Particle-in-cell}, + publisher = {Elsevier}, +} + +@article{lemons2009small, + title={Small-angle Coulomb collision model for particle-in-cell simulations}, + author={Lemons, Don S and Winske, Dan and Daughton, William and Albright, Brian}, + journal={Journal of Computational Physics}, + volume={228}, + number={5}, + pages={1391--1403}, + year={2009}, + publisher={Elsevier} +} + @Comment{jabref-meta: databaseType:bibtex;} diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index da9471c..c099b01 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 a55c7dd..fccf4a0 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -223,8 +223,15 @@ \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% \subsection{\acrlong{cs}} -% Although not yet implement, a first approach will be soon implemented using Ref.~\cite{higginson2020corrected} as a guideline. + \subsection{\acrlong{cs}} + A simple linearization of the Coulomb operator based on Ref.~\cite{sherlock2008monte} is implemented. + This method assumes that the species involved in the scattering process have a Maxwellian distribution. + The method is made to conserve momentum and kinetic energy based on the approach in Ref.~\cite{lemons2009small} for self (same species) and intra (different species) collisions. + + The user must specify the charged species that will interact together. + The Coulomb logarithm involved in these processes is currently set to a fix value of $10$. + + This method is not valid for situations in which the distribution functions are far from Maxwellian. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Reset of particle array} @@ -619,6 +626,8 @@ make Type of distribution function used to obtain injected particle velocity: \begin{itemize} \item \textbf{Maxwellian}: Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction. + \item \textbf{Half-Maxwellian}: Half-Maxwellian distribution of temperature \textbf{T} and mean \textbf{v} times the value of \textbf{n} in the specified direction. + Only takes into account the positive part of the half-Maxwellian. \item \textbf{Delta}: Dirac's delta distribution function. All particles are injected with velocity \textbf{v} times the value of \textbf{n} in the specified direction. \end{itemize} \item \textbf{T}: Real. @@ -757,6 +766,16 @@ make Only valid for \textbf{ionization} and \textbf{recombination} processes. \end{itemize} \end{itemize} + \item \textbf{Coulomb}: Array of objects. + Contains the information about which species must use the Coulomb linear scattering. + This method assumes a Maxwellian distribution for all species involved. + Each object in the array is defined by: + \begin{itemize} + \item \textbf{species\_i}, \textbf{species\_j}: Character. + Define the two species involved in the collision processes. + Order is indiferent. + + \end{itemize} \end{itemize} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{parallel} diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 3a8a033..e85f5f6 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -74,7 +74,7 @@ PROGRAM fpakc tColl = omp_get_wtime() !$OMP END SINGLE - IF (doMCC) THEN + IF (doMCCollisions) THEN CALL meshForMCC%doCollisions(t) END IF @@ -86,7 +86,7 @@ PROGRAM fpakc tCoul = omp_get_wTime() !$OMP END SINGLE - IF (ASSOCIATED(mesh%doCoulomb)) THEN + IF (doCoulombScattering) THEN CALL mesh%doCoulomb() END IF diff --git a/src/makefile b/src/makefile index c9c747d..247c7d2 100644 --- a/src/makefile +++ b/src/makefile @@ -4,7 +4,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \ - $(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o \ + $(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o $(OBJDIR)/moduleCoulomb.o \ $(OBJDIR)/moduleMeshInoutCommon.o \ $(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \ $(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \ diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index e68a732..9891806 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -155,7 +155,7 @@ MODULE moduleInput sigmaVrel_ref = PI*(r_ref+r_ref)**2*v_ref !reference cross section times velocity ELSE - sigmaVrel_ref = 0.D0 !Assume no collisions + sigmaVrel_ref = L_ref**2 * v_ref END IF @@ -412,7 +412,7 @@ MODULE moduleInput CALL partInitial%add(partNew) !Assign particle to list in volume - IF (doMCC) THEN + IF (listInCells) THEN cell => meshforMCC%cells(partNew%cellColl)%obj CALL OMP_SET_LOCK(cell%lock) CALL cell%listPart_in(sp)%add(partNew) @@ -613,6 +613,7 @@ MODULE moduleInput USE moduleSpecies USE moduleList USE moduleCollisions + USE moduleCoulomb USE moduleErrors USE moduleMesh USE moduleCaseParam @@ -640,10 +641,11 @@ MODULE moduleInput !Firstly, check if the object 'interactions' exists CALL config%info('interactions', found) IF (found) THEN - !Checks if MC collisions have been defined + !Check if MC collisions have been defined CALL config%info('interactions.collisions', found) IF (found) THEN - !Reads collision time step + doMCCollisions = .TRUE. + !Read collision time step CALL config%info('interactions.timeStep', found) IF (found) THEN CALL config%get('interactions.timeStep', tauColl, found) @@ -752,8 +754,35 @@ MODULE moduleInput END IF + !Check if Coulomb scattering is defined + CALL config%info('interactions.Coulomb', found) + IF (found) THEN + + CALL config%info('interactions.Coulomb', found, n_children = nPairs) + IF (nPairs > 0) THEN + nCoulombPairs = nPairs + doCoulombScattering = .TRUE. + ALLOCATE(coulombMatrix(1:nPairs)) + DO i = 1, nPairs + WRITE(iString, '(I2)') i + object = 'interactions.Coulomb(' // 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 coulombMatrix(i)%init(pt_i, pt_j) + + END DO + + END IF + + END IF + END IF + listInCells = doMCCollisions .OR. doCoulombScattering + END SUBROUTINE readInteractions !Reads boundary conditions for the mesh @@ -906,8 +935,6 @@ MODULE moduleInput END IF - doMCC = ASSOCIATED(meshForMCC) - !Get the dimension of the geometry CALL config%get(object // '.dimension', mesh%dimen, found) IF (.NOT. found) THEN @@ -1195,7 +1222,6 @@ MODULE moduleInput CHARACTER(:), ALLOCATABLE:: name REAL(8):: v REAL(8), ALLOCATABLE:: T(:), normal(:) - LOGICAL:: fixDirection REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface @@ -1216,10 +1242,7 @@ MODULE moduleInput CALL config%get(object // '.v', v, found) CALL config%get(object // '.T', T, found) CALL config%get(object // '.n', normal, found) - IF (found) THEN - fixDirection = .TRUE. - ELSE - fixDirection = .FALSE. + IF (.NOT. found) THEN ALLOCATE(normal(1:3)) normal = 0.D0 END IF @@ -1227,7 +1250,7 @@ MODULE moduleInput CALL config%get(object // '.units', units, found) CALL config%get(object // '.physicalSurface', physicalSurface, found) - CALL inject(i)%init(i, v, normal, fixDirection, T, flow, units, sp, physicalSurface) + CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface) CALL readVelDistr(config, inject(i), object) @@ -1306,6 +1329,10 @@ MODULE moduleInput T = inj%T(i) CALL initVelDistMaxwellian(inj%v(i)%obj, t, m) + CASE ("Half-Maxwellian") + T = inj%T(i) + CALL initVelDistHalfMaxwellian(inj%v(i)%obj, t, m) + CASE ("Delta") v = inj%vMod*inj%n(i) CALL initVelDistDelta(inj%v(i)%obj) diff --git a/src/modules/makefile b/src/modules/makefile index 743144f..508bb56 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -1,6 +1,6 @@ OBJS = common.o output.o mesh.o solver.o init.o \ moduleBoundary.o moduleCollisions.o moduleInject.o \ - moduleList.o moduleProbe.o \ + moduleList.o moduleProbe.o moduleCoulomb.o \ moduleSpecies.o all: $(OBJS) @@ -11,7 +11,7 @@ common.o: output.o: moduleSpecies.o common.o $(MAKE) -C output all -mesh.o: moduleCollisions.o moduleBoundary.o output.o common.o +mesh.o: moduleCollisions.o moduleCoulomb.o moduleBoundary.o output.o common.o $(MAKE) -C mesh all solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o @@ -26,6 +26,9 @@ moduleBoundary.o: common.o moduleBoundary.f90 moduleCollisions.o: moduleList.o moduleSpecies.o common.o moduleCollisions.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ +moduleCoulomb.o: moduleSpecies.o common.o moduleCoulomb.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + moduleList.o: common.o moduleSpecies.o moduleList.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index a2ffb7a..d4baedd 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -163,7 +163,7 @@ MODULE moduleMesh2DCyl r2 = self%n2%getCoordinates() self%z = (/r1(1), r2(1)/) self%r = (/r1(2), r2(2)/) - self%weight = SUM(self%r)*5.D-1 + self%weight = r2(2)**2 - r1(2)**2 !Normal vector self%normal = (/ -(self%r(2)-self%r(1)), & self%z(2)-self%z(1) , & diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 9a3c4ef..7e30326 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -300,7 +300,7 @@ MODULE moduleMesh END FUNCTION inside_interface - PURE FUNCTION phy2log_interface(self,r) RESULT(Xi) + PURE FUNCTION phy2log_interface(self,r) RESULT(Xi) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) @@ -387,17 +387,17 @@ MODULE moduleMesh !Array of boundary elements TYPE(meshEdgeCont), ALLOCATABLE:: edges(:) !Global stiffness matrix - REAL(8), ALLOCATABLE, DIMENSION(:,:):: K + REAL(8), ALLOCATABLE, DIMENSION(:,:):: K !Permutation matrix for P L U factorization INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV !PROCEDURES SPECIFIC OF FILE TYPE PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() - PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL() PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL() CONTAINS !GENERIC PROCEDURES PROCEDURE, PASS:: constructGlobalK + PROCEDURE, PASS:: doCoulomb END TYPE meshParticles @@ -424,7 +424,7 @@ MODULE moduleMesh CLASS(meshParticles), INTENT(inout):: self END SUBROUTINE doCoulomb_interface - + !Prints average values SUBROUTINE printAverage_interface(self) IMPORT meshParticles @@ -481,7 +481,11 @@ MODULE moduleMesh !Logical to indicate if an specific mesh for MC Collisions is used LOGICAL:: doubleMesh !Logical to indicate if MCC collisions are performed - LOGICAL:: doMCC + LOGICAL:: doMCCollisions = .FALSE. + !Logical to indicate if Coulomb scattering is performed + LOGICAL:: doCoulombScattering = .FALSE. + !Logica to indicate if particles have to be listed in list inside the cells + LOGICAL:: listInCells = .FALSE. !Complete path for the two meshes CHARACTER(:), ALLOCATABLE:: pathMeshColl, pathMeshParticle @@ -511,7 +515,7 @@ MODULE moduleMesh END DO END DO - + DEALLOCATE(n, localK) END DO @@ -616,7 +620,7 @@ MODULE moduleMesh tensorS = outerProduct(part%v, part%v) sp = part%species%n - + DO i = 1, nNodes node => mesh%nodes(cellNodes(i))%obj CALL OMP_SET_LOCK(node%lock) @@ -650,7 +654,7 @@ MODULE moduleMesh part%Xi = Xi part%n_in = .TRUE. !Assign particle to listPart_in - IF (doMCC) THEN + IF (listInCells) THEN CALL OMP_SET_LOCK(self%lock) sp = part%species%n CALL self%listPart_in(sp)%add(part) @@ -673,7 +677,7 @@ MODULE moduleMesh CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part) !If particle is still inside the domain, call findCell - IF (part%n_in) THEN + IF (part%n_in) THEN IF(PRESENT(oldCell)) THEN CALL self%findCell(part, oldCell) @@ -719,13 +723,13 @@ MODULE moduleMesh INTEGER:: sp found = .FALSE. - + cell => meshColl%cells(part%cellColl)%obj DO WHILE(.NOT. found) Xi = cell%phy2log(part%r) IF (cell%inside(Xi)) THEN part%cellColl = cell%n - IF (doMCC) THEN + IF (listInCells) THEN CALL OMP_SET_LOCK(cell%lock) sp = part%species%n CALL cell%listPart_in(sp)%add(part) @@ -923,7 +927,7 @@ MODULE moduleMesh END DO END IF - + !Deallocate arrays for next collision DEALLOCATE(sigmaVrel, probabilityColl) @@ -946,9 +950,308 @@ MODULE moduleMesh END SUBROUTINE doCollisions SUBROUTINE doCoulomb(self) + USE moduleCoulomb + USE moduleRandom + USE moduleOutput + USE moduleList + USE moduleMath + USE moduleRefParam + USE moduleConstParam IMPLICIT NONE - CLASS(meshParticles), INTENT(inout):: self + CLASS(meshParticles), INTENT(in), TARGET:: self + CLASS(meshCell), POINTER:: cell + TYPE(interactionsCoulomb):: pair + INTEGER:: e + INTEGER:: k + INTEGER:: i, j + INTEGER:: n + INTEGER:: p + TYPE(lNode), POINTER:: partTemp + INTEGER(8), ALLOCATABLE:: cellNodes(:) + CLASS(meshNode), POINTER:: node + TYPE(outputFormat):: output + REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node + REAL(8):: density, velocity(1:3), temperature!values at particle position + REAL(8):: C(1:3), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference + REAL(8):: l, lW, l2 + REAL(8):: GlW, HlW + REAL(8):: normC + REAL(8):: cosThe, sinThe + REAL(8):: cosPhi, sinPhi + REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame + REAL(8):: A, AW + REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W + REAL(8):: theta_per !Random angle for perpendicular direction + REAL(8):: eps = 1.D-12 + REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ij, p_ij + REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ij + REAL(8):: massSum_ij + REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ji, p_ji + REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ji + REAL(8):: massSum_ji + REAL(8):: alpha_num, alpha_den, alpha, beta(1:3) + + + !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp) + DO e = 1, self%numCells + cell => self%cells(e)%obj + cellNodes = cell%getNodes(cell%nNodes) + + ALLOCATE(densityNodes(1:cell%nNodes), & + velocityNodes(1:cell%nNodes, 1:3), & + temperatureNodes(1:cell%nNodes)) + + DO k=1, nCoulombPairs + pair = coulombMatrix(k) + i = pair%sp_i%n + j = pair%sp_j%n + + !Do scattering of particles from species_i due to species j + !Compute background properties of species_j + DO n = 1, cell%nNodes + node => self%nodes(cellNodes(n))%obj + CALL calculateOutput(node%output(j), output, node%v, pair%sp_j) + densityNodes(n) = output%density/n_ref + velocityNodes(n,1:3) = output%velocity(1:3)/v_ref + temperatureNodes(n) = output%temperature/T_ref + + END DO + + ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3)) + ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3)) + ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount)) + !Loop over particles of species_i + partTemp => cell%listPart_in(i)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) + velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) + temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) + + !If cell temperature is too low, skip particle to avoid division by zero + IF (temperature>eps) THEN + l2 = pair%l2_j/temperature + l = SQRT(l2) + + ELSE + partTemp => partTemp%next + + CYCLE + + END IF + + A = pair%A_i*density + + C = partTemp%part%v - velocity + normC = NORM2(C) + + !C_3 = z; C_1, C2 = x, y (per) + C_per = NORM2(C(1:2)) + cosPhi = C(1) / C_per + sinPhi = C(2) / C_per + cosThe = C(3) / normC + sinThe = C_per / normC + + !Rotation matrix to go from W to C + rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column + -sinPhi, cosPhi, 0.D0, & !Second column + sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column + (/ 3, 3 /)) + + !W at start is = (0, 0, normC), so normW = normC + lW = l * normC + GlW = G(lW) + HlW = H(lW) + AW = A / normC + + !Calculate changes in W due to collision process + deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin + deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() + deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() + + !Random angle to distribute perpendicular change in velocity + theta_per = PI2*random() + + !Change W + W(1) = deltaW_per_square * COS(theta_per) + W(2) = deltaW_per_square * SIN(theta_per) + W(3) = normC + deltaW_par + deltaW_par_square + + !Compute changes in velocity for each particle + deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v + mass_ij(p) = pair%sp_i%m*partTemp%part%weight + p_ij(p,1:3) = mass_ij(p)*partTemp%part%v + + !Move to the next particle in the list + partTemp => partTemp%next + p = p + 1 + + END DO + + !Do corresponding collisions + IF (i /= j) THEN + !Do scattering of particles from species_j due to species i + !Compute background properties of species_i + DO n = 1, cell%nNodes + node => self%nodes(cellNodes(n))%obj + CALL calculateOutput(node%output(i), output, node%v, pair%sp_i) + densityNodes(n) = output%density/n_ref + velocityNodes(n,1:3) = output%velocity(1:3)/v_ref + temperatureNodes(n) = output%temperature/T_ref + + END DO + + ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3)) + ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3)) + ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount)) + !Loop over particles of species_j + partTemp => cell%listPart_in(j)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + density = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes) + velocity = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes) + temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes) + + !If cell temperature is too low, skip particle to avoid division by zero + IF (temperature>eps) THEN + l2 = pair%l2_i/temperature + l = SQRT(l2) + + ELSE + partTemp => partTemp%next + + CYCLE + + END IF + A = pair%A_j*density + + C = partTemp%part%v - velocity + normC = NORM2(C) + + !C_3 = z; C_1, C2 = x, y (per) + C_per = NORM2(C(1:2)) + cosPhi = C(1) / C_per + sinPhi = C(2) / C_per + cosThe = C(3) / normC + sinThe = C_per / normC + + !Rotation matrix to go from W to C + rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe, & !First column + -sinPhi, cosPhi, 0.D0, & !Second column + sinThe*cosPhi, sinThe*sinPhi, cosThe /), & !Third column + (/ 3, 3 /)) + + !W at start is = (0, 0, normC), so normW = normC + lW = l * normC + GlW = G(lW) + HlW = H(lW) + AW = A / normC + + !Calculate changes in W due to collision process + deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin + deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian() + deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian() + + !Random angle to distribute perpendicular change in velocity + theta_per = PI2*random() + + !Change W + W(1) = deltaW_per_square * COS(theta_per) + W(2) = deltaW_per_square * SIN(theta_per) + W(3) = normC + deltaW_par + deltaW_par_square + + !Compute changes in velocity for each particle + deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v + mass_ji(p) = pair%sp_j%m*partTemp%part%weight + p_ji(p,1:3) = mass_ji(p)*partTemp%part%v + + !Move to the next particle in the list + partTemp => partTemp%next + p = p + 1 + + END DO + + END IF + + !Calculate correction + !Total mass + massSum_ij = SUM(mass_ij) + massSum_ji = 0.D0 + + !Beta + beta = 0.D0 + DO p = 1, cell%listPart_in(i)%amount + beta = beta + mass_ij(p) * deltaV_ij(p,1:3) + + END DO + + IF (i /= j) THEN + massSum_ji = SUM(mass_ji) + DO p = 1, cell%listPart_in(j)%amount + beta = beta + mass_ji(p) * deltaV_ji(p,1:3) + + END DO + + END IF + + beta = beta / (massSum_ij + massSum_ji) + + !Alpha + alpha_num = 0.D0 + alpha_den = 0.D0 + DO p =1, cell%listPart_in(i)%amount + alpha_num = alpha_num + DOT_PRODUCT(p_ij(p,1:3), deltav_ij(p,1:3) - beta(1:3)) + alpha_den = alpha_den + mass_ij(p) * NORM2(deltav_ij(p,1:3) - beta(1:3))**2 + + END DO + + IF (i /= j) THEN + DO p = 1, cell%listPart_in(j)%amount + alpha_num = alpha_num + DOT_PRODUCT(p_ji(p,1:3), deltav_ji(p,1:3) - beta(1:3)) + alpha_den = alpha_den + mass_ji(p) * NORM2(deltav_ji(p,1:3) - beta(1:3))**2 + + END DO + + END IF + + alpha = -2.D0*alpha_num / alpha_den + + !Apply correction to particles velocity + partTemp => cell%listPart_in(i)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + partTemp%part%v = partTemp%part%v + alpha * (deltaV_ij(p,1:3) - beta(1:3)) + partTemp => partTemp%next + p = p + 1 + + END DO + + IF (i /= j) THEN + partTemp => cell%listPart_in(j)%head + p = 1 + DO WHILE(ASSOCIATED(partTemp)) + partTemp%part%v = partTemp%part%v + alpha * (deltaV_ji(p,1:3) - beta(1:3)) + partTemp => partTemp%next + p = p + 1 + + END DO + + END IF + + DEALLOCATE(deltaV_ij, p_ij, mass_ij) + + IF (i /= j) THEN + DEALLOCATE(deltaV_ji, p_ji, mass_ji) + + END IF + + END DO + + DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes) + + END DO + !$OMP END DO END SUBROUTINE doCoulomb diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 new file mode 100644 index 0000000..166b6b2 --- /dev/null +++ b/src/modules/moduleCoulomb.f90 @@ -0,0 +1,98 @@ +MODULE moduleCoulomb + USE moduleSpecies + IMPLICIT NONE + + INTEGER:: nCoulombPairs = 0 + + !Type for Coulomb iteraction matrix + TYPE:: interactionsCoulomb + CLASS(speciesGeneric), POINTER:: sp_i + CLASS(speciesGeneric), POINTER:: sp_j + REAL(8):: one_plus_massRatio_ij + REAL(8):: lnCoulomb !This can be done a function in the future + REAL(8):: A_i + REAL(8):: A_j + REAL(8):: l2_j + REAL(8):: l2_i + CONTAINS + PROCEDURE, PASS:: init => initInteractionCoulomb + + END TYPE interactionsCoulomb + + !Coulomb collision 'matrix' + TYPE(interactionsCoulomb), ALLOCATABLE:: coulombMatrix(:) + + CONTAINS + PURE REAL(8) FUNCTION G(x) + USE moduleConstParam + IMPLICIT NONE + + REAL(8), INTENT(in):: x + + G = 0.D0 + IF (x /= 0.D0) THEN + G = (ERF(x) - x*2.D0/SQRT(PI)*EXP(-x**2))/(2.D0*x**2) + + END IF + + END FUNCTION G + + PURE REAL(8) FUNCTION H(x) + IMPLICIT NONE + + REAL(8), INTENT(in):: x + + H = ERF(x) - G(x) + + END FUNCTION H + + SUBROUTINE initInteractionCoulomb(self, i, j) + USE moduleSpecies + USE moduleErrors + USE moduleConstParam + USE moduleRefParam + IMPLICIT NONE + + CLASS(interactionsCoulomb), INTENT(out):: self + INTEGER, INTENT(in):: i, j + REAL(8):: Z_i, Z_j + REAL(8):: scaleFactor + + self%sp_i => species(i)%obj + self%sp_j => species(j)%obj + + self%one_plus_massRatio_ij = 1.D0 + self%sp_i%m/self%sp_j%m + + Z_i = 0.D0 + Z_j = 0.D0 + SELECT TYPE(sp => self%sp_i) + TYPE IS (speciesCharged) + Z_i = sp%q + + CLASS DEFAULT + CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb') + + END SELECT + + SELECT TYPE(sp => self%sp_j) + TYPE IS (speciesCharged) + Z_j = sp%q + + CLASS DEFAULT + CALL criticalError('Species ' // sp%name // ' for Coulomb scattering has no charge', 'initInteractionCoulomb') + + END SELECT + + self%lnCoulomb = 10.0 !Make this function dependent + + scaleFactor = (n_ref * qe**4 * ti_ref) / (eps_0**2 * m_ref**2 * v_ref**3) + + self%A_i = Z_i**2*Z_j**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_i%m**2) * scaleFactor !Missing density because it's cell dependent + self%A_j = Z_j**2*Z_i**2*self%lnCoulomb / (2.D0 * PI**2 * self%sp_j%m**2) * scaleFactor !Missing density because it's cell dependent + + self%l2_j = self%sp_j%m / 2.D0 !Missing temperature because it's cell dependent + self%l2_i = self%sp_i%m / 2.D0 !Missing temperature because it's cell dependent + + END SUBROUTINE initInteractionCoulomb + +END MODULE moduleCoulomb diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index e49fb04..4e57083 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -35,6 +35,13 @@ MODULE moduleInject END TYPE velDistMaxwellian + TYPE, EXTENDS(velDistGeneric):: velDistHalfMaxwellian + REAL(8):: vTh !Thermal Velocity + CONTAINS + PROCEDURE, PASS:: randomVel => randomVelHalfMaxwellian + + END TYPE velDistHalfMaxwellian + !Dirac's delta distribution function TYPE, EXTENDS(velDistGeneric):: velDistDelta CONTAINS @@ -68,7 +75,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, fixDirection, T, flow, units, sp, physicalSurface) + SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface) USE moduleMesh USE moduleRefParam USE moduleConstParam @@ -80,7 +87,6 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i REAL(8), INTENT(in):: v, n(1:3), T(1:3) - LOGICAL, INTENT(in):: fixDirection INTEGER, INTENT(in):: sp, physicalSurface REAL(8):: tauInject REAL(8), INTENT(in):: flow @@ -91,12 +97,7 @@ MODULE moduleInject self%id = i self%vMod = v / v_ref - IF (.NOT. fixDirection) THEN - self%n = n / NORM2(n) - ELSE - self%n = 0.D0 - END IF - self%fixDirection = fixDirection + self%n = n / NORM2(n) self%T = T / T_ref self%species => species(sp)%obj tauInject = tau(self%species%n) @@ -212,6 +213,16 @@ MODULE moduleInject END SUBROUTINE initVelDistMaxwellian + SUBROUTINE initVelDistHalfMaxwellian(velDist, T, m) + IMPLICIT NONE + + CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist + REAL(8), INTENT(in):: T, m + + velDist = velDistHalfMaxwellian(vTh = DSQRT(T/m)) + + END SUBROUTINE initVelDistHalfMaxwellian + SUBROUTINE initVelDistDelta(velDist) IMPLICIT NONE @@ -234,6 +245,22 @@ MODULE moduleInject END FUNCTION randomVelMaxwellian + !Random velocity from Half Maxwellian distribution + FUNCTION randomVelHalfMaxwellian(self) RESULT (v) + USE moduleRandom + IMPLICIT NONE + + CLASS(velDistHalfMaxwellian), INTENT(in):: self + REAL(8):: v + v = 0.D0 + + DO WHILE (v <= 0.D0) + v = self%vTh*randomMaxwellian() + + END DO + + END FUNCTION randomVelHalfMaxwellian + !Random velocity from Dirac's delta distribution PURE FUNCTION randomVelDelta(self) RESULT(v) IMPLICIT NONE @@ -305,18 +332,20 @@ MODULE moduleInject !Assign particle type partInj(n)%species => self%species - IF (self%fixDirection) THEN - direction = self%n + direction = self%n - ELSE - direction = randomEdge%normal - - END IF + 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() /) + !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 + !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 diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 69d4055..e539fed 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -70,7 +70,6 @@ MODULE moduleSolver CHARACTER(:), ALLOCATABLE:: pusherType REAL(8):: tau, tauSp - !TODO: Reorganize if Cart pushers are combined SELECT CASE(mesh%dimen) CASE(0) self%pushParticle => push0D @@ -322,7 +321,7 @@ MODULE moduleSolver !$OMP SECTION !Erase the list of particles inside the cell if particles have been pushed - IF (doMCC) THEN + IF (listInCells) THEN DO s = 1, nSpecies DO e = 1, mesh%numCells IF (solver%pusher(s)%pushSpecies) THEN @@ -456,7 +455,7 @@ MODULE moduleSolver CALL partWScheme%unsetLock() !Add particle to cell list sp = part%species%n - IF (doMCC) THEN + IF (listInCells) THEN CALL OMP_SET_lock(cell%lock) CALL cell%listPart_in(sp)%add(newPart) CALL OMP_UNSET_lock(cell%lock)