From 601103105f469573a453f6334a45db9f7c6cd2f9 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 24 Feb 2023 21:46:01 +0100 Subject: [PATCH] First attempt at Coulomb collisions First attemp for Coulomb collisions based on the moments distribtuions. Still the method is not done and far from being complete but input options and basic math are implemented. --- src/fpakc.f90 | 4 +- src/makefile | 2 +- src/modules/common/moduleRefParam.f90 | 2 +- src/modules/init/moduleInput.f90 | 38 +++++- src/modules/makefile | 7 +- src/modules/mesh/moduleMesh.f90 | 173 ++++++++++++++++++++++++-- src/modules/moduleCoulomb.f90 | 91 ++++++++++++++ src/modules/solver/moduleSolver.f90 | 4 +- 8 files changed, 295 insertions(+), 26 deletions(-) create mode 100644 src/modules/moduleCoulomb.f90 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/common/moduleRefParam.f90 b/src/modules/common/moduleRefParam.f90 index 051190b..e8819ce 100644 --- a/src/modules/common/moduleRefParam.f90 +++ b/src/modules/common/moduleRefParam.f90 @@ -3,7 +3,7 @@ MODULE moduleRefParam !Parameters that define the problem (inputs) REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigmaVrel_ref !Reference parameters for non-dimensional problem - REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref + REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref, e_ref END MODULE moduleRefParam diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index e68a732..7af67db 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -144,6 +144,7 @@ MODULE moduleInput EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field Volt_ref = EF_ref*L_ref !reference voltage B_ref = m_ref / (ti_ref * qe) !reference magnetic field + e_ref = qe/(m_ref*v_ref**2) !reference charge !If a reference cross section is given, it is used CALL config%get(object // '.sigmaVrel', sigmavRel_ref, found) @@ -412,7 +413,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 +614,7 @@ MODULE moduleInput USE moduleSpecies USE moduleList USE moduleCollisions + USE moduleCoulomb USE moduleErrors USE moduleMesh USE moduleCaseParam @@ -640,10 +642,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 +755,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 +936,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 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/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 9a3c4ef..6607fd0 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,152 @@ 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 + INTEGER:: e + INTEGER:: k + INTEGER:: i, j + INTEGER:: n + TYPE(lNode), POINTER:: partTemp + REAL(8):: W(3), dW(2) !Relative velocity between particle and species and its increment + 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), DIMENSION(1:3):: e1, e2, e3 + REAL(8):: delta_par, delta_par_square, delta_per, delta_per_square + REAL(8):: l, lW, AW + REAL(8):: rnd + + + 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 + i = coulombMatrix(k)%sp_i%n + j = coulombMatrix(k)%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, coulombMatrix(k)%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 + + !Loop over particles of species_i + partTemp => cell%listPart_in(i)%head + 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) + l = coulombMatrix(k)%l_j/SQRT(temperature) + + W = partTemp%part%v - velocity + lW = l * NORM2(W) + AW = coulombMatrix(k)%A_ij/NORM2(W) + !Axis of the relative velocity + !First one is parallel to the relative velocity + e1 = normalize(W) + !Second one is perpendicular to it + e2(1) = -e1(2) + e2(2) = e1(1) + e2(3) = 0.D0 + e2 = normalize(e2) + !Third one is perpendicular to the other two + e3 = crossProduct(e2, e1) + e3 = normalize(e3) + + delta_par = -coulombMatrix(k)%A_ij*coulombMatrix(k)%one_plus_massRatio_ij*density*l**2*G(lW) + + delta_par_square = AW*density*G(lW) + + delta_per_square = AW*density*H(lW) + + dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin) + dW(2) = DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) + + rnd = random() + partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + & + SIN(PI2*rnd)*e3) + + partTemp => partTemp%next + + END DO + + IF (i /= j) THEN + !Do scattering of particles from species_j due to species i + DO n = 1, cell%nNodes + node => self%nodes(cellNodes(n))%obj + CALL calculateOutput(node%output(i), output, node%v, coulombMatrix(k)%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 + + partTemp => cell%listPart_in(j)%head + 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) + l = coulombMatrix(k)%l_i/SQRT(temperature) + + W = partTemp%part%v - velocity + lW = l * NORM2(W) + AW = coulombMatrix(k)%A_ji/NORM2(W) + !Axis of the relative velocity + !First one is parallel to the relative velocity + e1 = normalize(W) + !Second one is perpendicular to it + e2(1) = -e1(2) + e2(2) = e1(1) + e2(3) = 0.D0 + e2 = normalize(e2) + !Third one is perpendicular to the other two + e3 = crossProduct(e2, e1) + e3 = normalize(e3) + + delta_par = -coulombMatrix(k)%A_ji*coulombMatrix(k)%one_plus_massRatio_ji*density*l**2*G(lW) + + delta_par_square = AW*density*G(lW) + + delta_per_square = AW*density*H(lW) + + dW(1) = delta_par*tauMin + randomMaxwellian()*SQRT(delta_par_square*tauMin) + dW(2) = DABS(randomMaxwellian()*SQRT(delta_per_square*tauMin)) + + rnd = random() + partTemp%part%v = partTemp%part%v + dW(1)*e1 + dW(2)*(COS(PI2*rnd)*e2 + & + SIN(PI2*rnd)*e3) + + partTemp => partTemp%next + + END DO + + END IF + + END DO + + END DO END SUBROUTINE doCoulomb diff --git a/src/modules/moduleCoulomb.f90 b/src/modules/moduleCoulomb.f90 new file mode 100644 index 0000000..734b898 --- /dev/null +++ b/src/modules/moduleCoulomb.f90 @@ -0,0 +1,91 @@ +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, one_plus_massRatio_ji + REAL(8):: lnCoulomb !This can be done a function in the future + REAL(8):: A_ij, A_ji + REAL(8):: l_j, l_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 + + self%sp_i => species(i)%obj + self%sp_j => species(j)%obj + + self%one_plus_massRatio_ij = 1.D0 + (self%sp_i%weight*self%sp_i%m)/(self%sp_j%weight*self%sp_j%m) + self%one_plus_massRatio_ji = 1.D0 + (self%sp_j%weight*self%sp_j%m)/(self%sp_i%weight*self%sp_i%m) + + 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 = 12.0 + self%A_ij = 8.D0*PI*Z_i**2*Z_j**2*e_ref**4*self%lnCoulomb / self%sp_i%m + self%A_ji = 8.D0*PI*Z_j**2*Z_i**2*e_ref**4*self%lnCoulomb / self%sp_j%m + + self%l_j = SQRT(self%sp_j%m / 2.D0) !Missing temperature because it's cell dependent + self%l_i = SQRT(self%sp_i%m / 2.D0) !Missing temperature because it's cell dependent + + END SUBROUTINE initInteractionCoulomb + +END MODULE moduleCoulomb diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 69d4055..12f9e9a 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -322,7 +322,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 +456,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)