From c82279f5c57673e8f2a855c88e44fc86562979f6 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 15 Nov 2020 21:16:02 +0100 Subject: [PATCH] First version with possibility for charged particles to be included. Now, the solver needs to be an input parameter of the case, to select if it is for charged or neutral particles. Resolution of Poisson equation with Dirichlet boundary conditions is possible. The source vector is the charge density. This resolution is done in two steps to save computational time: 1. When reading the mesh, the PLU factorization of the K matrix is computed. 2. In each iteration, the system K*u = f is solved, in which f is the source vector (charge density) and u is the solution (potential) in each node. No case has been added to the repository. This will be done in next commit. The 'non-analog' scheme has been commented. It still needs to split the particle to avoid 'overweight' particles. --- .gitignore | 2 +- makefile | 12 +- runs/cylFlow/input.json | 13 +- src/{DSMC_Neutrals.f95 => PPartiC.f95} | 29 ++- src/makefile | 13 +- src/modules/makefile | 13 +- src/modules/moduleCompTime.f95 | 1 + src/modules/moduleConstParam.f95 | 5 +- src/modules/moduleEM.f95 | 133 +++++++++++++ src/modules/moduleErrors.f95 | 1 - src/modules/moduleInject.f95 | 68 +++++-- src/modules/moduleInput.f95 | 168 +++++++++++++--- src/modules/moduleMesh.f95 | 64 +++++- src/modules/moduleMeshCyl.f95 | 264 +++++++++++++++---------- src/modules/moduleMeshCylBoundary.f95 | 2 - src/modules/moduleMeshCylRead.f95 | 118 +++++++++-- src/modules/moduleOutput.f95 | 9 + src/modules/moduleRefParam.f95 | 5 +- src/modules/moduleSolver.f95 | 161 +++++++++++++-- src/modules/moduleSpecies.f95 | 5 +- 20 files changed, 859 insertions(+), 227 deletions(-) rename src/{DSMC_Neutrals.f95 => PPartiC.f95} (77%) create mode 100644 src/modules/moduleEM.f95 diff --git a/.gitignore b/.gitignore index 0a9a292..ca181ee 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ -DSMC_Neutrals +PPartiC mod/ obj/ doc/user_manual/ diff --git a/makefile b/makefile index b10a4ac..3dfb28d 100644 --- a/makefile +++ b/makefile @@ -12,20 +12,24 @@ JSONINC := $(JSONDIR)/include INCDIR := INCDIR += $(JSONINC) +VAR := "" + # compile flags -FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -fno-stack-arrays -g +FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -g #Output file -OUTPUT = DSMC_Neutrals +OUTPUT = PPartiC export -$(OUTPUT): src.o +all: $(OUTPUT) + +$(OUTPUT): src.o src.o: @mkdir -p $(MODDIR) @mkdir -p $(OBJDIR) - $(MAKE) -C src all + $(MAKE) -C src $(OUTPUT) clean: rm -f $(OUTPUT) diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json index 04814f7..440fb58 100644 --- a/runs/cylFlow/input.json +++ b/runs/cylFlow/input.json @@ -24,14 +24,15 @@ {"name": "Nozzle", "species": "Argon", "flow": 10.0, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0], "n": [1, 0, 0], "physicalSurface": 1} ], "reference": { - "density": 1.0e20, - "mass": 6.633e-26, - "temperature": 300.0, - "radius": 1.88e-10 + "density": 1.0e20, + "mass": 6.633e-26, + "temperature": 300.0, + "radius": 1.88e-10 }, "case": { - "tau": 6.e-7, - "time": 3.e-3 + "tau": 6.0e-7, + "time": 3.0e-3, + "solver": "2DCylNeutral" }, "interactions": { "folderCollisions": "./data/collisions/", diff --git a/src/DSMC_Neutrals.f95 b/src/PPartiC.f95 similarity index 77% rename from src/DSMC_Neutrals.f95 rename to src/PPartiC.f95 index f71bdbb..e7db5c7 100644 --- a/src/DSMC_Neutrals.f95 +++ b/src/PPartiC.f95 @@ -1,5 +1,5 @@ -! PPartiC neutral DSMC main program. -PROGRAM DSMC_Neutrals +! PPartiC neutral PIC main program. +PROGRAM PPartiC USE moduleInput USE moduleErrors USE OMP_LIB @@ -10,28 +10,32 @@ PROGRAM DSMC_Neutrals USE moduleMesh IMPLICIT NONE - ! t = time step; n = number of particle; e = number of element; i = number of inject + ! t = time step INTEGER:: t = 0 + ! arg1 = Input argument 1 (input file) CHARACTER(200):: arg1 + ! inputFile = path+name of input file CHARACTER(:), ALLOCATABLE:: inputFile - tStep = omp_get_wtime() !Gets the input file CALL getarg(1, arg1) inputFile = TRIM(arg1) !If no input file is declared, a critical error is called - IF (inputFile == "") CALL criticalError("No input file provided", "DSMC_Neutrals") + IF (inputFile == "") CALL criticalError("No input file provided", "PPartiC") !Reads the json configuration file CALL readConfig(inputFile) + CALL verboseError("Calculating initial EM field...") + CALL doEMField() + tStep = omp_get_wtime() - tStep !Output initial state CALL doOutput(t) CALL verboseError('Starting main loop...') - !$OMP PARALLEL PRIVATE(t) DEFAULT(SHARED) + !$OMP PARALLEL DEFAULT(SHARED) DO t = 1, tmax !Insert new particles and push them !$OMP SINGLE @@ -46,18 +50,18 @@ PROGRAM DSMC_Neutrals !$OMP SINGLE tPush = omp_get_wtime() - tPush + !Collisions tColl = omp_get_wtime() !$OMP END SINGLE - !Collisions CALL doCollisions() !$OMP SINGLE tColl = omp_get_wtime() - tColl + !Reset particles tReset = omp_get_wtime() !$OMP END SINGLE - !Reset particles CALL doReset() !$OMP SINGLE @@ -71,6 +75,13 @@ PROGRAM DSMC_Neutrals !$OMP SINGLE tWeight = omp_get_wtime() - tWeight + tEMField = omp_get_wtime() + !$OMP END SINGLE + + CALL doEMField() + + !$OMP SINGLE + tEMField = omp_get_wtime() - tEMField tStep = omp_get_wtime() - tStep !Output data CALL doOutput(t) @@ -79,5 +90,5 @@ PROGRAM DSMC_Neutrals END DO !$OMP END PARALLEL -END PROGRAM DSMC_Neutrals +END PROGRAM PPartiC diff --git a/src/makefile b/src/makefile index b0dfff0..f9b4928 100644 --- a/src/makefile +++ b/src/makefile @@ -1,17 +1,16 @@ -OBJECTS = $(OBJDIR)/$(OUTPUT).o \ - $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \ +OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \ $(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \ $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ - $(OBJDIR)/moduleMeshCyl.o $(OBJDIR)/moduleMeshCylRead.o $(OBJDIR)/moduleMeshCylBoundary.o \ - $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ - $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o - + $(OBJDIR)/moduleMeshCyl.o $(OBJDIR)/moduleMeshCylRead.o $(OBJDIR)/moduleMeshCylBoundary.o \ + $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ + $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ + $(OBJDIR)/moduleEM.o all: $(OUTPUT) - $(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(JSONLIB) $(OUTPUT): modules.o $(OUTPUT).f95 $(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f95 + $(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(OBJDIR)/$(OUTPUT).o $(JSONLIB) -L/usr/local/lib -llapack -lopenblas modules.o: $(MAKE) -C modules all diff --git a/src/modules/makefile b/src/modules/makefile index 0cf8c24..e091015 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -1,9 +1,9 @@ OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o\ - moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\ - moduleMeshCylRead.o moduleOutput.o moduleInput.o \ - moduleSolver.o moduleCollisions.o moduleTable.o \ - moduleParallel.o + moduleMesh.o moduleMeshCyl.o moduleMeshCylBoundary.o\ + moduleMeshCylRead.o moduleOutput.o moduleInput.o \ + moduleSolver.o moduleCollisions.o moduleTable.o \ + moduleParallel.o moduleEM.o all: $(OBJS) @@ -40,11 +40,14 @@ moduleRefParam.o: moduleConstParam.o moduleRefParam.f95 moduleSpecies.o: moduleCaseParam.o moduleSpecies.f95 $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ -moduleSolver.o: moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95 +moduleSolver.o: moduleEM.o moduleSpecies.o moduleRefParam.o moduleMesh.o moduleOutput.o moduleSolver.f95 $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ moduleTable.o: moduleErrors.o moduleTable.f95 $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleEM.o: moduleSpecies.o moduleMesh.o moduleEM.f95 + $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ + %.o: %.f95 $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/moduleCompTime.f95 b/src/modules/moduleCompTime.f95 index a45a07d..28273ce 100644 --- a/src/modules/moduleCompTime.f95 +++ b/src/modules/moduleCompTime.f95 @@ -5,6 +5,7 @@ MODULE moduleCompTime REAL(8):: tReset=0.D0 REAL(8):: tColl=0.D0 REAL(8):: tWeight=0.D0 + REAL(8):: tEMField=0.D0 END MODULE moduleCompTime diff --git a/src/modules/moduleConstParam.f95 b/src/modules/moduleConstParam.f95 index a335a53..634f845 100644 --- a/src/modules/moduleConstParam.f95 +++ b/src/modules/moduleConstParam.f95 @@ -6,9 +6,10 @@ MODULE moduleConstParam REAL(8), PARAMETER:: PI = 4.D0*DATAN(1.D0) !number pi REAL(8), PARAMETER:: sccm2atomPerS = 4.5D17 !sccm to atom s^-1 - REAL(8), PARAMETER:: eV2J = 1.60217662D-19 !Electron volt to Joule (elementary charge) + REAL(8), PARAMETER:: qe = 1.60217662D-19 !Elementary charge REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI - ! REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0 (Not used in Neutrals) + REAL(8), PARAMETER:: eV2J = qe !Electron volt to Joule conversion + REAL(8), PARAMETER:: eps_0 = 8.8542D-12 !Epsilon_0 END MODULE moduleConstParam diff --git a/src/modules/moduleEM.f95 b/src/modules/moduleEM.f95 new file mode 100644 index 0000000..996bd33 --- /dev/null +++ b/src/modules/moduleEM.f95 @@ -0,0 +1,133 @@ +!Module to solve the electromagnetic field +MODULE moduleEM + IMPLICIT NONE + + TYPE:: boundaryEM + CHARACTER(:), ALLOCATABLE:: typeEM + INTEGER:: physicalSurface + REAL(8):: potential + + CONTAINS + PROCEDURE, PASS:: apply + + END TYPE boundaryEM + + INTEGER:: nBoundaryEM + TYPE(boundaryEM), ALLOCATABLE:: boundEM(:) + + !Information of charge and reference parameters for rho vector + REAL(8), ALLOCATABLE:: qSpecies(:) + + CONTAINS + SUBROUTINE apply(self, edge) + USE moduleMesh + IMPLICIT NONE + + CLASS(boundaryEM), INTENT(in):: self + CLASS(meshEdge):: edge + INTEGER, ALLOCATABLE:: nodes(:) + INTEGER:: n + + nodes = edge%getNodes() + + DO n = 1, SIZE(nodes) + SELECT CASE(self%typeEM) + CASE ("dirichlet") + mesh%K(nodes(n), :) = 0.D0 + mesh%K(nodes(n), nodes(n)) = 1.D0 + + mesh%nodes(nodes(n))%obj%emData%type = self%typeEM + mesh%nodes(nodes(n))%obj%emData%phi = self%potential + + END SELECT + + END DO + + END SUBROUTINE + + PURE FUNCTION gatherElecField(part) RESULT(elField) + USE moduleSpecies + USE moduleMesh + IMPLICIT NONE + + TYPE(particle), INTENT(in):: part + REAl(8):: xi(1:3) !Logical coordinates of particle in element + REAL(8):: elField(1:3) + + elField = 0.D0 + + xi = part%xLog + + elField = mesh%vols(part%e_p)%obj%gatherEF(xi) + + END FUNCTION gatherElecField + + SUBROUTINE assembleSourceVector(vectorF) + USE moduleMesh + USE moduleRefParam + IMPLICIT NONE + + REAL(8), INTENT(out):: vectorF(1:mesh%numNodes) + REAL(8), ALLOCATABLE:: localF(:) + INTEGER, ALLOCATABLE:: nodes(:) + REAL(8), ALLOCATABLE:: rho(:) + INTEGER:: nNodes + INTEGER:: e, i, ni + CLASS(meshNode), POINTER:: node + + !$OMP SINGLE + vectorF = 0.D0 + !$OMP END SINGLE + + !$OMP SINGLE + DO e = 1, mesh%numVols + nodes = mesh%vols(e)%obj%getNodes() + nNodes = SIZE(nodes) + !Calculates charge density (rho) in element nodes + ALLOCATE(rho(1:nNodes)) + rho = 0.D0 + DO i = 1, nNodes + ni = nodes(i) + node => mesh%nodes(ni)%obj + rho(i) = DOT_PRODUCT(qSpecies(:), node%output(:)%den/(vol_ref*node%v*n_ref)) + + END DO + + !If there is charge in the nodes, compute the local F vector + IF (ANY(rho > 0.D0)) THEN + !Calculates local F vector + localF = mesh%vols(e)%obj%elemF(rho) + + !Assign local F to global F + DO i = 1, nNodes + ni = nodes(i) + vectorF(ni) = vectorF(ni) + localF(i) + + END DO + + DEALLOCATE(localF) + + END IF + + DEALLOCATE(nodes, rho) + + END DO + !$OMP END SINGLE + + !Apply boundary conditions + !$OMP DO + DO i = 1, mesh%numNodes + node => mesh%nodes(i)%obj + + SELECT CASE(node%emData%type) + CASE ("dirichlet") + vectorF(i) = node%emData%phi + + END SELECT + + END DO + !$OMP END DO + + END SUBROUTINE assembleSourceVector + +END MODULE moduleEM diff --git a/src/modules/moduleErrors.f95 b/src/modules/moduleErrors.f95 index 1d61142..0a95db8 100644 --- a/src/modules/moduleErrors.f95 +++ b/src/modules/moduleErrors.f95 @@ -37,7 +37,6 @@ MODULE moduleErrors errorMsg = msg WRITE (*, '(A)') errorMsg - WRITE (*, *) END SUBROUTINE verboseError diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95 index 9b87133..2fe2e6f 100644 --- a/src/modules/moduleInject.f95 +++ b/src/modules/moduleInject.f95 @@ -25,24 +25,14 @@ MODULE moduleInject TYPE(injectGeneric), ALLOCATABLE:: inject(:) CONTAINS - !Injection of particles - SUBROUTINE doInjects() - IMPLICIT NONE - - INTEGER:: i - - DO i=1, nInject - CALL inject(i)%addParticles() - END DO - - END SUBROUTINE doInjects - - SUBROUTINE initInject(self, i, v, n, T, flow, pt, physicalSurface) + !Initialize an injection of particles + SUBROUTINE initInject(self, i, v, n, T, flow, units, pt, physicalSurface) USE moduleMesh USE moduleMeshCyl USE moduleRefParam USE moduleConstParam USE moduleSpecies + USE moduleErrors IMPLICIT NONE CLASS(injectGeneric), INTENT(inout):: self @@ -50,6 +40,7 @@ MODULE moduleInject REAL(8), INTENT(in):: v, n(1:3), T(1:3) INTEGER, INTENT(in):: pt, physicalSurface REAL(8), INTENT(in):: flow + CHARACTER(:), ALLOCATABLE, INTENT(in):: units INTEGER:: e, et INTEGER:: phSurface(1:mesh%numEdges) @@ -57,8 +48,21 @@ MODULE moduleInject self%vMod = v/v_ref self%n = n self%T = T/T_ref - self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight) - self%pt = pt + SELECT CASE(units) + CASE ("sccm") + !Standard cubic centimeter per minute + self%nParticles = INT(flow*sccm2atomPerS*tau*ti_ref/species(pt)%obj%weight) + + CASE ("A") + !Input current in Ampers + self%nParticles = INT(flow*tau*ti_ref/(qe*species(pt)%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%v = self%v * self%n self%vTh = DSQRT(self%T/species(self%pt)%obj%m) @@ -92,6 +96,19 @@ MODULE moduleInject END SUBROUTINE + !Injection of particles + SUBROUTINE doInjects() + + IMPLICIT NONE + + INTEGER:: i + + DO i=1, nInject + CALL inject(i)%addParticles() + END DO + + END SUBROUTINE doInjects + !Random velocity from maxwellian distribution REAL(8) FUNCTION vBC(u, vth) USE moduleConstParam, ONLY: PI @@ -119,18 +136,28 @@ MODULE moduleInject CLASS(injectGeneric), INTENT(in):: self REAL(8):: randomX INTEGER:: j - REAL(8):: randomPos - INTEGER:: nMin, nMax !Min and Max index in partInj array + INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array INTEGER:: n CLASS(meshEdge), POINTER:: randomEdge !Insert particles + !$OMP SINGLE nMin = SUM(inject(1:(self%id-1))%nParticles) + 1 - nMax = nMin + self%nParticles + nMax = nMin + self%nParticles - 1 !Assign particle type partInj(nMin:nMax)%pt = self%pt !Assign weight to particle. partInj(nMin:nMax)%weight = species(self%pt)%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) + TYPE IS(speciesCharged) + partInj(nMin:nMax)%qm = sp%qm + + END SELECT + !$OMP END SINGLE + !$OMP DO DO n = nMin, nMax !Select edge randomly from which inject particle @@ -143,8 +170,7 @@ MODULE moduleInject randomEdge => mesh%edges(self%edges(j))%obj !Random position in edge - randomPos = RAND() - partInj(n)%r = randomEdge%randPos(randomPos) + 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 @@ -159,7 +185,7 @@ MODULE moduleInject vBC(self%v(3), self%vTh(3)) /) !Push new particle - CALL push(partInj(n)) + CALL solver%push(partInj(n)) !Assign cell to new particle CALL mesh%vols(partInj(n)%e_p)%obj%findCell(partInj(n)) diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95 index 926147b..e489a75 100644 --- a/src/modules/moduleInput.f95 +++ b/src/modules/moduleInput.f95 @@ -44,6 +44,9 @@ MODULE moduleInput CALL verboseError('Reading Geometry...') CALL readGeometry(config) + !Read boundary for EM field + CALL readEMBoundary(config) + !Read injection of particles CALL verboseError('Reading Interactions between species...') CALL readInject(config) @@ -63,7 +66,7 @@ MODULE moduleInput IMPLICIT NONE TYPE(json_file), INTENT(inout):: config - LOGICAL:: found + LOGICAL:: found, found_r CHARACTER(:), ALLOCATABLE:: object object = 'reference' @@ -77,15 +80,25 @@ MODULE moduleInput CALL config%get(object // '.temperature', T_ref, found) IF (.NOT. found) CALL criticalError('Reference temperature not found','readReference') - CALL config%get(object // '.radius', r_ref, found) - IF (.NOT. found) CALL criticalError('Reference radius not found','readReference') - + CALL config%get(object // '.radius', r_ref, found_r) + !Derived parameters v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity - sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section - L_ref = 1.D0/(sigma_ref*n_ref) !mean free path + !TODO: Make this solver dependent + IF (found_r) THEN + L_ref = 1.D0/(sigma_ref*n_ref) !mean free path + sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section + + ELSE + L_ref = DSQRT(kb*T_ref*eps_0/n_ref)/qe !Debye length + !TODO: Obtain right sigma_ref for PIC case + sigma_ref = PI*(4.D-10)**2 !reference cross section + + END IF ti_ref = L_ref/v_ref !reference time Vol_ref = L_ref**3 !reference volume + EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field + Volt_ref = EF_ref*L_ref !reference voltage END SUBROUTINE readReference @@ -94,6 +107,7 @@ MODULE moduleInput USE moduleRefParam USE moduleErrors USE moduleCaseParam + USE moduleSolver USE json_module IMPLICIT NONE @@ -101,6 +115,7 @@ MODULE moduleInput LOGICAL:: found CHARACTER(:), ALLOCATABLE:: object REAL(8):: time !simulation time in [t] + CHARACTER(:), ALLOCATABLE:: solverType object = 'case' !Time parameters @@ -110,6 +125,21 @@ MODULE moduleInput 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()) + + END SELECT + !Convert simulation time to number of iterations tmax = INT(time/tau) @@ -149,6 +179,7 @@ MODULE moduleInput CALL config%get(object // '.cpuTime', timeOutput, found) CALL config%get(object // '.numColl', collOutput, found) + CALL config%get(object // '.EMField', emOutput, found) END SUBROUTINE readOutput @@ -164,12 +195,12 @@ MODULE moduleInput CHARACTER(2):: iString CHARACTER(:), ALLOCATABLE:: object CHARACTER(:), ALLOCATABLE:: speciesType - REAL(8):: mass + REAL(8):: mass, charge LOGICAL:: found INTEGER:: i !Gets the number of species - CALL config%info('species', n_children = nSpecies) + CALL config%info('species', found, n_children = nSpecies) !Zero species means critical error IF (nSpecies == 0) CALL criticalError("No species found", "configRead") @@ -180,21 +211,28 @@ MODULE moduleInput WRITE(iString, '(I2)') i object = 'species(' // TRIM(iString) // ')' CALL config%get(object // '.type', speciesType, found) + CALL config%get(object // '.mass', mass, found) + mass = mass/m_ref + !Allocate species depending on type and assign specific parameters SELECT CASE(speciesType) CASE ("neutral") ALLOCATE(species(i)%obj, source=speciesNeutral()) - !TODO: move to subroutine - CALL config%get(object // '.name', species(i)%obj%name, found) - CALL config%get(object // '.mass', mass, found) - CALL config%get(object // '.weight', species(i)%obj%weight, found) - species(i)%obj%pt = i - species(i)%obj%m = mass/m_ref + + CASE ("charged") + CALL config%get(object // '.charge', charge, found) + ALLOCATE(species(i)%obj, source=speciesCharged(q = charge, & + qm = charge/mass)) CASE DEFAULT CALL warningError("Species " // speciesType // " not supported yet") END SELECT + !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%m = mass END DO @@ -217,29 +255,30 @@ MODULE moduleInput CHARACTER(:), ALLOCATABLE:: species_i, species_j CHARACTER(:), ALLOCATABLE:: crossSecFile CHARACTER(:), ALLOCATABLE:: crossSecFilePath + LOGICAL:: found INTEGER:: nInteractions, nCollisions INTEGER:: i, k, ij INTEGER:: pt_i, pt_j CALL initInteractionMatrix(interactionMatrix) - CALL config%get('interactions.folderCollisions', pathCollisions) + CALL config%get('interactions.folderCollisions', pathCollisions, found) - CALL config%info('interactions.collisions', n_children = nInteractions) + CALL config%info('interactions.collisions', found, n_children = nInteractions) DO i = 1, nInteractions WRITE(iString, '(I2)') i object = 'interactions.collisions(' // TRIM(iString) // ')' - CALL config%get(object // '.species_i', species_i) + CALL config%get(object // '.species_i', species_i, found) pt_i = speciesName2Index(species_i) - CALL config%get(object // '.species_j', species_j) + CALL config%get(object // '.species_j', species_j, found) pt_j = speciesName2Index(species_j) - CALL config%info(object // '.crossSections', n_children = nCollisions) + CALL config%info(object // '.crossSections', found, n_children = nCollisions) ij = interactionIndex(pt_i,pt_j) CALL interactionMatrix(ij)%init(nCollisions) DO k = 1, nCollisions WRITE (kString, '(I2)') k - CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile) + CALL config%get(object // '.crossSections(' // TRIM(kString)// ')', crossSecFile, found) crossSecFilePath = pathCollisions // crossSecFile CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) @@ -247,7 +286,6 @@ MODULE moduleInput END DO - END SUBROUTINE readInteractions !Reads boundary conditions for the mesh @@ -263,7 +301,7 @@ MODULE moduleInput LOGICAL:: found INTEGER:: i - CALL config%info('boundary', n_children = nBoundary) + CALL config%info('boundary', found, n_children = nBoundary) ALLOCATE(boundary(1:nBoundary)) DO i = 1, nBoundary WRITE(istring, '(i2)') i @@ -322,6 +360,85 @@ MODULE moduleInput END SUBROUTINE readGeometry + SUBROUTINE readEMBoundary(config) + USE moduleMesh + USE moduleOutput + USE moduleErrors + USE moduleEM + USE moduleRefParam + USE moduleSpecies + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + CHARACTER(:), ALLOCATABLE:: object + LOGICAL:: found + CHARACTER(2):: istring + INTEGER:: i, e, s + INTEGER:: info + EXTERNAL:: dgetrf + + CALL config%info('boundaryEM', found, n_children = nBoundaryEM) + + IF (found) ALLOCATE(boundEM(1:nBoundaryEM)) + + DO i = 1, nBoundaryEM + WRITE(istring, '(i2)') i + object = 'boundaryEM(' // trim(istring) // ')' + + CALL config%get(object // '.type', boundEM(i)%typeEM, found) + + SELECT CASE(boundEM(i)%typeEM) + CASE ("dirichlet") + CALL config%get(object // '.potential', boundEM(i)%potential, found) + IF (.NOT. found) & + CALL criticalError('Required parameter "potential" for Dirichlet boundary condition not found', 'readEMBoundary') + boundEM(i)%potential = boundEM(i)%potential/Volt_ref + + CALL config%get(object // '.physicalSurface', boundEM(i)%physicalSurface, found) + IF (.NOT. found) & + CALL criticalError('Required parameter "physicalSurface" for Dirichlet boundary condition not found', 'readEMBoundary') + + CASE DEFAULT + CALL criticalError('Boundary type ' // boundEM(i)%typeEM // ' not yet supported', 'readEMBoundary') + + END SELECT + + END DO + + ALLOCATE(qSpecies(1:nSpecies)) + DO s = 1, nSpecies + SELECT TYPE(sp => species(s)%obj) + TYPE IS (speciesCharged) + qSpecies(s) = sp%q + + CLASS DEFAULT + qSpecies(s) = 0.D0 + + END SELECT + + END DO + + !TODO: Improve this + IF (ALLOCATED(boundEM)) THEN + DO e = 1, mesh%numEdges + IF (ANY(mesh%edges(e)%obj%physicalSurface == boundEM(:)%physicalSurface)) THEN + DO i = 1, nBoundaryEM + IF (mesh%edges(e)%obj%physicalSurface == boundEM(i)%physicalSurface) THEN + CALL boundEM(i)%apply(mesh%edges(e)%obj) + + END IF + END DO + END IF + END DO + END IF + + !Compute the PLU factorization of K once boundary conditions have been read + CALL dgetrf(mesh%numNodes, mesh%numNodes, mesh%K, mesh%numNodes, mesh%IPIV, info) + IF (info /= 0) CALL criticalError('Factorization of K matrix failed', 'readEMBoundary') + + END SUBROUTINE readEMBoundary + !Reads the injection of particles from the boundaries SUBROUTINE readInject(config) USE moduleSpecies @@ -340,10 +457,11 @@ MODULE moduleInput REAL(8):: v REAL(8), ALLOCATABLE:: T(:), normal(:) REAL(8):: flow + CHARACTER(:), ALLOCATABLE:: units INTEGER:: physicalSurface INTEGER:: pt - CALL config%info('inject', n_children = nInject) + CALL config%info('inject', found, n_children = nInject) ALLOCATE(inject(1:nInject)) nPartInj = 0 DO i = 1, nInject @@ -359,11 +477,13 @@ MODULE moduleInput CALL config%get(object // '.T', T, found) CALL config%get(object // '.n', normal, found) CALL config%get(object // '.flow', flow, found) + CALL config%get(object // '.units', units, found) CALL config%get(object // '.physicalSurface', physicalSurface, found) - CALL inject(i)%init(i, v, normal, T, flow, pt, physicalSurface) + CALL inject(i)%init(i, v, normal, T, flow, units, pt, physicalSurface) END DO + PRINT *, inject(:)%nParticles !Allocate array for injected particles IF (nPartInj > 0) THEN diff --git a/src/modules/moduleMesh.f95 b/src/modules/moduleMesh.f95 index 85048ba..0b6e37f 100644 --- a/src/modules/moduleMesh.f95 +++ b/src/modules/moduleMesh.f95 @@ -2,7 +2,6 @@ MODULE moduleMesh USE moduleList USE moduleOutput - USE OMP_LIB IMPLICIT NONE !Parent of Node element @@ -13,6 +12,7 @@ MODULE moduleMesh REAL(8):: v = 0.D0 !Output values TYPE(outputNode), ALLOCATABLE:: output(:) + TYPE(emNode):: emData CONTAINS PROCEDURE(initNode_interface), DEFERRED, PASS:: init PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates @@ -58,15 +58,17 @@ MODULE moduleMesh !id for boundary condition INTEGER:: bt = 0 CONTAINS - PROCEDURE(initEdge_interface), DEFERRED, PASS:: init - PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary - PROCEDURE(randPos_interface), DEFERRED, PASS:: randPos + PROCEDURE(initEdge_interface), DEFERRED, PASS:: init + PROCEDURE(boundary_interface), DEFERRED, PASS:: fBoundary + PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes + PROCEDURE(randPos_interface), DEFERRED, PASS:: randPos END TYPE meshEdge ABSTRACT INTERFACE SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface) IMPORT:: meshEdge + CLASS(meshEdge), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) @@ -84,10 +86,16 @@ MODULE moduleMesh END SUBROUTINE - PURE FUNCTION randPos_interface(self, rnd) RESULT(r) + PURE FUNCTION getNodesEdge_interface(self) RESULT(n) + IMPORT:: meshEdge + CLASS(meshEdge), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + END FUNCTION + + FUNCTION randPos_interface(self) RESULT(r) IMPORT:: meshEdge CLASS(meshEdge), INTENT(in):: self - REAL(8), INTENT(in):: rnd REAL(8):: r(1:3) END FUNCTION randPos_interface @@ -119,6 +127,9 @@ MODULE moduleMesh CONTAINS PROCEDURE(initVol_interface), DEFERRED, PASS:: init 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(collision_interface), DEFERRED, PASS:: collision PROCEDURE(findCell_interface), DEFERRED, PASS:: findCell PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput @@ -143,6 +154,32 @@ MODULE moduleMesh END SUBROUTINE scatter_interface + PURE FUNCTION gatherEF_interface(self, xi) RESULT(EF) + + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: EF(1:3) + + END FUNCTION gatherEF_interface + + PURE FUNCTION getNodesVol_interface(self) RESULT(n) + + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + END FUNCTION getNodesVol_interface + + PURE FUNCTION elemF_interface(self, source) RESULT(localF) + + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + REAL(8), INTENT(in):: source(1:) + REAL(8), ALLOCATABLE:: localF(:) + + END FUNCTION elemF_interface + SUBROUTINE collision_interface(self) IMPORT:: meshVol CLASS(meshVol), INTENT(inout):: self @@ -184,13 +221,14 @@ MODULE moduleMesh TYPE(meshVolCont), ALLOCATABLE:: vols(:) !Global stiffness matrix REAL(8), ALLOCATABLE, DIMENSION(:,:):: K - !Global load vector - REAL(8), ALLOCATABLE, DIMENSION(:):: F + !Permutation matrix for P L U factorization + INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV CONTAINS PROCEDURE(readMesh_interface), PASS, DEFERRED:: readMesh PROCEDURE(printOutput_interface), PASS, DEFERRED:: printOutput PROCEDURE(printColl_interface), PASS, DEFERRED:: printColl + PROCEDURE(printEM_interface), PASS, DEFERRED:: printEM END TYPE meshGeneric @@ -213,7 +251,7 @@ MODULE moduleMesh END SUBROUTINE printOutput_interface - !Prints unmber of collisions + !Prints number of collisions SUBROUTINE printColl_interface(self, t) IMPORT meshGeneric @@ -222,6 +260,14 @@ MODULE moduleMesh END SUBROUTINE printColl_interface + !Prints EM info + SUBROUTINE printEM_interface(self, t) + IMPORT meshGeneric + + CLASS(meshGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t + + END SUBROUTINE printEM_interface END INTERFACE !Generic mesh diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95 index 69c9bc6..ac912c8 100644 --- a/src/modules/moduleMeshCyl.f95 +++ b/src/modules/moduleMeshCyl.f95 @@ -25,8 +25,9 @@ MODULE moduleMeshCyl !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS - PROCEDURE, PASS:: init => initEdgeCyl - PROCEDURE, PASS:: randPos => randPosCyl + PROCEDURE, PASS:: init => initEdgeCyl + PROCEDURE, PASS:: getNodes => getNodesCyl + PROCEDURE, PASS:: randPos => randPosCyl END TYPE meshEdgeCyl @@ -43,22 +44,20 @@ MODULE moduleMeshCyl CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL() !Connectivity to adjacent elements CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() - !Maximum collision rate - !TODO: Move to other more appropiate section - REAL(8):: phi(1:4) = 0.D0 - REAL(8):: Ke(1:4,1:4) = 0.D0 REAL(8):: arNodes(1:4) = 0.D0 CONTAINS PROCEDURE, PASS:: init => initVolQuadCyl - PROCEDURE, PASS:: locKe => localKeQuad + PROCEDURE, PASS:: elemK => elemKQuad + PROCEDURE, PASS:: elemF => elemFQuad PROCEDURE, PASS:: detJac => detJQuad PROCEDURE, PASS:: invJ => invJQuad PROCEDURE, PASS:: area => areaQuad PROCEDURE, NOPASS:: weight => weightQuad PROCEDURE, NOPASS:: inside => insideQuad - PROCEDURE, PASS:: dVal => dValQuad PROCEDURE, PASS:: scatter => scatterQuad + PROCEDURE, PASS:: gatherEF => gatherEFQuad + PROCEDURE, PASS:: getNodes => getNodesQuad PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: findCell => findCellCylQuad PROCEDURE, PASS:: resetOutput => resetOutputQuad @@ -78,8 +77,8 @@ MODULE moduleMeshCyl REAL(8), INTENT(in):: r(1:3) self%n = n - self%r = r(1)/L_ref self%z = r(2)/L_ref + self%r = r(1)/L_ref !Node volume, to be determined in mesh self%v = 0.D0 @@ -128,13 +127,25 @@ MODULE moduleMeshCyl END SUBROUTINE initEdgeCyl + PURE FUNCTION getNodesCyl(self) RESULT(n) + IMPLICIT NONE - PURE FUNCTION randPosCyl(self, rnd) RESULT(r) CLASS(meshEdgeCyl), INTENT(in):: self - REAL(8), INTENT(in):: rnd + INTEGER, ALLOCATABLE:: n(:) + + ALLOCATE(n(1:2)) + n = (/self%n1%n, self%n2%n /) + + END FUNCTION getNodesCyl + + FUNCTION randPosCyl(self) RESULT(r) + CLASS(meshEdgeCyl), INTENT(in):: self + REAL(8):: rnd REAL(8):: r(1:3) REAL(8):: p1(1:2), p2(1:2) + rnd = RAND() + p1 = (/self%z(1), self%r(1) /) p2 = (/self%z(2), self%r(2) /) r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 @@ -173,15 +184,14 @@ MODULE moduleMeshCyl self%n3%v = self%n3%v + self%arNodes(3) self%n4%v = self%n4%v + self%arNodes(4) - CALL self%locKe() - self%sigmaVrelMax = sigma_ref/L_ref**2 CALL OMP_INIT_LOCK(self%lock) END SUBROUTINE initVolQuadCyl - FUNCTION fpsi(xi1,xi2) RESULT(psi) + !TODO: MAKE ELEMENT DEPENDENT + PURE FUNCTION fpsi(xi1,xi2) RESULT(psi) IMPLICIT NONE REAL(8),INTENT(in):: xi1,xi2 @@ -195,7 +205,8 @@ MODULE moduleMeshCyl END FUNCTION fpsi !Derivative element function (xi1) - FUNCTION dpsiXi1(xi2) RESULT(psiXi1) + !TODO: MAKE ELEMENT DEPENDENT + PURE FUNCTION dpsiXi1(xi2) RESULT(psiXi1) IMPLICIT NONE REAL(8),INTENT(in):: xi2 @@ -209,7 +220,8 @@ MODULE moduleMeshCyl END FUNCTION dpsiXi1 !Derivative element function (xi2) - FUNCTION dpsiXi2(xi1) RESULT(psiXi2) + !TODO: MAKE ELEMENT DEPENDENT + PURE FUNCTION dpsiXi2(xi1) RESULT(psiXi2) IMPLICIT NONE REAL(8),INTENT(in):: xi1 @@ -223,12 +235,13 @@ MODULE moduleMeshCyl END FUNCTION dpsiXi2 !Transforms physical coordinates to element coordinates - FUNCTION phy2logQuad(this,r) RESULT(xN) + PURE FUNCTION phy2logQuad(self,r) RESULT(xN) IMPLICIT NONE - CLASS(meshVolCylQuad):: this - REAL(8):: r(1:2) - REAL(8):: xN(1:2), xO(1:2), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2) + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: r(1:3) + REAL(8):: xN(1:3) + REAL(8):: xO(1:3), dPsi(1:2,1:4), detJ, j(1:2,1:2), psi(1:4), f(1:2) REAL(8):: conv !Iterative newton method to transform coordinates @@ -238,12 +251,12 @@ MODULE moduleMeshCyl DO WHILE(conv>1.D-4) dPsi(1,:) = dpsiXi1(xO(2)) dPsi(2,:) = dpsiXi2(xO(1)) - j = this%invJ(xO(1),xO(2),dPsi) - psi = fpsi(xO(1), xO(2)) - f(1) = DOT_PRODUCT(psi,this%z)-r(1) - f(2) = DOT_PRODUCT(psi,this%r)-r(2) - detJ = this%detJac(xO(1),xO(2),dPsi) - xN=xO - MATMUL(j, f)/detJ + j = self%invJ(xO, dPsi) + psi = fpsi(xO(1), xO(2)) + f(1) = DOT_PRODUCT(psi,self%z)-r(1) + f(2) = DOT_PRODUCT(psi,self%r)-r(2) + detJ = self%detJac(xO,dPsi) + xN(1:2)=xO(1:2) - MATMUL(j, f)/detJ conv=MAXVAL(DABS(xN-xO),1) xO=xN @@ -252,106 +265,139 @@ MODULE moduleMeshCyl END FUNCTION phy2logQuad !Computes element local stiffness matrix - SUBROUTINE localKeQuad(this) + PURE FUNCTION elemKQuad(self) RESULT(ke) USE moduleConstParam IMPLICIT NONE - CLASS(meshVolCylQuad):: this - REAL(8):: r, xi1, xi2, dPsi(1:2,1:4) + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8):: r, xi(1:3), dPsi(1:2,1:4) + REAL(8):: ke(1:4,1:4) REAL(8):: j(1:2,1:2), detJ INTEGER:: l, m - this%Ke=0.D0 + ke=0.D0 !Start 2D Gauss Quad Integral DO l=1, 3 + xi(1) = cor(l) + dPsi(2,:) = dpsiXi2(xi(1)) DO m = 1, 3 - xi1 = cor(l) - xi2 = cor(m) - dPsi(1,:) = dpsiXi1(xi2) - dPsi(2,:) = dpsiXi2(xi1) - detJ = this%detJac(xi1,xi2,dPsi) - j = this%invJ(xi1,xi2,dPsi) - r = DOT_PRODUCT(fpsi(xi1,xi2),this%r) - this%Ke = this%Ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ + xi(2) = cor(m) + dPsi(1,:) = dpsiXi1(xi(2)) + detJ = self%detJac(xi,dPsi) + j = self%invJ(xi,dPsi) + r = DOT_PRODUCT(fpsi(xi(1),xi(2)),self%r) + ke = ke + MATMUL(TRANSPOSE(MATMUL(j,dPsi)),MATMUL(j,dPsi))*r*w(l)*w(m)/detJ END DO END DO - this%Ke = this%Ke*2.D0*PI + ke = ke*2.D0*PI - END SUBROUTINE localKeQuad + END FUNCTION elemKQuad - !Computes element Jacobian determinant - FUNCTION detJQuad(this,xi1,xi2,dPsi_in) RESULT(dJ) + !Computes the local source vector for a force f + PURE FUNCTION elemFQuad(self, source) RESULT(localF) + USE moduleConstParam IMPLICIT NONE - REAL(8),OPTIONAL:: dPsi_in(1:2,1:4) + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: source(1:) + REAL(8), ALLOCATABLE:: localF(:) + REAL(8):: r, xi(1:3), psi(1:4) + REAL(8):: detJ, f + INTEGER:: l, m + + ALLOCATE(localF(1:4)) + localF = 0.D0 + xi = 0.D0 + DO l=1, 3 + DO m = 1, 3 + xi(1) = cor(l) + xi(2) = cor(m) + detJ = self%detJac(xi) + psi = fpsi(xi(1),xi(2)) + r = DOT_PRODUCT(psi,self%r) + f = DOT_PRODUCT(psi,source) + localF = localF + r*f*psi*w(l)*w(m)*detJ + + END DO + END DO + localF = localF*2.D0*PI + + END FUNCTION elemFQuad + + !Computes element Jacobian determinant + PURE FUNCTION detJQuad(self,xi,dPsi_in) RESULT(dJ) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:2,1:4) + REAL(8):: dJ REAL(8):: dPsi(1:2,1:4) REAL(8):: dz(1:2), dr(1:2) - REAL(8):: xi1,xi2, dJ - CLASS(meshVolCylQuad):: this IF(PRESENT(dPsi_in)) THEN dPsi=dPsi_in ELSE - dPsi(1,:) = dpsiXi1(xi2) - dPsi(2,:) = dpsiXi2(xi1) + dPsi(1,:) = dpsiXi1(xi(2)) + dPsi(2,:) = dpsiXi2(xi(1)) END IF - dz(1) = DOT_PRODUCT(dPsi(1,:),this%z) - dz(2) = DOT_PRODUCT(dPsi(2,:),this%z) - dr(1) = DOT_PRODUCT(dPsi(1,:),this%r) - dr(2) = DOT_PRODUCT(dPsi(2,:),this%r) + dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) dJ = dz(1)*dr(2)-dz(2)*dr(1) END FUNCTION detJQuad - FUNCTION invJQuad(this,xi1,xi2,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJQuad(self,xi,dPsi_in) RESULT(j) !Computes element Jacobian inverse matrix (without determinant) IMPLICIT NONE - REAL(8),OPTIONAL:: dpsi_in(1:2,1:4) + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in), OPTIONAL:: dpsi_in(1:2,1:4) REAL(8):: dPsi(1:2,1:4) REAL(8):: dz(1:2), dr(1:2) - REAL(8):: xi1,xi2, j(1:2,1:2) - CLASS(meshVolCylQuad):: this + REAL(8):: j(1:2,1:2) IF(PRESENT(dPsi_in)) THEN dPsi=dPsi_in ELSE - dPsi(1,:) = dpsiXi1(xi2) - dPsi(2,:) = dpsiXi2(xi1) + dPsi(1,:) = dpsiXi1(xi(2)) + dPsi(2,:) = dpsiXi2(xi(1)) END IF - dz(1) = DOT_PRODUCT(dPsi(1,:),this%z) - dz(2) = DOT_PRODUCT(dPsi(2,:),this%z) - dr(1) = DOT_PRODUCT(dPsi(1,:),this%r) - dr(2) = DOT_PRODUCT(dPsi(2,:),this%r) + dz(1) = DOT_PRODUCT(dPsi(1,:),self%z) + dz(2) = DOT_PRODUCT(dPsi(2,:),self%z) + dr(1) = DOT_PRODUCT(dPsi(1,:),self%r) + dr(2) = DOT_PRODUCT(dPsi(2,:),self%r) j(1,:) = (/ dr(2), -dz(2) /) j(2,:) = (/ -dr(1), dz(1) /) END FUNCTION invJQuad - SUBROUTINE areaQuad(this)!Computes element area + SUBROUTINE areaQuad(self)!Computes element area USE moduleConstParam IMPLICIT NONE - CLASS(meshVolCylQuad):: this - REAL(8):: r, xi1, xi2 + CLASS(meshVolCylQuad):: self + REAL(8):: r, xi(1:3) REAL(8):: detJ, fpsi0(1:4) - this%volume = 0.D0 - this%arNodes = 0.D0 + self%volume = 0.D0 + self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral - xi1 = 0.D0 - xi2 = 0.D0 - detJ = this%detJac(xi1,xi2)*8.D0*PI !4*2*pi - fpsi0 = fpsi(xi1,xi2) - r = DOT_PRODUCT(fpsi0,this%r) - this%volume = r*detJ - this%arNodes = fpsi0*r*detJ + xi = 0.D0 + detJ = self%detJac(xi)*8.D0*PI !4*2*pi + fpsi0 = fpsi(xi(1),xi(2)) + r = DOT_PRODUCT(fpsi0,self%r) + self%volume = r*detJ + self%arNodes = fpsi0*r*detJ END SUBROUTINE areaQuad @@ -376,24 +422,6 @@ MODULE moduleMeshCyl END FUNCTION insideQuad - FUNCTION dValQuad(this,xi1,xi2) RESULT(EF)!Computes electric field in xi1,xi2 - IMPLICIT NONE - - CLASS(meshVolCylQuad):: this - REAL(8):: xi1, xi2, dPsi(1:2,1:4) - REAL(8):: j(1:2,1:2), detJ - REAL(8):: EF(1:3) - - dPsi(1,:) = dpsiXi1(xi2) - dPsi(2,:) = dpsiXi2(xi1) - detJ = this%detJac(xi1,xi2,dPsi) - j = this%invJ(xi1,xi2,dPsi) - EF(1) = -DOT_PRODUCT(dPsi(1,:), this%phi)*j(2,2)/detJ - EF(2) = -DOT_PRODUCT(dPsi(2,:), this%phi)*j(1,1)/detJ - EF(3) = 0.D0 - - END FUNCTION dValQuad - SUBROUTINE scatterQuad(self, part) USE moduleOutput USE moduleSpecies @@ -428,29 +456,65 @@ MODULE moduleMeshCyl vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:) vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS - END SUBROUTINE scatterQuad + PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF)!Computes electric field inside element + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: dPsi(1:2,1:4) + REAL(8):: j(1:2,1:2), detJ + REAL(8):: phi(1:4) + REAL(8):: EF(1:3) + + phi = (/self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%emData%phi, & + self%n4%emData%phi /) + + dPsi(1,:) = dpsiXi1(xi(2)) + dPsi(2,:) = dpsiXi2(xi(1)) + detJ = self%detJac(xi,dPsi) + j = self%invJ(xi,dPsi) + EF(1) = -DOT_PRODUCT(dPsi(1,:), phi)*j(2,2)/detJ + EF(2) = -DOT_PRODUCT(dPsi(2,:), phi)*j(1,1)/detJ + EF(3) = 0.D0 + + END FUNCTION gatherEFQuad + + PURE FUNCTION getNodesQuad(self) RESULT(n) + IMPLICIT NONE + + CLASS(meshVolCylQuad), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + ALLOCATE(n(1:4)) + n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) + + END FUNCTION getNodesQuad + RECURSIVE SUBROUTINE findCellCylQuad(self, part, oldCell) USE moduleSpecies + USE OMP_LIB IMPLICIT NONE CLASS(meshVolCylQuad), INTENT(inout):: self CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell CLASS(particle), INTENT(inout), TARGET:: part - REAL(8):: xLog(1:2) + REAL(8):: xLog(1:3) REAL(8):: xLogArray(1:4) CLASS(*), POINTER:: nextElement INTEGER:: nextInt - xLog = self%phy2log(part%r(1:2)) + xLog = self%phy2log(part%r) + !Checks if particle is inside of current cell IF (self%inside(xLog(1), xLog(2))) THEN - !Checks if particle is inside of current cell - IF (PRESENT(oldCell)) THEN - !If oldCell, recalculate particle weight, as particle has entered a new cell. - part%weight = part%weight*oldCell%volume/self%volume - - END IF + ! IF (PRESENT(oldCell)) THEN + ! !If oldCell, recalculate particle weight, as particle has entered a new cell. + ! part%weight = part%weight*oldCell%volume/self%volume + ! + ! END IF part%e_p = self%n part%xLog = xLog part%n_in = .TRUE. diff --git a/src/modules/moduleMeshCylBoundary.f95 b/src/modules/moduleMeshCylBoundary.f95 index 4298193..7c29115 100644 --- a/src/modules/moduleMeshCylBoundary.f95 +++ b/src/modules/moduleMeshCylBoundary.f95 @@ -74,8 +74,6 @@ MODULE moduleMeshCylBoundary CLASS(meshEdgeCylAxis), INTENT(inout):: self CLASS(particle), INTENT(inout):: part - !Do to the 2D Cyl pusher, this function should never be called. - END SUBROUTINE symmetryAxis diff --git a/src/modules/moduleMeshCylRead.f95 b/src/modules/moduleMeshCylRead.f95 index df6e51b..9a0f5ff 100644 --- a/src/modules/moduleMeshCylRead.f95 +++ b/src/modules/moduleMeshCylRead.f95 @@ -8,6 +8,7 @@ MODULE moduleMeshCylRead PROCEDURE, PASS:: readMesh => readMeshCyl PROCEDURE, PASS:: printOutput => printOutputCyl PROCEDURE, PASS:: printColl => printCollisionsCyl + PROCEDURE, PASS:: printEM => printEMCyl END TYPE @@ -41,8 +42,10 @@ MODULE moduleMeshCylRead READ(10, *) self%numNodes !Allocate required matrices and vectors ALLOCATE(self%nodes(1:self%numNodes)) - ! ALLOCATE(self%K(1:numNodes,1:numNodes)) - ! ALLOCATE(self%F(1:numNodes)) + ALLOCATE(self%K(1:self%numNodes,1:self%numNodes)) + ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes)) + self%K = 0.D0 + self%IPIV = 0 !Read nodes cartesian coordinates (x=z, y=r, z=null) DO e=1, self%numNodes READ(10, *) n, z, r @@ -123,24 +126,10 @@ MODULE moduleMeshCylRead END DO - END DO + !Constructs the global K matrix + CALL constructGlobalK(self%K, self%vols(e)%obj) - ! !Compute global stiffness matrix - ! GlobalK=0.D0 - ! DO e=1, numElem - ! DO i=1, 4 - ! DO j=1, 4 - ! GlobalK(elems(e)%p(i),elems(e)%p(j)) = GlobalK(elems(e)%p(i),elems(e)%p(j)) + elems(e)%Ke(i,j) - ! END DO - ! END DO - ! END DO - ! ! Apply Dirichlet boundary conditions to GlobalK - ! DO n=1, numNodes - ! IF (nodes(n)%bound == 1) THEN - ! GlobalK(n,:)=0.D0 - ! GlobalK(n,n)=1.D0 - ! END IF - ! END DO + END DO END SUBROUTINE @@ -298,6 +287,41 @@ MODULE moduleMeshCylRead END SUBROUTINE connectedQuadEdge + SUBROUTINE constructGlobalK(K, elem) + IMPLICIT NONE + + REAL(8), INTENT(inout):: K(1:,1:) + CLASS(meshVol), INTENT(in):: elem + REAL(8), ALLOCATABLE:: localK(:,:) + INTEGER:: nNodes, i, j + INTEGER, ALLOCATABLE:: n(:) + + SELECT TYPE(elem) + TYPE IS(meshVolCylQuad) + nNodes = 4 + ALLOCATE(localK(1:nNodes,1:nNodes)) + localK = elem%elemK() + ALLOCATE(n(1:nNodes)) + n = (/ elem%n1%n, elem%n2%n, & + elem%n3%n, elem%n4%n /) + + CLASS DEFAULT + nNodes = 0 + ALLOCATE(localK(1:1, 1:1)) + localK = 0.D0 + ALLOCATE(n(1:1)) + n = 0 + + END SELECT + + DO i = 1, nNodes + DO j = 1, nNodes + K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j) + END DO + END DO + + END SUBROUTINE constructGlobalK + SUBROUTINE printOutputCyl(self, t) USE moduleRefParam USE moduleSpecies @@ -425,4 +449,60 @@ MODULE moduleMeshCylRead END SUBROUTINE printCollisionsCyl + SUBROUTINE printEMCyl(self, t) + USE moduleRefParam + USE moduleCaseParam + USE moduleOutput + IMPLICIT NONE + + CLASS(meshCylGeneric), INTENT(in):: self + INTEGER, INTENT(in):: t + INTEGER:: n, e + REAL(8):: time + CHARACTER(:), ALLOCATABLE:: fileName + CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations + + IF (emOutput) THEN + time = DBLE(t)*tau*ti_ref + WRITE(tstring, '(I6.6)') t + + fileName='OUTPUT_' // tstring// '_EMField.msh' + WRITE(*, "(6X,A15,A)") "Creating file: ", fileName + OPEN (20, file = path // folder // '/' // fileName) + WRITE(20, "(A)") '$MeshFormat' + WRITE(20, "(A)") '2.2 0 8' + WRITE(20, "(A)") '$EndMeshFormat' + WRITE(20, "(A)") '$NodeData' + WRITE(20, "(A)") '1' + WRITE(20, "(A)") '"Potential (V)"' + WRITE(20, *) 1 + WRITE(20, *) time + WRITE(20, *) 3 + WRITE(20, *) t + WRITE(20, *) 1 + WRITE(20, *) self%numNodes + DO n=1, self%numNodes + WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref + END DO + WRITE(20, "(A)") '$EndNodeData' + + WRITE(20, "(A)") '$ElementData' + WRITE(20, "(A)") '1' + WRITE(20, "(A)") '"Electric Field (V/m)"' + WRITE(20, *) 1 + WRITE(20, *) time + WRITE(20, *) 3 + WRITE(20, *) t + WRITE(20, *) 3 + WRITE(20, *) self%numVols + DO e=1, self%numVols + WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF((/0.D0, 0.D0, 0.D0/))*EF_ref + END DO + WRITE(20, "(A)") '$EndElementData' + CLOSE(20) + + END IF + + END SUBROUTINE printEMCyl + END MODULE moduleMeshCylRead diff --git a/src/modules/moduleOutput.f95 b/src/modules/moduleOutput.f95 index dfe3808..fe1df6f 100644 --- a/src/modules/moduleOutput.f95 +++ b/src/modules/moduleOutput.f95 @@ -7,6 +7,13 @@ MODULE moduleOutput END TYPE + !Type for EM data in node + TYPE emNode + CHARACTER(:), ALLOCATABLE:: type + REAL(8):: phi + + END TYPE emNode + !Output in dimensional units to print TYPE outputFormat REAL(8):: density, velocity(1:3), pressure, temperature @@ -18,6 +25,7 @@ MODULE moduleOutput INTEGER:: triggerOutput, counterOutput LOGICAL:: timeOutput = .FALSE. LOGICAL:: collOutput = .FALSE. + LOGICAL:: emOutput = .FALSE. CONTAINS FUNCTION outerProduct(a,b) RESULT(s) @@ -42,6 +50,7 @@ MODULE moduleOutput END FUNCTION tensorTrace SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn) + USE moduleConstParam USE moduleRefParam USE moduleSpecies IMPLICIT NONE diff --git a/src/modules/moduleRefParam.f95 b/src/modules/moduleRefParam.f95 index 152f103..9d2b8d9 100644 --- a/src/modules/moduleRefParam.f95 +++ b/src/modules/moduleRefParam.f95 @@ -1,10 +1,9 @@ !Reference parameters MODULE moduleRefParam - USE moduleConstParam !Parameters that define the problem (inputs) - REAL(8):: n_ref, m_ref, T_ref, r_ref, sigma_ref!, q_ref=1.6022D-19 + REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref !Reference parameters for non-dimensional problem - REAL(8):: L_ref, v_ref, ti_ref, Vol_ref + REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref END MODULE moduleRefParam diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95 index ec05cd7..21b3766 100644 --- a/src/modules/moduleSolver.f95 +++ b/src/modules/moduleSolver.f95 @@ -1,5 +1,44 @@ MODULE moduleSolver + TYPE, PUBLIC, ABSTRACT:: solverGeneric + CONTAINS + PROCEDURE(push_interafece), DEFERRED, NOPASS:: push + PROCEDURE(solveEM_interface), DEFERRED, NOPASS:: solveEMField + + END TYPE solverGeneric + + ABSTRACT INTERFACE + !Push a particle + PURE SUBROUTINE push_interafece(part) + USE moduleSpecies + + TYPE(particle), INTENT(inout):: part + + END SUBROUTINE push_interafece + + !Solve the electromagnetic field + SUBROUTINE solveEM_interface() + + END SUBROUTINE solveEM_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 + CONTAINS SUBROUTINE doPushes() USE moduleSpecies @@ -10,7 +49,9 @@ MODULE moduleSolver !$OMP DO DO n=1, nPartOld - CALL push(partOld(n)) + !Push particle + CALL solver%push(partOld(n)) + !Find cell in wich particle reside CALL mesh%vols(partOld(n)%e_p)%obj%findCell(partOld(n)) END DO @@ -18,16 +59,15 @@ MODULE moduleSolver END SUBROUTINE doPushes - !Push one particle. Boris pusher for 2D Cyl - !TODO: make general for any pusher - PURE SUBROUTINE push(part) + !Push one particle. Boris pusher for 2D Cyl Neutral particle + PURE SUBROUTINE pushCylNeutral(part) USE moduleSpecies - USE moduleMesh IMPLICIT NONE + TYPE(particle), INTENT(inout):: part REAL(8):: v_p_oh_star(2:3) TYPE(particle):: part_temp - REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha, alpha + REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha part_temp = part !z @@ -40,8 +80,6 @@ MODULE moduleSolver y_new = v_p_oh_star(3)*tau r = DSQRT(x_new**2+y_new**2) part_temp%r(2) = r - alpha = 0.D0!ATAN2(y_new,x_new) !0 in 2D problem - part_temp%r(3) = part%r(3) + alpha IF (r > 0.D0) THEN sin_alpha = y_new/r cos_alpha = x_new/r @@ -57,7 +95,47 @@ MODULE moduleSolver !Copy temporal particle to particle part=part_temp - END SUBROUTINE push + END SUBROUTINE pushCylNeutral + + !Push one particle. Boris pusher for 2D Cyl Charged particle + PURE SUBROUTINE pushCylCharged(part) + USE moduleSpecies + USE moduleEM + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + REAL(8):: v_p_oh_star(2:3) + TYPE(particle):: part_temp + REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha + REAL(8):: qmEFt(1:3)!charge*tau*EF/mass + + part_temp = part + !Get electric field at particle position + qmEFt = part_temp%qm*gatherElecField(part_temp)*tau + !z + part_temp%v(1) = part%v(1) + qmEFt(1) + part_temp%r(1) = part%r(1) + part_temp%v(1)*tau + !r,theta + v_p_oh_star(2) = part%v(2) + qmEFt(2) + x_new = part%r(2) + v_p_oh_star(2)*tau + v_p_oh_star(3) = part%v(3) + qmEFt(3) + y_new = v_p_oh_star(3)*tau + r = DSQRT(x_new**2+y_new**2) + part_temp%r(2) = r + IF (r > 0.D0) THEN + sin_alpha = y_new/r + cos_alpha = x_new/r + ELSE + sin_alpha = 0.D0 + cos_alpha = 1.D0 + 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%n_in = .FALSE. !Assume particle is outside until cell is found + !Copy temporal particle to particle + part=part_temp + + END SUBROUTINE pushCylCharged !Do the collisions in all the cells SUBROUTINE doCollisions() @@ -153,6 +231,64 @@ MODULE moduleSolver END SUBROUTINE doScatter + SUBROUTINE doEMField() + IMPLICIT NONE + + CALL solver%solveEMField() + + END SUBROUTINE doEMField + + !Solving the Poisson equation for electrostatic potential + SUBROUTINE solveElecField() + USE moduleEM + USE moduleMesh + USE moduleErrors + IMPLICIT NONE + + INTEGER, SAVE:: INFO + INTEGER:: n + REAL(8), ALLOCATABLE, SAVE:: tempF(:) + EXTERNAL:: dgetrs + + !$OMP SINGLE + ALLOCATE(tempF(1:mesh%numNodes)) + !$OMP END SINGLE + + CALL assembleSourceVector(tempF) + + !$OMP SINGLE + CALL dgetrs('N', mesh%numNodes, 1, mesh%K, mesh%numNodes, & + mesh%IPIV, tempF, mesh%numNodes, info) + !$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 = tempF(n) + + END DO + !$OMP END DO + + ELSE + !$OMP SINGLE + CALL criticalError('Poisson equation failed', 'solveElecField') + !$OMP END SINGLE + + END IF + + !$OMP SINGLE + DEALLOCATE(tempF) + !$OMP END SINGLE + + END SUBROUTINE solveElecField + + !Empty module that does no computation of EM field for neutral cases + SUBROUTINE noEMField() + IMPLICIT NONE + + END SUBROUTINE noEMField + SUBROUTINE doOutput(t) USE moduleMesh USE moduleOutput @@ -172,8 +308,9 @@ MODULE moduleSolver CALL mesh%printOutput(t) CALL printTime(t, t == 0) CALL mesh%printColl(t) - WRITE(*, "(5X,A21,I10,A1,I8)") "t/tmax: ", t, "/", tmax - WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld + CALL mesh%printEM(t) + WRITE(*, "(5X,A21,I10,A1,I10)") "t/tmax: ", t, "/", tmax + WRITE(*, "(5X,A21,I10)") "Particles: ", nPartOld IF (t == 0) THEN WRITE(*, "(5X,A21,F8.1,A2)") " init time: ", 1.D3*tStep, "ms" @@ -183,7 +320,7 @@ MODULE moduleSolver END IF IF (nPartOld > 0) THEN - WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*(tPush+tReset+tColl+tWeight)/DBLE(nPartOld), "ns" + WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns" END IF WRITE(*,*) diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95 index 9cf9875..dce36a4 100644 --- a/src/modules/moduleSpecies.f95 +++ b/src/modules/moduleSpecies.f95 @@ -31,9 +31,10 @@ MODULE moduleSpecies REAL(8):: v(1:3) !Velocity INTEGER:: pt !Particle species id INTEGER:: e_p !Index of element in which the particle is located - REAL(8):: xLog(1:2) !Logical coordinates of particle in element e_p. + REAL(8):: xLog(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 + REAL(8):: weight=0.D0 !weight of particle + REAL(8):: qm = 0.D0 !charge over mass END TYPE particle