From 9af3429395b9cb14c58408d2336f2802a06c8941 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Mon, 24 May 2021 12:37:16 +0200 Subject: [PATCH] Issue with random position in volumes Fixed an issue in which the position in triangular an thetrahedron elements were not correctly being computed. Other minor issues fixed: - Units in input file now do not use '/'. - Collisions accuratly conserve momentum. - Minor improvements in mass calculation in collisions. --- runs/0D_Argon/curve_Temperature.gp | 64 +++++++++++++++++++ src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 2 +- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 3 +- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 10 +-- .../inout/gmsh2/moduleMeshOutputGmsh2.f90 | 4 +- src/modules/mesh/moduleMesh.f90 | 12 +--- src/modules/moduleCollisions.f90 | 35 +++++----- 7 files changed, 95 insertions(+), 35 deletions(-) create mode 100644 runs/0D_Argon/curve_Temperature.gp diff --git a/runs/0D_Argon/curve_Temperature.gp b/runs/0D_Argon/curve_Temperature.gp new file mode 100644 index 0000000..4f3eb76 --- /dev/null +++ b/runs/0D_Argon/curve_Temperature.gp @@ -0,0 +1,64 @@ +reset +set macros +term = 'postscript eps color enhanced' +font = 'font "CMUSerif-Roman,20"' +linew = 'lw 2.0' +times='"{/Symbol \264}"' +size_x=8.5 +size_y=5.2 + +tmar=0.998 +bmar=0.145 + +lmarLabel = 0.070 + +lmar=lmarLabel+0.085 +rmar=0.97 + +space_x=(rmar-lmar) + +xt_off_0=0.5 +yt_off_0=0.6 + +set xtics nomirror offset 0.0,xt_off_0 +set mxtics 2 + +set ytics nomirror offset yt_off_0,0.0 +set mytics 2 + +set pointsize 1.5 + +set style line 1 pt 4 lc rgb "#B50427" #Squares red +set style line 2 pt 6 lc rgb "#3B4CC1" #Circles blue +set style line 3 pt 1 lc rgb "#2CA02C" #Crosses green +set style line 4 pt 2 lc rgb "#FE7F0E" #Exes orange +set style line 5 pt 8 lc rgb "#D6696B" #Triangles light red +set style line 10 lt 1 lw 2.0 lc rgb "black" #Black solid line + +set terminal @term size size_x cm, size_y cm @linew @font + +set output "comp_temp.eps" + +#files +filename1 = "OUTPUT_Argon.dat" +filename2 = "OUTPUT_Argon+.dat" + +set lmargin at screen lmar +set rmargin at screen rmar + +set xrange [-0.01:1.01] +set xtics 0.2 +set xlabel "Time (s)" offset 0.0,1.2 + +set yrange [250:3550] +set format y "%3.0f" +set ytics 500 +set ylabel "Temperature (K)" offset 1.4,0.0 + +set key width 0.5 height 0.1 spacing 1.3 samplen 0.2 box opaque font ",16" +set key at graph 0.95, graph 0.9 right top + +plot filename1 u ($1):($7) t "Ar" ls 1 with lines, \ + filename2 u ($1):($7) t "Ar^{+}" ls 2 with lines, \ + '< paste OUTPUT_Argon.dat OUTPUT_Argon+.dat' u ($1):($7 + $14) t "Sum" ls 3 with lines + diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index ff17688..bc60a69 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -622,7 +622,7 @@ MODULE moduleMesh2DCart REAL(8), ALLOCATABLE:: fPsi(:) xii(1) = random( 0.D0, 1.D0) - xii(2) = random( 0.D0, 1.D0) + xii(2) = random( 0.D0, 1.D0 - xii(1)) xii(3) = 0.D0 fPsi = self%fPsi(xii) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index dcd285d..62933d8 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -643,10 +643,9 @@ MODULE moduleMesh2DCyl REAL(8), ALLOCATABLE:: fPsi(:) xii(1) = random( 0.D0, 1.D0) - xii(2) = random( 0.D0, 1.D0) + xii(2) = random( 0.D0, 1.D0 - xii(1)) xii(3) = 0.D0 - ALLOCATE(fPsi(1:3)) fPsi = self%fPsi(xii) r(1) = DOT_PRODUCT(fPsi, self%z) diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 15aed4d..477373f 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -215,7 +215,9 @@ MODULE moduleMesh3DCart REAL(8):: xii(1:3) REAL(8):: fPsi(1:3) - xii = (/random(), random(), 0.D0 /) + xii(1) = random( 0.D0, 1.D0) + xii(2) = random( 0.D0, 1.D0 - xii(1)) + xii(3) = 0.D0 fPsi = self%fPsi(xii) r = (/DOT_PRODUCT(fPsi, self%x), & @@ -294,9 +296,9 @@ MODULE moduleMesh3DCart 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) = random(0.D0, 1.D0) + xii(1) = random( 0.D0, 1.D0) + xii(2) = random( 0.D0, 1.D0 - xii(1)) + xii(3) = random( 0.D0, 1.D0 - xii(1) - xii(2)) ALLOCATE(fPsi(1:4)) fPsi = self%fPsi(xii) diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 92d4a5b..4159873 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -43,7 +43,7 @@ MODULE moduleMeshOutputGmsh2 WRITE(60, "(A)") '$EndNodeData' WRITE(60, "(A)") '$NodeData' WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"' // species(i)%obj%name // ' velocity (m/s)"' + WRITE(60, "(A)") '"' // species(i)%obj%name // ' velocity (m s^-1)"' WRITE(60, *) 1 WRITE(60, *) time WRITE(60, *) 3 @@ -189,7 +189,7 @@ MODULE moduleMeshOutputGmsh2 WRITE(20, "(A)") '$ElementData' WRITE(20, "(A)") '1' - WRITE(20, "(A)") '"Electric Field (V/m)"' + WRITE(20, "(A)") '"Electric Field (V m^-1)"' WRITE(20, *) 1 WRITE(20, *) time WRITE(20, *) 3 diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 030f725..86343f0 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -186,15 +186,6 @@ MODULE moduleMesh END SUBROUTINE initVol_interface - SUBROUTINE scatter_interface(self, part) - USE moduleSpecies - - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - CLASS(particle), INTENT(in):: part - - END SUBROUTINE scatter_interface - PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF) IMPORT:: meshVol CLASS(meshVol), INTENT(in):: self @@ -574,6 +565,7 @@ MODULE moduleMesh !TODO: try to combine this with the findCell for a regular mesh !Find the volume in which particle reside in the mesh for collisions + !No boundary interaction taken into account SUBROUTINE findCellCollMesh(part) USE moduleSpecies IMPLICIT NONE @@ -669,7 +661,7 @@ MODULE moduleMesh DO e=1, self%numVols vol => self%vols(e)%obj nPart = vol%listPart_in%amount - !Computes iterations if there is more than one particle in the cell + !Calculates number of collisions if there is more than one particle in the cell IF (nPart > 1) THEN !Probability of collision pMax = vol%totalWeight*vol%sigmaVrelMax*tauMin/vol%volume diff --git a/src/modules/moduleCollisions.f90 b/src/modules/moduleCollisions.f90 index a2c48ca..6557fa7 100644 --- a/src/modules/moduleCollisions.f90 +++ b/src/modules/moduleCollisions.f90 @@ -5,7 +5,7 @@ MODULE moduleCollisions !Abstract type for collision between two particles TYPE, ABSTRACT:: collisionBinary REAL(8):: rMass !Reduced mass - REAL(8):: sMass !Summed mass + REAL(8):: sMassInv !Summed mass TYPE(table1D):: crossSec !cross section of collision CONTAINS PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide @@ -169,8 +169,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%sMass = mass_i+mass_j - collision%rMass = (mass_i*mass_j)/collision%sMass + collision%sMassInv = 1.D0/(mass_i+mass_j) + collision%rMass = (mass_i*mass_j)*collision%sMassInv END SUBROUTINE initBinaryElastic @@ -180,6 +180,7 @@ MODULE moduleCollisions USE moduleSpecies USE moduleConstParam USE moduleRandom + USE moduleMath IMPLICIT NONE CLASS(collisionBinaryElastic), INTENT(in):: self @@ -193,7 +194,7 @@ MODULE moduleCollisions REAL(8), DIMENSION(1:3):: vCM REAL(8):: vp(1:3) - vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1 + vRel = NORM2(part_i%v-part_j%v) eRel = self%rMass*vRel**2 sigmaVrel = self%crossSec%get(eRel)*vRel sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel @@ -205,8 +206,8 @@ MODULE moduleCollisions vp = vRel*randomDirectionVHS() !Assign velocities to particles - part_i%v = vCM + m_j*vp/self%sMass - part_j%v = vCM - m_i*vp/self%sMass + part_i%v = vCM + m_j*vp*self%sMassInv + part_j%v = vCM - m_i*vp*self%sMassInv END IF @@ -238,8 +239,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%sMass = mass_i+mass_j - collision%rMass = (mass_i*mass_j)/collision%sMass + collision%sMassInv = 1.D0/(mass_i+mass_j) + collision%rMass = (mass_i*mass_j)*collision%sMassInv !Specific parameters for ionization collision SELECT TYPE(collision) @@ -269,6 +270,7 @@ MODULE moduleCollisions USE moduleErrors USE moduleList USE moduleRandom + USE moduleMath USE OMP_LIB IMPLICIT NONE @@ -283,7 +285,7 @@ MODULE moduleCollisions REAL(8), DIMENSION(1:3):: vp_e, vp_n !eRel (in units of [m][L]^2[t]^-2 - vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1 + vRel = NORM2(part_i%v-part_j%v) eRel = self%rMass*vRel**2 !Relative energy must be higher than threshold IF (eRel > self%eThreshold) THEN @@ -370,8 +372,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%sMass = mass_i+mass_j - collision%rMass = (mass_i*mass_j)/collision%sMass + collision%sMassInv = 1.D0/(mass_i+mass_j) + collision%rMass = (mass_i*mass_j)*collision%sMassInv !Specific parameters for ionization collision SELECT TYPE(collision) @@ -401,7 +403,7 @@ MODULE moduleCollisions USE moduleErrors USE moduleList USE moduleRandom - USE OMP_LIB + USE moduleMath IMPLICIT NONE CLASS(collisionBinaryRecombination), INTENT(in):: self @@ -414,7 +416,7 @@ MODULE moduleCollisions REAL(8), DIMENSION(1:3):: vp_i !eRel (in units of [m][L]^2[t]^-2 - vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1 + vRel = NORM2(part_i%v-part_j%v) eRel = self%rMass*vRel**2 !Relative energy must be higher than threshold sigmaVrel = self%crossSec%get(eRel)*vRel @@ -476,8 +478,8 @@ MODULE moduleCollisions CALL collision%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) !Calculates reduced mass - collision%sMass = mass_i+mass_j - collision%rMass = (mass_i*mass_j)/collision%sMass + collision%sMassInv = 1.D0/(mass_i+mass_j) + collision%rMass = (mass_i*mass_j)/collision%sMassInv END SUBROUTINE initBinaryChargeExchange @@ -485,6 +487,7 @@ MODULE moduleCollisions part_i, part_j) USE moduleSpecies USE moduleRandom + USE moduleMath IMPLICIT NONE CLASS(collisionBinaryChargeExchange), INTENT(in):: self @@ -496,7 +499,7 @@ MODULE moduleCollisions REAL(8):: eRel !relative energy !eRel (in units of [m][L]^2[t]^-2 - vRel = SUM(DABS(part_i%v-part_j%v)) !TODO make function of norm1 + vRel = NORM2(part_i%v-part_j%v) eRel = self%rMass*vRel**2 sigmaVrel = self%crossSec%get(eRel)*vRel sigmaVrelMaxNew = sigmaVrelMaxNew + sigmaVrel