From 9e0d1a7cc75bb79f3634704f699e1ff2b8408672 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sat, 26 Dec 2020 22:45:55 +0100 Subject: [PATCH] Final implementation of ionization process by electron impact. Possibility to input initial species distributions (density, velocity and temperature) via an input file for each species. New moduleRandom includes function to generate random numbers in different ways (still uses) the implicit RANDOM_NUMBER(). --- src/fpakc.f90 | 15 +- src/makefile | 2 +- src/modules/makefile | 9 +- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 26 +++- .../mesh/1DCart/moduleMesh1DCartBoundary.f90 | 13 +- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 21 ++- .../mesh/1DRad/moduleMesh1DRadBoundary.f90 | 13 +- src/modules/mesh/2DCyl/moduleMeshCyl.f90 | 57 ++++++- .../mesh/2DCyl/moduleMeshCylBoundary.f90 | 13 +- src/modules/mesh/moduleMesh.f90 | 24 +-- src/modules/moduleCollisions.f90 | 23 +-- src/modules/moduleInject.f90 | 16 +- src/modules/moduleInput.f90 | 146 +++++++++++++++++- src/modules/moduleList.f90 | 1 + src/modules/moduleRandom.f90 | 69 +++++++++ src/modules/moduleTable.f90 | 1 + 16 files changed, 363 insertions(+), 86 deletions(-) create mode 100644 src/modules/moduleRandom.f90 diff --git a/src/fpakc.f90 b/src/fpakc.f90 index 0320d57..62ab291 100644 --- a/src/fpakc.f90 +++ b/src/fpakc.f90 @@ -2,13 +2,11 @@ PROGRAM fpakc USE moduleInput USE moduleErrors - USE OMP_LIB USE moduleInject USE moduleSolver - USE moduleOutput USE moduleCompTime - USE moduleCollisions - USE moduleMesh + USE moduleCaseParam + USE OMP_LIB IMPLICIT NONE ! t = time step @@ -28,8 +26,17 @@ PROGRAM fpakc !Reads the json configuration file CALL readConfig(inputFile) + !$OMP PARALLEL DEFAULT(SHARED) + !$OMP SINGLE + CALL verboseError("Initial scatter of particles...") + !$OMP END SINGLE + CALL doScatter() + + !$OMP SINGLE CALL verboseError("Calculating initial EM field...") + !$OMP END SINGLE CALL doEMField() + !$OMP END PARALLEL tStep = omp_get_wtime() - tStep diff --git a/src/makefile b/src/makefile index 715fe52..3e151a3 100644 --- a/src/makefile +++ b/src/makefile @@ -3,7 +3,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolv $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ - $(OBJDIR)/moduleEM.o \ + $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \ $(OBJDIR)/moduleMeshCyl.o $(OBJDIR)/moduleMeshCylRead.o $(OBJDIR)/moduleMeshCylBoundary.o \ $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o $(OBJDIR)/moduleMesh1DCartBoundary.o \ $(OBJDIR)/moduleMesh1DRad.o $(OBJDIR)/moduleMesh1DRadRead.o $(OBJDIR)/moduleMesh1DRadBoundary.o diff --git a/src/modules/makefile b/src/modules/makefile index 6c8895b..28fc4ec 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -2,20 +2,20 @@ OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \ moduleOutput.o moduleInput.o moduleSolver.o \ moduleCollisions.o moduleTable.o moduleParallel.o \ - moduleEM.o + moduleEM.o moduleRandom.o all: $(OBJS) mesh.o mesh.o: moduleCollisions.o moduleBoundary.o $(MAKE) -C mesh all -moduleCollisions.o: moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90 +moduleCollisions.o: moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ moduleInput.o: moduleParallel.o moduleRefParam.o moduleCaseParam.o moduleSolver.o moduleInject.o moduleBoundary.o moduleErrors.o moduleSpecies.o moduleInput.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleInject.o: moduleSpecies.o moduleSolver.o moduleInject.f90 +moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90 @@ -24,6 +24,9 @@ moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90 moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ +moduleRandom.o: moduleConstParam.o moduleRandom.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + moduleRefParam.o: moduleConstParam.o moduleRefParam.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 659cbd7..3bd7d05 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -23,7 +23,7 @@ MODULE moduleMesh1DCart CONTAINS PROCEDURE, PASS:: init => initEdge1DCart PROCEDURE, PASS:: getNodes => getNodes1DCart - PROCEDURE, PASS:: randPos => randPos1DCart + PROCEDURE, PASS:: randPos => randPosEdge END TYPE meshEdge1DCart @@ -110,6 +110,7 @@ MODULE moduleMesh1DCart REAL(8):: arNodes(1:2) CONTAINS PROCEDURE, PASS:: init => initVol1DCartSegm + PROCEDURE, PASS:: randPos => randPos1DCartSeg PROCEDURE, PASS:: area => areaSegm PROCEDURE, NOPASS:: fPsi => fPsiSegm PROCEDURE, NOPASS:: dPsi => dPsiSegm @@ -226,13 +227,13 @@ MODULE moduleMesh1DCart END FUNCTION getNodes1DCart !Calculates a 'random' position in edge - FUNCTION randPos1DCart(self) RESULT(r) + FUNCTION randPosEdge(self) RESULT(r) CLASS(meshEdge1DCart), INTENT(in):: self REAL(8):: r(1:3) r = (/ self%x, 0.D0, 0.D0 /) - END FUNCTION randPos1DCart + END FUNCTION randPosEdge !VOLUME FUNCTIONS !SEGMENT FUNCTIONS @@ -265,6 +266,24 @@ MODULE moduleMesh1DCart END SUBROUTINE initVol1DCartSegm + !Calculates a random position in 1D volume + FUNCTION randPos1DCartSeg(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshVol1DCartSegm), INTENT(in):: self + REAL(8):: r(1:3) + REAL(8):: xii(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + xii(1) = random(-1.D0, 1.D0) + xii(2:3) = 0.D0 + + fPsi = self%fPsi(xii) + r(1) = DOT_PRODUCT(fPsi, self%x) + + END FUNCTION randPos1DCartSeg + !Computes element area PURE SUBROUTINE areaSegm(self) IMPLICIT NONE @@ -479,6 +498,7 @@ MODULE moduleMesh1DCart END SUBROUTINE nextElementSegm !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS + !Calculates a random position in 1D volume !Computes the element Jacobian determinant PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ) IMPLICIT NONE diff --git a/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 b/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 index 443c7cc..ce2f229 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90 @@ -66,25 +66,16 @@ SUBMODULE (moduleMesh1DCart) moduleMesh1DCartBoundary SUBROUTINE wallTemperature(edge, part) USE moduleSpecies USE moduleBoundary - USE moduleConstParam, ONLY: PI + USE moduleRandom IMPLICIT NONE CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - REAL(8):: x, y !Modifies particle velocity according to wall temperature SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) TYPE IS(boundaryWallTemperature) - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - - END DO - - CALL RANDOM_NUMBER(y) - - part%v(1) = part%v(1) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + part%v(1) = part%v(1) + bound%vTh*randomMaxwellian() END SELECT diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 9693319..eb78eb3 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -75,7 +75,7 @@ MODULE moduleMesh1DRad PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - END TYPE meshVol1Drad + END TYPE meshVol1DRad ABSTRACT INTERFACE PURE FUNCTION fPsi_interface(xi) RESULT(fPsi) @@ -111,6 +111,7 @@ MODULE moduleMesh1DRad REAL(8):: arNodes(1:2) CONTAINS PROCEDURE, PASS:: init => initVol1DRadSegm + PROCEDURE, PASS:: randPos => randPos1DRadSeg PROCEDURE, PASS:: area => areaRad PROCEDURE, NOPASS:: fPsi => fPsiRad PROCEDURE, NOPASS:: dPsi => dPsiRad @@ -266,6 +267,24 @@ MODULE moduleMesh1DRad END SUBROUTINE initVol1DRadSegm + !Calculates a random position in 1D volume + FUNCTION randPos1DRadSeg(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshVol1DRadSegm), INTENT(in):: self + REAL(8):: r(1:3) + REAL(8):: xii(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + xii(1) = random(-1.D0, 1.D0) + xii(2:3) = 0.D0 + + fPsi = self%fPsi(xii) + r(1) = DOT_PRODUCT(fPsi, self%r) + + END FUNCTION randPos1DRadSeg + !Computes element area PURE SUBROUTINE areaRad(self) IMPLICIT NONE diff --git a/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 b/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 index a9cdca8..75cf27d 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90 @@ -68,25 +68,16 @@ SUBMODULE (moduleMesh1DRad) moduleMesh1DRadBoundary SUBROUTINE wallTemperature(edge, part) USE moduleSpecies USE moduleBoundary - USE moduleConstParam, ONLY: PI + USE moduleRandom IMPLICIT NONE CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part - REAL(8):: x, y !Modifies particle velocity according to wall temperature SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) TYPE IS(boundaryWallTemperature) - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - - END DO - - CALL RANDOM_NUMBER(y) - - part%v(1) = part%v(1) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + part%v(1) = part%v(1) + bound%vTh*randomMaxwellian() END SELECT diff --git a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 index c4f2b55..251573e 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMeshCyl.f90 @@ -31,7 +31,7 @@ MODULE moduleMeshCyl CONTAINS PROCEDURE, PASS:: init => initEdgeCyl PROCEDURE, PASS:: getNodes => getNodesCyl - PROCEDURE, PASS:: randPos => randPosCyl + PROCEDURE, PASS:: randPos => randPosEdge END TYPE meshEdgeCyl @@ -129,6 +129,7 @@ MODULE moduleMeshCyl CONTAINS PROCEDURE, PASS:: init => initVolQuadCyl + PROCEDURE, PASS:: randPos => randPosVolQuad PROCEDURE, PASS:: area => areaQuad PROCEDURE, NOPASS:: fPsi => fPsiQuad PROCEDURE, NOPASS:: dPsi => dPsiQuad @@ -161,6 +162,7 @@ MODULE moduleMeshCyl CONTAINS PROCEDURE, PASS:: init => initVolTriaCyl + PROCEDURE, PASS:: randPos => randPosVolTria PROCEDURE, PASS:: area => areaTria PROCEDURE, NOPASS:: fPsi => fPsiTria PROCEDURE, NOPASS:: dPsi => dPsiTria @@ -274,6 +276,28 @@ MODULE moduleMeshCyl END SUBROUTINE initEdgeCyl + !Random position in quadrilateral volume + FUNCTION randPosVolQuad(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8):: r(1:3) + REAL(8):: xii(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + xii(1) = random(-1.D0, 1.D0) + xii(2) = random(-1.D0, 1.D0) + xii(3) = 0.D0 + + fPsi = self%fPsi(xii) + + r(1) = DOT_PRODUCT(fPsi, self%z) + r(2) = DOT_PRODUCT(fPsi, self%r) + r(3) = 0.D0 + + END FUNCTION randposVolQuad + !Get nodes from edge PURE FUNCTION getNodesCyl(self) RESULT(n) IMPLICIT NONE @@ -287,20 +311,23 @@ MODULE moduleMeshCyl END FUNCTION getNodesCyl !Calculates a random position in edge - FUNCTION randPosCyl(self) RESULT(r) + FUNCTION randPosEdge(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + CLASS(meshEdgeCyl), INTENT(in):: self REAL(8):: rnd REAL(8):: r(1:3) REAL(8):: p1(1:2), p2(1:2) - CALL RANDOM_NUMBER(rnd) + rnd = random() p1 = (/self%z(1), self%r(1) /) p2 = (/self%z(2), self%r(2) /) r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 r(3) = 0.D0 - END FUNCTION randPosCyl + END FUNCTION randPosEdge !VOLUME FUNCTIONS !QUAD FUNCTIONS @@ -709,6 +736,28 @@ MODULE moduleMeshCyl END SUBROUTINE initVolTriaCyl + !Random position in quadrilateral volume + FUNCTION randPosVolTria(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshVolCylTria), INTENT(in):: self + REAL(8):: r(1:3) + REAL(8):: xii(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + xii(1) = random( 0.D0, 1.D0) + xii(2) = random( 0.D0, 1.D0) + xii(3) = 0.D0 + + fPsi = self%fPsi(xii) + + r(1) = DOT_PRODUCT(fPsi, self%z) + r(2) = DOT_PRODUCT(fPsi, self%r) + r(3) = 0.D0 + + END FUNCTION randposVolTria + !Calculates area for triangular element PURE SUBROUTINE areaTria(self) USE moduleConstParam diff --git a/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 b/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 index 556abf7..5410007 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 +++ b/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 @@ -104,28 +104,19 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary SUBROUTINE wallTemperature(edge, part) USE moduleSpecies USE moduleBoundary - USE moduleConstParam, ONLY: PI + USE moduleRandom IMPLICIT NONE CLASS(meshEdge), INTENT(inout):: edge CLASS(particle), INTENT(inout):: part REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2) INTEGER:: i - REAL(8):: x, y !Modifies particle velocity according to wall temperature SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj) TYPE IS(boundaryWallTemperature) DO i = 1, 3 - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - - END DO - - CALL RANDOM_NUMBER(y) - - part%v(i) = part%v(i) + bound%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + part%v(i) = part%v(i) + bound%vTh*randomMaxwellian() END DO diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index a75c939..b0b2bd2 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -70,7 +70,7 @@ MODULE moduleMesh CONTAINS PROCEDURE(initEdge_interface), DEFERRED, PASS:: init PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes - PROCEDURE(randPos_interface), DEFERRED, PASS:: randPos + PROCEDURE(randPosEdge_interface), DEFERRED, PASS:: randPos END TYPE meshEdge @@ -93,12 +93,12 @@ MODULE moduleMesh END FUNCTION - FUNCTION randPos_interface(self) RESULT(r) + FUNCTION randPosEdge_interface(self) RESULT(r) IMPORT:: meshEdge CLASS(meshEdge), INTENT(in):: self REAL(8):: r(1:3) - END FUNCTION randPos_interface + END FUNCTION randPosEdge_interface END INTERFACE @@ -138,9 +138,10 @@ MODULE moduleMesh REAL(8):: totalWeight = 0.D0 CONTAINS PROCEDURE(initVol_interface), DEFERRED, PASS:: init + PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes + PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF - PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF PROCEDURE, PASS:: findCell PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log @@ -215,6 +216,13 @@ MODULE moduleMesh END FUNCTION inside_interface + FUNCTION randPosVol_interface(self) RESULT(r) + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + REAL(8):: r(1:3) + + END FUNCTION randPosVol_interface + END INTERFACE !Containers for volumes in the mesh @@ -376,13 +384,13 @@ MODULE moduleMesh USE moduleSpecies USE moduleList use moduleRefParam + USE moduleRandom IMPLICIT NONE CLASS(meshVol), INTENT(inout):: self INTEGER:: nPart !Number of particles inside the cell REAL(8):: pMax !Maximum probability of collision INTEGER:: rnd !random index - REAL(8):: rndReal TYPE(particle), POINTER:: part_i, part_j INTEGER:: n !collision INTEGER:: ij, k @@ -403,11 +411,9 @@ MODULE moduleMesh DO n = 1, self%nColl !Select random numbers - CALL RANDOM_NUMBER(rndReal) - rnd = 1 + FLOOR(nPart*rndReal) + rnd = random(1, nPart) part_i => partTemp(rnd)%part - CALL RANDOM_NUMBER(rndReal) - rnd = 1 + FLOOR(nPart*rndReal) + rnd = random(1, nPart) part_j => partTemp(rnd)%part ij = interactionIndex(part_i%sp, part_j%sp) sigmaVrelMaxNew = 0.D0 diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index 7816075..a53c7c3 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -148,6 +148,7 @@ MODULE moduleCollisions part_i, part_j) USE moduleSpecies USE moduleConstParam + USE moduleRandom IMPLICIT NONE CLASS(collisionBinaryElastic), INTENT(in):: self @@ -160,27 +161,23 @@ MODULE moduleCollisions REAL(8):: vp_i, vp_j, v_i, v_j !post and pre-collision velocities REAL(8):: v_ij !sum of velocities modules REAL(8):: alpha !random angle of scattering - REAL(8):: rnd !eRel (in units of [m][L]^2[t]^-2 vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1 eRel = self%rMass*vRel**2 sigmaVrel = self%crossSec%get(eRel)*vRel sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel - CALL RANDOM_NUMBER(rnd) - IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN + IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN !Applies the collision v_i = NORM2(part_i%v) v_j = NORM2(part_j%v) v_ij = v_i+v_j vp_j = v_ij*self%w_i vp_i = v_ij*self%w_j - CALL RANDOM_NUMBER(rnd) - alpha = PI*rnd + alpha = PI*random() part_i%v(1) = vp_i*DCOS(alpha) part_i%v(2) = vp_i*DSIN(alpha) - CALL RANDOM_NUMBER(rnd) - alpha = PI*rnd + alpha = PI*random() part_j%v(1) = vp_j*DCOS(alpha) part_j%v(2) = vp_j*DSIN(alpha) @@ -243,8 +240,7 @@ MODULE moduleCollisions USE moduleSpecies USE moduleErrors USE moduleList - USE moduleRefParam !TODO: DELETE - USE moduleConstParam !TODO: DELETE + USE moduleRandom USE OMP_LIB IMPLICIT NONE @@ -256,7 +252,6 @@ MODULE moduleCollisions TYPE(particle), POINTER:: newElectron REAL(8):: vRel, eRel REAL(8):: sigmaVrel - REAL(8):: rnd REAL(8), DIMENSION(1:3):: vp_e, vp_n !eRel (in units of [m][L]^2[t]^-2 @@ -266,8 +261,7 @@ MODULE moduleCollisions IF (eRel > self%eThreshold) THEN sigmaVrel = self%crossSec%get(eRel)*vRel sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel - CALL RANDOM_NUMBER(rnd) - IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN + IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN !Find which particle is the ionizing electron IF (part_i%sp == self%electron%sp) THEN electron => part_i @@ -350,6 +344,7 @@ MODULE moduleCollisions SUBROUTINE collideBinaryChargeExchange(self, sigmaVrelMax, sigmaVrelMaxNew, & part_i, part_j) USE moduleSpecies + USE moduleRandom IMPLICIT NONE CLASS(collisionBinaryChargeExchange), INTENT(in):: self @@ -359,15 +354,13 @@ MODULE moduleCollisions REAL(8):: sigmaVrel REAL(8):: vRel !relative velocity REAL(8):: eRel !relative energy - REAL(8):: rnd !eRel (in units of [m][L]^2[t]^-2 vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1 eRel = self%rMass*vRel**2 sigmaVrel = self%crossSec%get(eRel)*vRel sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel - CALL RANDOM_NUMBER(rnd) - IF (sigmaVrelMaxNew/sigmaVrelMax > rnd) THEN + IF (sigmaVrelMaxNew/sigmaVrelMax > random()) THEN SELECT TYPE(sp => species(part_i%sp)%obj) TYPE IS (speciesNeutral) !Species i is neutral, ionize particle i diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index d222cac..f239fbc 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -174,21 +174,14 @@ MODULE moduleInject !Random velocity from Maxwellian distribution FUNCTION randomVelMaxwellian(self) RESULT (v) - USE moduleConstParam, ONLY: PI + USE moduleRandom IMPLICIT NONE CLASS(velDistMaxwellian), INTENT(in):: self REAL(8):: v - REAL(8):: x, y v = 0.D0 - x = 0.D0 - DO WHILE (x == 0.D0) - CALL RANDOM_NUMBER(x) - END DO - CALL RANDOM_NUMBER(y) - - v = self%v + self%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + v = self%v + self%vTh*randomMaxwellian() END FUNCTION randomVelMaxwellian @@ -208,6 +201,7 @@ MODULE moduleInject USE moduleSpecies USE moduleSolver USE moduleMesh + USE moduleRandom IMPLICIT NONE CLASS(injectGeneric), INTENT(in):: self @@ -215,7 +209,6 @@ MODULE moduleInject INTEGER:: i!, j INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array INTEGER:: n - REAL(8):: rnd CLASS(meshEdge), POINTER:: randomEdge !Insert particles @@ -243,8 +236,7 @@ MODULE moduleInject !$OMP DO DO n = nMin, nMax - CALL RANDOM_NUMBER(rnd) - randomX = INT(DBLE(self%nEdges-1)*rnd) + 1 + randomX = random(1, self%nEdges) randomEdge => mesh%edges(self%edges(randomX))%obj !Random position in edge diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 4fd7188..2243ad5 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -21,42 +21,67 @@ MODULE moduleInput !Loads the config file CALL verboseError('Loading input file...') CALL config%load(filename = inputFile) + CALL checkStatus(config, "load") !Reads reference parameters + CALL verboseError('Reading Reference parameters...') CALL readReference(config) + CALL checkStatus(config, "readReference") !Reads output parameters + CALL verboseError('Reading Output parameters...') CALL readOutput(config) + CALL checkStatus(config, "readOutput") !Read species CALL verboseError('Reading species information...') CALL readSpecies(config) + CALL checkStatus(config, "readSpecies") !Read interactions between species CALL verboseError('Reading interaction between species...') CALL readInteractions(config) + CALL checkStatus(config, "readInteractions") !Read boundaries CALL verboseError('Reading boundary conditions...') CALL readBoundary(config) + CALL checkStatus(config, "readBoundary") !Read Geometry CALL verboseError('Reading Geometry...') CALL readGeometry(config) + CALL checkStatus(config, "readGeometry") !Reads case parameters - CALL verboseError('Reading Case Parameters...') + CALL verboseError('Reading Case parameters...') CALL readCase(config) + CALL checkStatus(config, "readCase") !Read injection of particles CALL verboseError('Reading Interactions between species...') CALL readInject(config) + CALL checkStatus(config, "readInject") !Read parallel parameters CALL verboseError('Reading Parallel configuration...') CALL readParallel(config) + CALL checkStatus(config, "readParallel") END SUBROUTINE readConfig + + !Checks the status of the JSON case file and, if failed, exits the execution. + SUBROUTINE checkStatus(config, step) + USE moduleErrors + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + CHARACTER(LEN=*), INTENT(in):: step + + IF (config%failed()) CALL criticalError("Error reading the JSON input file", TRIM(step)) + + END SUBROUTINE checkStatus !Reads the reference parameters SUBROUTINE readReference(config) @@ -203,8 +228,127 @@ MODULE moduleInput WRITE(tString, '(I1)') iterationDigits iterationFormat = "(I" // tString // "." // tString // ")" + !Read initial state for species + CALL verboseError('Reading Initial state...') + CALL readInitial(config) + CALL checkStatus(config, "readInitial") + END SUBROUTINE readCase + !Reads the initial information for the species + SUBROUTINE readInitial(config) + USE moduleSpecies + USE moduleMesh + USE moduleOutput + USE moduleRefParam + USE moduleRandom + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + LOGICAL:: found + CHARACTER(:), ALLOCATABLE:: object + INTEGER:: nInitial + INTEGER:: i, p, e + CHARACTER(LEN=2):: iString + CHARACTER(:), ALLOCATABLE:: spName + INTEGER:: sp + CHARACTER(:), ALLOCATABLE:: spFile + INTEGER:: stat + CHARACTER(100):: dummy + REAL(8):: density, velocity(1:3), temperature + INTEGER:: nNewPart = 0.D0 + TYPE(particle), POINTER:: partNew + REAL(8):: vTh + TYPE(lNode), POINTER:: partCurr, partNext + + CALL config%info('case.initial', found, n_children = nInitial) + + IF (found) THEN + !Reads the information from initial species + DO i = 1, nInitial + WRITE(iString, '(I2)') i + object = 'case.initial(' // iString // ')' + CALL config%get(object // '.speciesName', spName, found) + sp = speciesName2Index(spName) + CALL config%get(object // '.initialState', spFile, found) + OPEN (10, FILE = path // spFile, ACTION = 'READ') + DO + READ(10, '(A)', IOSTAT = stat) dummy + !If EoF, exit reading + IF (stat /= 0) EXIT + !If comment, skip + IF (INDEX(dummy,'#') /= 0) CYCLE + !Go up one line + BACKSPACE(10) + !Read information + READ(10, *) e, density, velocity, temperature + !Scale variables + !Particles in cell volume + nNewPart = INT(density * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) + !Non-dimensional velocity + velocity = velocity / v_ref + !Non-dimensional temperature + temperature = temperature / T_ref + !Non-dimensional thermal temperature + vTh = DSQRT(temperature/species(sp)%obj%m) + !Allocate new particles + DO p = 1, nNewPart + ALLOCATE(partNew) + partNew%sp = sp + partNew%v(1) = velocity(1) + vTh*randomMaxwellian() + partNew%v(2) = velocity(2) + vTh*randomMaxwellian() + partNew%v(3) = velocity(3) + vTh*randomMaxwellian() + partNew%vol = e + partNew%r = mesh%vols(e)%obj%randPos() + partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r) + partNew%n_in = .TRUE. + partNew%weight = species(sp)%obj%weight + !If charged species, add qm to particle + SELECT TYPE(sp => species(sp)%obj) + TYPE IS (speciesCharged) + partNew%qm = sp%qm + + CLASS DEFAULT + partNew%qm = 0.D0 + + END SELECT + + !Assign particle to temporal list of particles + CALL partInitial%add(partNew) + + END DO + + END DO + + END DO + + !Convert temporal list of particles into initial partOld array + !Deallocates the list of initial particles + nNewPart = partInitial%amount + IF (nNewPart > 0) THEN + ALLOCATE(partOld(1:nNewPart)) + partCurr => partInitial%head + DO p = 1, nNewPart + partNext => partCurr%next + partOld(p) = partCurr%part + DEALLOCATE(partCurr) + partCurr => partNext + + END DO + + IF (ASSOCIATED(partInitial%head)) NULLIFY(partInitial%head) + IF (ASSOCIATED(partInitial%tail)) NULLIFY(partInitial%tail) + partInitial%amount = 0 + + END IF + + nPartOld = SIZE(partOld) + + END IF + + END SUBROUTINE readInitial + !Reads configuration for the output files SUBROUTINE readOutput(config) USE moduleErrors diff --git a/src/modules/moduleList.f90 b/src/modules/moduleList.f90 index 9ef1316..15c4b23 100644 --- a/src/modules/moduleList.f90 +++ b/src/modules/moduleList.f90 @@ -24,6 +24,7 @@ MODULE moduleList INTEGER(KIND=OMP_LOCK_KIND):: lockWScheme !Lock for the NA list of particles TYPE(listNode):: partCollisions !Particles created in collisional process INTEGER(KIND=OMP_LOCK_KIND):: lockCollisions !Lock for the NA list of particles + TYPE(listNode):: partInitial !Initial distribution of particles TYPE pointerArray TYPE(particle), POINTER:: part diff --git a/src/modules/moduleRandom.f90 b/src/modules/moduleRandom.f90 new file mode 100644 index 0000000..ec5eddd --- /dev/null +++ b/src/modules/moduleRandom.f90 @@ -0,0 +1,69 @@ +MODULE moduleRandom + !Interface for random number generator + INTERFACE random + MODULE PROCEDURE randomReal, randomRealAB, randomIntAB + + END INTERFACE random + + CONTAINS + !Returns a Real random number between 0 and 1 + FUNCTION randomReal() RESULT(rnd) + IMPLICIT NONE + + REAL(8):: rnd + + rnd = 0.D0 + CALL RANDOM_NUMBER(rnd) + + END FUNCTION randomReal + + !Returns a Real random number between a and b + FUNCTION randomRealAB(a, b) RESULT(rnd) + IMPLICIT NONE + + REAL(8), INTENT(in):: a, b + REAL(8):: rnd + REAL(8):: rnd01 !random real between 0 and 1 + + rnd = 0.D0 + CALL RANDOM_NUMBER(rnd01) + + rnd = (b - a) * rnd01 + a + + END FUNCTION randomRealAB + + !Returns an Integer random numnber between a and b + FUNCTION randomIntAB(a, b) RESULT(rnd) + IMPLICIT NONE + + INTEGER, INTENT(in):: a, b + INTEGER:: rnd + REAL(8):: rnd01 + + rnd = 0.D0 + CALL RANDOM_NUMBER(rnd01) + + rnd = INT(REAL(b - a) * rnd01) + 1 + + END FUNCTION randomIntAB + + !Returns a random number in a Maxwellian distribution of mean 0 and width 1 + FUNCTION randomMaxwellian() RESULT(rnd) + USE moduleConstParam, ONLY: PI + IMPLICIT NONE + + REAL(8):: rnd + REAL(8):: x, y + + 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) + + END FUNCTION randomMaxwellian + +END MODULE moduleRandom diff --git a/src/modules/moduleTable.f90 b/src/modules/moduleTable.f90 index 39097c0..e4da335 100644 --- a/src/modules/moduleTable.f90 +++ b/src/modules/moduleTable.f90 @@ -52,6 +52,7 @@ MODULE moduleTable i = 0 DO READ(id, '(A)', iostat = stat) dummy + !TODO: Make this a function IF (INDEX(dummy,'#') /= 0) CYCLE IF (stat /= 0) EXIT !Add data