diff --git a/data/collisions/IO_e-Kr.dat b/data/collisions/IO_e-Kr.dat new file mode 100644 index 0000000..532ca4d --- /dev/null +++ b/data/collisions/IO_e-Kr.dat @@ -0,0 +1,52 @@ +# D108525 "refs": {"B56": {"note": "CLM-R294 (1989)"}} +# Relative energy (eV) cross section (m^2) +1.40E+01 0 +1.62E+01 7.249E-21 +1.88E+01 1.199E-20 +2.18E+01 1.644E-20 +2.53E+01 2.1E-20 +2.94E+01 2.542E-20 +3.41E+01 2.937E-20 +3.95E+01 3.26E-20 +4.58E+01 3.499E-20 +5.32E+01 3.653E-20 +6.17E+01 3.726E-20 +7.15E+01 3.728E-20 +8.29E+01 3.671E-20 +9.62E+01 3.566E-20 +1.12E+02 3.426E-20 +1.29E+02 3.259E-20 +1.50E+02 3.075E-20 +1.74E+02 2.881E-20 +2.02E+02 2.682E-20 +2.34E+02 2.484E-20 +2.72E+02 2.289E-20 +3.15E+02 2.101E-20 +3.65E+02 1.922E-20 +4.24E+02 1.751E-20 +4.91E+02 1.592E-20 +5.70E+02 1.443E-20 +6.61E+02 1.305E-20 +7.67E+02 1.177E-20 +8.89E+02 1.06E-20 +1.03E+03 9.526E-21 +1.20E+03 8.547E-21 +1.39E+03 7.658E-21 +1.61E+03 6.851E-21 +1.87E+03 6.121E-21 +2.16E+03 5.462E-21 +2.51E+03 4.868E-21 +2.91E+03 4.334E-21 +3.38E+03 3.855E-21 +3.92E+03 3.426E-21 +4.54E+03 3.041E-21 +5.27E+03 2.698E-21 +6.11E+03 2.391E-21 +7.09E+03 2.118E-21 +8.22E+03 1.875E-21 +9.53E+03 1.658E-21 +1.11E+04 1.466E-21 +1.28E+04 1.295E-21 +1.49E+04 1.143E-21 +1.72E+04 1.009E-21 +2.00E+04 8.898E-22 diff --git a/data/collisions/IO_e-Xe.dat b/data/collisions/IO_e-Xe.dat new file mode 100644 index 0000000..3067578 --- /dev/null +++ b/data/collisions/IO_e-Xe.dat @@ -0,0 +1,52 @@ +# EL cross sections extracted from PROGRAM MAGBOLTZ, VERSION 7.1 JUNE 2004 www.lxcat.net/Biagi-v7.1 +# Relative energy (eV) cross section (m^2) +1.21E+01 0 +1.41E+01 3.923E-21 +1.64E+01 1.194E-20 +1.91E+01 2.1E-20 +2.22E+01 2.946E-20 +2.58E+01 3.65E-20 +3.00E+01 4.185E-20 +3.49E+01 4.552E-20 +4.06E+01 4.766E-20 +4.72E+01 4.85E-20 +5.49E+01 4.828E-20 +6.39E+01 5.031E-20 +7.43E+01 5.1E-20 +8.64E+01 5.1E-20 +1.01E+02 5.032E-20 +1.17E+02 4.906E-20 +1.36E+02 4.732E-20 +1.58E+02 4.521E-20 +1.84E+02 4.283E-20 +2.14E+02 4.029E-20 +2.49E+02 3.764E-20 +2.90E+02 3.497E-20 +3.37E+02 3.233E-20 +3.92E+02 2.975E-20 +4.56E+02 2.726E-20 +5.31E+02 2.489E-20 +6.17E+02 2.266E-20 +7.18E+02 2.056E-20 +8.35E+02 1.861E-20 +9.72E+02 1.68E-20 +1.13E+03 1.514E-20 +1.32E+03 1.361E-20 +1.53E+03 1.221E-20 +1.78E+03 1.094E-20 +2.07E+03 9.781E-21 +2.41E+03 8.735E-21 +2.80E+03 7.789E-21 +3.26E+03 6.938E-21 +3.79E+03 6.171E-21 +4.41E+03 5.484E-21 +5.13E+03 4.868E-21 +5.97E+03 4.316E-21 +6.94E+03 3.824E-21 +8.07E+03 3.385E-21 +9.39E+03 2.994E-21 +1.09E+04 2.646E-21 +1.27E+04 2.336E-21 +1.48E+04 2.062E-21 +1.72E+04 1.818E-21 +2.00E+04 1.602E-21 diff --git a/src/modules/common/moduleRandom.f90 b/src/modules/common/moduleRandom.f90 index ae5c548..156f5e4 100644 --- a/src/modules/common/moduleRandom.f90 +++ b/src/modules/common/moduleRandom.f90 @@ -47,24 +47,42 @@ MODULE moduleRandom END FUNCTION randomIntAB - !Returns a random number in a Maxwellian distribution of mean 0 and width 1 + !Returns a random number in a Maxwellian distribution of mean 0 and width 1 with the Box-Muller Method FUNCTION randomMaxwellian() RESULT(rnd) USE moduleConstParam, ONLY: PI IMPLICIT NONE REAL(8):: rnd - REAL(8):: x, y + REAL(8):: v1, v2, Rsquare + + Rsquare = 1.D0 + do while (Rsquare >= 1.D0 .and. Rsquare > 0.D0) + v1 = 2.D0 * random() - 1.D0 + v2 = 2.D0 * random() - 1.D0 + Rsquare = v1**2 + v2**2 + + end do + + rnd = v2 * sqrt(-2.D0 * log(Rsquare) / Rsquare) + + END FUNCTION randomMaxwellian + + !Returns a random number in a Maxwellian distribution of mean 0 and width 1 + FUNCTION randomHalfMaxwellian() RESULT(rnd) + IMPLICIT NONE + + REAL(8):: rnd + REAL(8):: x rnd = 0.D0 x = 0.D0 DO WHILE (x == 0.D0) CALL RANDOM_NUMBER(x) END DO - CALL RANDOM_NUMBER(y) - rnd = DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + rnd = DSQRT(-DLOG(x)) - END FUNCTION randomMaxwellian + END FUNCTION randomHalfMaxwellian !Returns a random number weighted with the cumWeight array FUNCTION randomWeighted(cumWeight, sumWeight) RESULT(rnd) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index b5f6881..18c6e65 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -259,6 +259,10 @@ MODULE moduleInput !Read BC CALL readEMBoundary(config) + CASE("ElectrostaticBoltzmann") + !Read BC + CALL readEMBoundary(config) + CASE("ConstantB") !Read BC CALL readEMBoundary(config) @@ -1327,6 +1331,7 @@ MODULE moduleInput USE moduleCaseParam, ONLY: tauMin USE moduleMesh, ONLY: mesh USE moduleSpecies, ONLY: nSpecies + USE moduleRefParam, ONLY: ti_ref IMPLICIT NONE TYPE(json_file), INTENT(inout):: config @@ -1340,8 +1345,11 @@ MODULE moduleInput CALL config%get('average.startTime', tStart, found) IF (found) THEN - tAverageStart = INT(tStart / tauMin) + tAverageStart = INT(tStart / ti_ref / tauMin) + ELSE + tAverageStart = 0 + END IF ALLOCATE(averageScheme(1:mesh%numNodes)) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index dbc8b25..126d40c 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -495,34 +495,36 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 0.D0 - f(3) = 0.D0 DO WHILE(conv > 1.D-4) - dPsi = self%dPsi(XiO, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - invJ = self%invJac(pDer) - fPsi = self%fPsi(XiO, 4) - f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), & - DOT_PRODUCT(fPsi,self%y) /) - r(1:2) - Xi = XiO - MATMUL(invJ, f)/detJ - conv = MAXVAL(DABS(Xi-XiO),1) - XiO = Xi + fPsi = self%fPsi(Xi0, 4) + x0(1) = dot_product(fPsi, self%x) + x0(2) = dot_product(fPsi, self%y) + deltaR = r(1:2) - x0 + dPsi = self%dPsi(Xi0, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ + conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) + Xi0(1:2) = Xi(1:2) END DO @@ -680,8 +682,8 @@ MODULE moduleMesh2DCart dPsi = 0.D0 - dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -826,19 +828,19 @@ MODULE moduleMesh2DCart CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi = 0.D0 - deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - invJ = self%invJac(pDer) - detJ = self%detJac(pDer) - Xi = MATMUL(invJ,deltaR)/detJ + Xi(3) = 0.D0 + deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ END FUNCTION phy2logTria @@ -913,8 +915,8 @@ MODULE moduleMesh2DCart invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) + invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCart diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index ae1eb92..1f5f33c 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -510,34 +510,36 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 0.D0 - f(3) = 0.D0 DO WHILE(conv > 1.D-4) - dPsi = self%dPsi(XiO, 4) - pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer) - invJ = self%invJac(pDer) - fPsi = self%fPsi(XiO, 4) - f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), & - DOT_PRODUCT(fPsi,self%r) /) - r(1:2) - Xi = XiO - MATMUL(invJ, f)/detJ - conv = MAXVAL(DABS(Xi-XiO),1) - XiO = Xi + fPsi = self%fPsi(Xi0, 4) + x0(1) = dot_product(fPsi, self%z) + x0(2) = dot_product(fPsi, self%r) + deltaR = r(1:2) - x0 + dPsi = self%dPsi(Xi0, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ + conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) + Xi0(1:2) = Xi(1:2) END DO @@ -705,8 +707,8 @@ MODULE moduleMesh2DCyl dPsi = 0.D0 - dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -858,19 +860,19 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi = 0.D0 - deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - invJ = self%invJac(pDer) - detJ = self%detJac(pDer) - Xi = MATMUL(invJ,deltaR)/detJ + Xi(3) = 0.D0 + deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ END FUNCTION phy2logTria @@ -948,8 +950,8 @@ MODULE moduleMesh2DCyl invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) + invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCyl diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index b4e0ed2..528cb91 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -199,17 +199,21 @@ MODULE moduleInject END DO + self%nParticles = SUM(self%particlesPerEdge) + 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) - + self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight)) END DO - END IF + self%nParticles = SUM(self%particlesPerEdge) - self%nParticles = SUM(self%particlesPerEdge) + !Rescale weight to match flux + self%weightPerEdge = fluxPerStep * self%surface / (real(self%nParticles)) + + END IF !Scale particles for different species steps IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') @@ -285,7 +289,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - v = self%vTh*randomMaxwellian() + v = sqrt(2.0)*self%vTh*randomMaxwellian() END FUNCTION randomVelMaxwellian @@ -298,10 +302,7 @@ MODULE moduleInject REAL(8):: v v = 0.D0 - DO WHILE (v <= 0.D0) - v = self%vTh*randomMaxwellian() - - END DO + v = sqrt(2.0)*self%vTh*randomHalfMaxwellian() END FUNCTION randomVelHalfMaxwellian @@ -376,7 +377,13 @@ MODULE moduleInject !Assign particle type partInj(n)%species => self%species - direction = self%n + if (all(self%n == 0.D0)) then + direction = randomEdge%normal + + else + direction = self%n + + end if partInj(n)%v = 0.D0 @@ -384,11 +391,12 @@ MODULE moduleInject 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 + !If injecting a no-drift distribution and velocity is negative, reflect + if ((self%vMod == 0.D0) .and. & + (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) diff --git a/src/modules/output/moduleOutput.f90 b/src/modules/output/moduleOutput.f90 index e6dc91f..2af3c70 100644 --- a/src/modules/output/moduleOutput.f90 +++ b/src/modules/output/moduleOutput.f90 @@ -22,7 +22,6 @@ MODULE moduleOutput !Type for EM data in node TYPE emNode - CHARACTER(:), ALLOCATABLE:: type REAL(8):: phi REAL(8):: B(1:3) diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index 126f9a6..7f6891c 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -11,7 +11,6 @@ MODULE moduleEM CONTAINS 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 boundaryEMGeneric @@ -168,8 +167,8 @@ MODULE moduleEM INTEGER:: n, ni DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential - vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi + self%nodes(n)%obj%emData%phi = self%potential + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi END DO @@ -189,15 +188,15 @@ MODULE moduleEM timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin) DO n = 1, self%nNodes - self%nodes(n)%obj%emData%phi = self%potential * timeFactor - vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi + self%nodes(n)%obj%emData%phi = self%potential * timeFactor + vectorF(self%nodes(n)%obj%n) = self%nodes(n)%obj%emData%phi END DO END SUBROUTINE applyDirichletTime !Assemble the source vector based on the charge density to solve Poisson's equation - SUBROUTINE assembleSourceVector(vectorF) + SUBROUTINE assembleSourceVector(vectorF, n_e) USE moduleMesh USE moduleRefParam IMPLICIT NONE @@ -206,6 +205,7 @@ MODULE moduleEM REAL(8), ALLOCATABLE:: localF(:) INTEGER, ALLOCATABLE:: nodes(:) REAL(8), ALLOCATABLE:: rho(:) + REAL(8), INTENT(in), OPTIONAL:: n_e(1:mesh%numNodes) INTEGER:: nNodes INTEGER:: e, i, ni, b CLASS(meshNode), POINTER:: node @@ -214,7 +214,6 @@ MODULE moduleEM 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 @@ -226,6 +225,10 @@ MODULE moduleEM ni = nodes(i) node => mesh%nodes(ni)%obj rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref)) + IF (PRESENT(n_e)) THEN + rho(i) = rho(i) - n_e(i) + + END IF END DO @@ -247,10 +250,10 @@ MODULE moduleEM !Apply boundary conditions !$OMP SINGLE - DO b = 1, nBoundaryEM - CALL boundaryEM(b)%obj%apply(vectorF) + do b = 1, nBoundaryEM + call boundaryEM(b)%obj%apply(vectorF) - END DO + end do !$OMP END SINGLE END SUBROUTINE assembleSourceVector @@ -299,4 +302,86 @@ MODULE moduleEM END SUBROUTINE solveElecField + FUNCTION BoltzmannElectron(phi, n) RESULT(n_e) + USE moduleRefParam + USE moduleConstParam + IMPLICIT NONE + + INTEGER, INTENT(in):: n + REAL(8), INTENT(in):: phi(1:n) + REAL(8):: n_e(1:n) + REAL(8):: n_e0 = 1.0D16, phi_0 = -500.0D0, T_e = 11604.0 + INTEGER:: i + + n_e = n_e0 / n_ref * exp(qe * (phi*Volt_ref - phi_0) / (kb * T_e)) + + RETURN + + END FUNCTION BoltzmannElectron + + SUBROUTINE solveElecFieldBoltzmann() + USE moduleMesh + USE moduleErrors + IMPLICIT NONE + + INTEGER, SAVE:: INFO + INTEGER:: n + REAL(8), ALLOCATABLE, SAVE:: tempF(:) + REAL(8), ALLOCATABLE, SAVE:: n_e(:), phi_old(:), phi(:) + INTEGER:: k + EXTERNAL:: dgetrs + + !$OMP SINGLE + ALLOCATE(tempF(1:mesh%numNodes)) + ALLOCATE(n_e(1:mesh%numNodes)) + ALLOCATE(phi_old(1:mesh%numNodes)) + ALLOCATE(phi(1:mesh%numNodes)) + !$OMP END SINGLE + + !$OMP DO + DO n = 1, mesh%numNodes + phi_old(n) = mesh%nodes(n)%obj%emData%phi + + END DO + !$OMP END DO + + !$OMP SINGLE + DO k = 1, 100 + n_e = BoltzmannElectron(phi_old, mesh%numNodes) + CALL assembleSourceVector(tempF, n_e) + + CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, & + mesh%IPIV, tempF, mesh%numNodes, info) + phi = tempF + + PRINT *, MAXVAL(n_e), MINVAL(n_e) + PRINT *, MAXVAL(phi), MINVAL(phi) + PRINT*, k, "diff = ", MAXVAL(ABS(phi - phi_old)) + phi_old = phi + + END DO + !$OMP END SINGLE + + IF (info == 0) THEN + !Suscessful resolution of Poission equation + !$OMP DO + DO n = 1, mesh%numNodes + mesh%nodes(n)%obj%emData%phi = phi_old(n) + + END DO + !$OMP END DO + + ELSE + !$OMP SINGLE + CALL criticalError('Poisson equation failed', 'solveElecFieldBoltzmann') + !$OMP END SINGLE + + END IF + + !$OMP SINGLE + DEALLOCATE(tempF, n_e, phi_old, phi) + !$OMP END SINGLE + + END SUBROUTINE solveElecFieldBoltzmann + END MODULE moduleEM diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 5928508..df8d5e8 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -138,6 +138,9 @@ MODULE moduleSolver CASE('Electrostatic','ConstantB') self%solveEM => solveElecField + CASE('ElectrostaticBoltzmann') + self%solveEM => solveElecFieldBoltzmann + END SELECT END SUBROUTINE initEM