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