diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json index e0d30dd..bf48e8e 100644 --- a/runs/cylFlow/input.json +++ b/runs/cylFlow/input.json @@ -3,12 +3,12 @@ "path": "./runs/cylFlow/", "trigger": 10, "cpuTime": true, - "numColl": true + "numColl": false }, "geometry": { "type": "2DCyl", "meshType": "gmsh", - "meshFile": "meshTria.msh" + "meshFile": "mesh.msh" }, "species": [ {"name": "Argon", "type": "neutral", "mass": 6.633e-26, "weight": 5.0e7} @@ -30,9 +30,9 @@ "radius": 1.88e-10 }, "case": { - "tau": 6.0e-7, - "time": 3.0e-3, - "solver": "2DCylNeutral" + "tau": [6.0e-7], + "time": 3.0e-3, + "pusher": ["2DCylNeutral"] }, "interactions": { "folderCollisions": "./data/collisions/", diff --git a/src/modules/moduleEM.f95 b/src/modules/moduleEM.f95 index 141b1a8..d013064 100644 --- a/src/modules/moduleEM.f95 +++ b/src/modules/moduleEM.f95 @@ -58,7 +58,7 @@ MODULE moduleEM xi = part%xi - elField = mesh%vols(part%e_p)%obj%gatherEF(xi) + elField = mesh%vols(part%vol)%obj%gatherEF(xi) END FUNCTION gatherElecField diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95 index 2fe2e6f..90ad22f 100644 --- a/src/modules/moduleInject.f95 +++ b/src/modules/moduleInject.f95 @@ -10,7 +10,7 @@ MODULE moduleInject REAL(8):: vTh(1:3) !Thermal Velocity REAL(8):: n(1:3) !Direction of injection INTEGER:: nParticles !Number of particles to introduce each time step - INTEGER:: pt !Species of injection + INTEGER:: sp !Species of injection INTEGER:: nEdges INTEGER, ALLOCATABLE:: edges(:) !Array with edges REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection @@ -26,7 +26,7 @@ MODULE moduleInject CONTAINS !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, T, flow, units, pt, physicalSurface) + SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface) USE moduleMesh USE moduleMeshCyl USE moduleRefParam @@ -38,7 +38,7 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(inout):: self INTEGER, INTENT(in):: i REAL(8), INTENT(in):: v, n(1:3), T(1:3) - INTEGER, INTENT(in):: pt, physicalSurface + INTEGER, INTENT(in):: sp, physicalSurface REAL(8), INTENT(in):: flow CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et @@ -51,21 +51,21 @@ MODULE moduleInject SELECT CASE(units) CASE ("sccm") !Standard cubic centimeter per minute - self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight) + self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(sp)%obj%weight) CASE ("A") !Input current in Ampers - self%nParticles = INT(flow*tau*ti_ref/(qe*species(pt)%obj%weight)) + self%nParticles = INT(flow*tau*ti_ref/(qe*species(sp)%obj%weight)) CASE DEFAULT CALL criticalError("No support for units: " // units, 'initInject') END SELECT IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject') - self%pt = pt + self%sp = sp self%v = self%v * self%n - self%vTh = DSQRT(self%T/species(self%pt)%obj%m) + self%vTh = DSQRT(self%T/species(self%sp)%obj%m) !Gets the edge elements from which particles are injected !TODO: Improve this A LOT @@ -145,13 +145,13 @@ MODULE moduleInject nMin = SUM(inject(1:(self%id-1))%nParticles) + 1 nMax = nMin + self%nParticles - 1 !Assign particle type - partInj(nMin:nMax)%pt = self%pt + partInj(nMin:nMax)%sp = self%sp !Assign weight to particle. - partInj(nMin:nMax)%weight = species(self%pt)%obj%weight + partInj(nMin:nMax)%weight = species(self%sp)%obj%weight !Particle is considered to be outside the domain partInj(nMin:nMax)%n_in = .FALSE. !Assign charge/mass to charged particle. - SELECT TYPE(sp => species(self%pt)%obj) + SELECT TYPE(sp => species(self%sp)%obj) TYPE IS(speciesCharged) partInj(nMin:nMax)%qm = sp%qm @@ -173,10 +173,10 @@ MODULE moduleInject partInj(n)%r = randomEdge%randPos() !Volume associated to the edge: IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN - partInj(n)%e_p = randomEdge%e1%n + partInj(n)%vol = randomEdge%e1%n ELSE - partInj(n)%e_p = randomEdge%e2%n + partInj(n)%vol = randomEdge%e2%n END IF @@ -185,9 +185,9 @@ MODULE moduleInject vBC(self%v(3), self%vTh(3)) /) !Push new particle - CALL solver%push(partInj(n)) + CALL solver%pusher(self%sp)%pushParticle(partInj(n)) !Assign cell to new particle - CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n)) + CALL solver%updateParticleCell(partInj(n)) END DO !$OMP END DO diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95 index 7c24466..5831c7a 100644 --- a/src/modules/moduleInput.f95 +++ b/src/modules/moduleInput.f95 @@ -25,15 +25,15 @@ MODULE moduleInput !Reads reference parameters CALL readReference(config) - !Reads case parameters - CALL readCase(config) - !Reads output parameters CALL readOutput(config) !Read species CALL readSpecies(config) + !Reads case parameters + CALL readCase(config) + !Read interactions between species CALL readInteractions(config) @@ -108,44 +108,69 @@ MODULE moduleInput USE moduleErrors USE moduleCaseParam USE moduleSolver + USE moduleSpecies USE json_module IMPLICIT NONE TYPE(json_file), INTENT(inout):: config LOGICAL:: found + REAL(8), ALLOCATABLE:: tauSpecies(:) CHARACTER(:), ALLOCATABLE:: object REAL(8):: time !simulation time in [t] - CHARACTER(:), ALLOCATABLE:: solverType + CHARACTER(:), ALLOCATABLE:: pusherType, EMType, NAType + INTEGER:: nTau, nSolver + INTEGER:: i + CHARACTER(2):: iString object = 'case' - !Time parameters - CALL config%get(object // '.tau', tau, found) - IF (.NOT. found) CALL criticalError('Required parameter tau not found','readCase') + !Time parameters + CALL config%info(object // '.tau', found, n_children = nTau) + IF (.NOT. found .OR. nTau == 0) CALL criticalError('Required parameter tau not found','readCase') + ALLOCATE(tauSpecies(1:nSpecies)) + DO i = 1, nTau + WRITE(iString, '(I2)') i + CALL config%get(object // '.tau(' // TRIM(iString) // ')', tauSpecies(i), found) + + END DO + IF (nTau < nSpecies) THEN + CALL warningError('Using minimum time step for some species') + tauSpecies(nTau+1:nSpecies) = MINVAL(tauSpecies(1:nTau)) + + END IF + !Selects the minimum tau as reference value + tau = MINVAL(tauSpecies) + + !Gets the simulation time CALL config%get(object // '.time', time, found) IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase') - - CALL config%get(object // '.solver', solverType, found) - IF (.NOT. found) CALL criticalError('Required parameter solver not found','readCase') - - !Allocate the type of solver - SELECT CASE(solverType) - CASE('2DCylNeutral') - !Allocate the solver for neutral pusher 2D Cylindrical - ALLOCATE(solver, source=solverCylNeutral()) - - CASE('2DCylCharged') - !Allocate the solver for charged pusher 2D Cylindrical - ALLOCATE(solver, source=solverCylCharged()) - - CASE DEFAULT - CALL criticalError('Solver ' // solverType // ' not found','readCase') - - END SELECT - !Convert simulation time to number of iterations tmax = INT(time/tau) + !Gest the pusher for each species + CALL config%info(object // '.pusher', found, n_children = nSolver) + IF (.NOT. found .OR. nSolver /= nSpecies) CALL criticalError('Required parameter pusher not found','readCase') + !Allocates all the pushers for particles + ALLOCATE(solver%pusher(1:nSpecies)) + !Initialize pushers + DO i = 1, nSolver + WRITE(iString, '(I2)') i + CALL config%get(object // '.pusher(' // TRIM(iString) // ')', pusherType, found) + + !Associate the type of solver + CALL solver%pusher(i)%init(pusherType, tau, tauSpecies(i)) + + END DO + + !Gets the solver for the electromagnetic field + CALL config%get(object // '.EMSolver', EMType, found) + CALL solver%initEM(EMType) + + !Gest the non-analogue scheme + CALL config%get(object // '.NAScheme', NAType, found) + CALL solver%initNA(NAType) + + !Makes tau non-dimensional tau = tau / ti_ref @@ -234,7 +259,7 @@ MODULE moduleInput !Assign shared parameters for all species CALL config%get(object // '.name', species(i)%obj%name, found) CALL config%get(object // '.weight', species(i)%obj%weight, found) - species(i)%obj%pt = i + species(i)%obj%sp = i species(i)%obj%m = mass END DO @@ -462,7 +487,7 @@ MODULE moduleInput REAL(8):: flow CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface - INTEGER:: pt + INTEGER:: sp CALL config%info('inject', found, n_children = nInject) ALLOCATE(inject(1:nInject)) @@ -473,7 +498,7 @@ MODULE moduleInput !Find species CALL config%get(object // '.species', speciesName, found) - pt = speciesName2Index(speciesName) + sp = speciesName2Index(speciesName) CALL config%get(object // '.name', name, found) CALL config%get(object // '.v', v, found) @@ -483,7 +508,7 @@ MODULE moduleInput CALL config%get(object // '.units', units, found) CALL config%get(object // '.physicalSurface', physicalSurface, found) - CALL inject(i)%init(i, v, normal, T, flow, units, pt, physicalSurface) + CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface) END DO diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95 index d5c464f..a99d096 100644 --- a/src/modules/moduleMeshCyl.f95 +++ b/src/modules/moduleMeshCyl.f95 @@ -463,22 +463,22 @@ MODULE moduleMeshCyl w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%pt) + vertex => self%n1%output(part%sp) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%pt) + vertex => self%n2%output(part%sp) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS - vertex => self%n3%output(part%pt) + vertex => self%n3%output(part%sp) vertex%den = vertex%den + part%weight*w_p(3) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS - vertex => self%n4%output(part%pt) + vertex => self%n4%output(part%sp) vertex%den = vertex%den + part%weight*w_p(4) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS @@ -852,17 +852,17 @@ MODULE moduleMeshCyl w_p = self%weight(part%xi) tensorS = outerProduct(part%v, part%v) - vertex => self%n1%output(part%pt) + vertex => self%n1%output(part%sp) vertex%den = vertex%den + part%weight*w_p(1) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS - vertex => self%n2%output(part%pt) + vertex => self%n2%output(part%sp) vertex%den = vertex%den + part%weight*w_p(2) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS - vertex => self%n3%output(part%pt) + vertex => self%n3%output(part%sp) vertex%den = vertex%den + part%weight*w_p(3) vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS @@ -1041,7 +1041,7 @@ MODULE moduleMeshCyl ! part%weight = part%weight*oldCell%volume/self%volume ! ! END IF - part%e_p = self%n + part%vol = self%n part%xi = xi part%n_in = .TRUE. !Assign particle to listPart_in @@ -1117,7 +1117,7 @@ MODULE moduleMeshCyl part_i => partTemp(rnd)%part rnd = 1 + FLOOR(nPart*RAND()) part_j => partTemp(rnd)%part - ij = interactionIndex(part_i%pt, part_j%pt) + ij = interactionIndex(part_i%sp, part_j%sp) sigmaVrelMaxNew = 0.D0 DO k = 1, interactionMatrix(ij)%amount CALL interactionMatrix(ij)%collisions(k)%obj%collide(self%sigmaVrelMax, sigmaVrelMaxNew, part_i, part_j) diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95 index 21b3766..76e197c 100644 --- a/src/modules/moduleSolver.f95 +++ b/src/modules/moduleSolver.f95 @@ -1,13 +1,31 @@ MODULE moduleSolver - TYPE, PUBLIC, ABSTRACT:: solverGeneric + !Generic type for pusher of particles + TYPE, PUBLIC:: pusherGeneric + PROCEDURE(push_interafece), POINTER, NOPASS:: pushParticle => NULL() + !Boolean to indicate if the species is moved in the iteration + LOGICAL:: pushSpecies + !How many interations between advancing the species + INTEGER:: every CONTAINS - PROCEDURE(push_interafece), DEFERRED, NOPASS:: push - PROCEDURE(solveEM_interface), DEFERRED, NOPASS:: solveEMField + PROCEDURE, PASS:: init => initPusher + + END TYPE pusherGeneric + + !Generic type for solver + TYPE, PUBLIC:: solverGeneric + TYPE(pusherGeneric), ALLOCATABLE:: pusher(:) + PROCEDURE(solveEM_interface), POINTER, NOPASS:: solveEM => NULL() + PROCEDURE(nonAnalogue_interface), POINTER, NOPASS:: nonAnalogue => NULL() + CONTAINS + PROCEDURE, PASS:: initEM + PROCEDURE, PASS:: initNA + PROCEDURE, PASS:: updateParticleCell + PROCEDURE, PASS:: updatePushSpecies END TYPE solverGeneric - ABSTRACT INTERFACE + INTERFACE !Push a particle PURE SUBROUTINE push_interafece(part) USE moduleSpecies @@ -21,38 +39,94 @@ MODULE moduleSolver END SUBROUTINE solveEM_interface + !Apply nonAnalogue scheme to a particle + SUBROUTINE nonAnalogue_interface(part) + USE moduleSpecies + + TYPE(particle), INTENT(inout):: part + + END SUBROUTINE nonAnalogue_interface + END INTERFACE - TYPE, PUBLIC, EXTENDS(solverGeneric):: solverCylNeutral - CONTAINS - PROCEDURE, NOPASS:: push => pushCylNeutral - PROCEDURE, NOPASS:: solveEMField => noEMField - - END TYPE solverCylNeutral - - TYPE, PUBLIC, EXTENDS(solverGeneric):: solverCylCharged - CONTAINS - PROCEDURE, NOPASS:: push => pushCylCharged - PROCEDURE, NOPASS:: solveEMField => solveElecField - - END TYPE solverCylCharged - - CLASS(solverGeneric), ALLOCATABLE:: solver + TYPE(solverGeneric):: solver CONTAINS + !Init Pusher + SUBROUTINE initPusher(self, pusherType, tau, tauSp) + USE moduleErrors + IMPLICIT NONE + + CLASS(pusherGeneric), INTENT(out):: self + CHARACTER(:), ALLOCATABLE:: pusherType + REAL(8):: tau, tauSp + + SELECT CASE(pusherType) + CASE('2DCylNeutral') + self%pushParticle => pushCylNeutral + + CASE('2DCylCharged') + self%pushParticle => pushCylCharged + + CASE DEFAULT + CALL criticalError('Solver ' // pusherType // ' not found','readCase') + + END SELECT + + self%pushSpecies = .FALSE. + self%every = INT(tau/tauSp) + + END SUBROUTINE initPusher + + SUBROUTINE initEM(self, EMType) + IMPLICIT NONE + + CLASS(solverGeneric), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE:: EMType + + SELECT CASE(EMType) + CASE('Electrostatic') + self%solveEM => solveElecField + + CASE DEFAULT + self%solveEM => noEMField + + END SELECT + + END SUBROUTINE initEM + + !Initialize the non-analogue scheme + SUBROUTINE initNA(self, NAType) + IMPLICIT NONE + + CLASS(solverGeneric), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE:: NAType + + SELECT CASE(NAType) + CASE DEFAULT + self%nonAnalogue => noNAScheme + + END SELECT + + END SUBROUTINE initNA + + !Do all pushes for particles SUBROUTINE doPushes() USE moduleSpecies USE moduleMesh IMPLICIT NONE INTEGER:: n + INTEGER:: sp !$OMP DO DO n=1, nPartOld + !Select species type + sp = partOld(n)%sp !Push particle - CALL solver%push(partOld(n)) + CALL solver%pusher(sp)%pushParticle(partOld(n)) !Find cell in wich particle reside - CALL mesh%vols(partOld(n)%e_p)%obj%findCell(partOld(n)) + CALL solver%updateParticleCell(partOld(n)) END DO !$OMP END DO @@ -89,8 +163,8 @@ MODULE moduleSolver END IF part_temp%v(2) = cos_alpha*v_p_oh_star(2)+sin_alpha*v_p_oh_star(3) part_temp%v(3) = -sin_alpha*v_p_oh_star(2)+cos_alpha*v_p_oh_star(3) - part_temp%pt = part%pt - part_temp%e_p = part%e_p + part_temp%sp = part%sp + part_temp%vol = part%vol part_temp%n_in = .FALSE. !Assume particle is outside until cell is found !Copy temporal particle to particle part=part_temp @@ -224,7 +298,7 @@ MODULE moduleSolver !Loops over the particles to scatter them !$OMP DO DO n=1, nPartOld - CALL mesh%vols(partOld(n)%e_p)%obj%scatter(partOld(n)) + CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n)) END DO !$OMP END DO @@ -234,7 +308,7 @@ MODULE moduleSolver SUBROUTINE doEMField() IMPLICIT NONE - CALL solver%solveEMField() + CALL solver%solveEM() END SUBROUTINE doEMField @@ -283,12 +357,59 @@ MODULE moduleSolver END SUBROUTINE solveElecField - !Empty module that does no computation of EM field for neutral cases + !Empty procedure that does no computation of EM field for neutral cases SUBROUTINE noEMField() IMPLICIT NONE END SUBROUTINE noEMField + !Empty procedure that does no computation of EM field for neutral cases + SUBROUTINE noNAScheme(part) + USE moduleSpecies + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + + END SUBROUTINE noNAScheme + + SUBROUTINE updateParticleCell(self, part) + USE moduleSpecies + USE moduleMesh + IMPLICIT NONE + + CLASS(solverGeneric), INTENT(in):: self + TYPE(particle), INTENT(inout):: part + CLASS(meshVol), POINTER:: vol + INTEGER:: volOld !Old cell for particle + + volOld = part%vol + vol => mesh%vols(volOld)%obj + CALL vol%findCell(part) + !If particle has change cell, call nonAnalogue scheme + IF (volOld /= part%vol) THEN + CALL self%nonAnalogue(part) + + END IF + + END SUBROUTINE updateParticleCell + + !Update the information about if a species needs to be moved this iteration + SUBROUTINE updatePushSpecies(self, t) + USE moduleSpecies + IMPLICIT NONE + + CLASS(solverGeneric), INTENT(inout):: self + INTEGER, INTENT(in):: t + INTEGER:: s + + DO s=1, nSpecies + self%pusher(s)%pushSpecies = MOD(t, self%pusher(s)%every) == 0 + + END DO + + END SUBROUTINE updatePushSpecies + + !Output the different data and information SUBROUTINE doOutput(t) USE moduleMesh USE moduleOutput diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95 index 8ba5509..0e02dff 100644 --- a/src/modules/moduleSpecies.f95 +++ b/src/modules/moduleSpecies.f95 @@ -6,7 +6,7 @@ MODULE moduleSpecies TYPE, ABSTRACT:: speciesGeneric CHARACTER(:), ALLOCATABLE:: name REAL(8):: m=0.D0, weight=0.D0 - INTEGER:: pt=0 + INTEGER:: sp=0 END TYPE speciesGeneric TYPE, EXTENDS(speciesGeneric):: speciesNeutral @@ -29,8 +29,8 @@ MODULE moduleSpecies TYPE particle REAL(8):: r(1:3) !Position REAL(8):: v(1:3) !Velocity - INTEGER:: pt !Particle species id - INTEGER:: e_p !Index of element in which the particle is located + INTEGER:: sp !Particle species id + INTEGER:: vol !Index of element in which the particle is located REAL(8):: xi(1:3) !Logical coordinates of particle in element e_p. LOGICAL:: n_in !Flag that indicates if a particle is in the domain REAL(8):: weight=0.D0 !weight of particle @@ -46,17 +46,17 @@ MODULE moduleSpecies TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles CONTAINS - FUNCTION speciesName2Index(speciesName) RESULT(pt) + FUNCTION speciesName2Index(speciesName) RESULT(sp) IMPLICIT NONE CHARACTER(:), ALLOCATABLE:: speciesName - INTEGER:: pt + INTEGER:: sp INTEGER:: n - pt = 0 + sp = 0 DO n = 1, nSpecies IF (speciesName == species(n)%obj%name) THEN - pt = species(n)%obj%pt + sp = species(n)%obj%sp EXIT END IF END DO