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)