From 37b0139b1f4508ce649a3b675b1f59565a2f9729 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Sun, 13 Dec 2020 13:56:48 +0100 Subject: [PATCH] Implementation of different distribution functions for velocities. Maxwellian and Diract Delta distributions have been implemented. The input for injection of particles should be rewritten to allow more clear input file. --- makefile | 21 +- runs/1D_Cathode/input.json | 3 +- runs/ALPHIE_Grid/input.json | 10 +- runs/ALPHIE_Grid/inputDiffTau.json | 55 + runs/cylFlow/input.json | 3 +- src/fpakc.f95 | 103 -- src/makefile | 4 +- src/modules/makefile | 42 +- src/modules/mesh/1DCart/makefile | 12 +- src/modules/mesh/1DCart/moduleMesh1D.f95 | 489 -------- .../mesh/1DCart/moduleMesh1DBoundary.f95 | 40 - src/modules/mesh/1DCart/moduleMesh1DRead.f95 | 268 ----- src/modules/mesh/Cyl/makefile | 11 - src/modules/mesh/Cyl/moduleMeshCyl.f95 | 1022 ----------------- .../mesh/Cyl/moduleMeshCylBoundary.f95 | 80 -- src/modules/mesh/Cyl/moduleMeshCylRead.f95 | 610 ---------- src/modules/mesh/makefile | 10 +- src/modules/mesh/moduleMesh.f90 | 1 - src/modules/mesh/moduleMesh.f95 | 606 ---------- src/modules/moduleBoundary.f95 | 36 - src/modules/moduleCaseParam.f95 | 9 - src/modules/moduleCollisions.f95 | 171 --- src/modules/moduleCompTime.f95 | 11 - src/modules/moduleConstParam.f95 | 15 - src/modules/moduleEM.f95 | 128 --- src/modules/moduleErrors.f95 | 43 - src/modules/moduleInject.f90 | 110 +- src/modules/moduleInject.f95 | 202 ---- src/modules/moduleInput.f90 | 48 + src/modules/moduleInput.f95 | 559 --------- src/modules/moduleList.f95 | 91 -- src/modules/moduleOutput.f95 | 131 --- src/modules/moduleParallel.f95 | 20 - src/modules/moduleRefParam.f95 | 9 - src/modules/moduleSolver.f95 | 584 ---------- src/modules/moduleSpecies.f95 | 71 -- src/modules/moduleTable.f95 | 121 -- 37 files changed, 252 insertions(+), 5497 deletions(-) create mode 100644 runs/ALPHIE_Grid/inputDiffTau.json delete mode 100644 src/fpakc.f95 delete mode 100644 src/modules/mesh/1DCart/moduleMesh1D.f95 delete mode 100644 src/modules/mesh/1DCart/moduleMesh1DBoundary.f95 delete mode 100644 src/modules/mesh/1DCart/moduleMesh1DRead.f95 delete mode 100644 src/modules/mesh/Cyl/makefile delete mode 100644 src/modules/mesh/Cyl/moduleMeshCyl.f95 delete mode 100644 src/modules/mesh/Cyl/moduleMeshCylBoundary.f95 delete mode 100644 src/modules/mesh/Cyl/moduleMeshCylRead.f95 delete mode 100644 src/modules/mesh/moduleMesh.f95 delete mode 100644 src/modules/moduleBoundary.f95 delete mode 100644 src/modules/moduleCaseParam.f95 delete mode 100644 src/modules/moduleCollisions.f95 delete mode 100644 src/modules/moduleCompTime.f95 delete mode 100644 src/modules/moduleConstParam.f95 delete mode 100644 src/modules/moduleEM.f95 delete mode 100644 src/modules/moduleErrors.f95 delete mode 100644 src/modules/moduleInject.f95 delete mode 100644 src/modules/moduleInput.f95 delete mode 100644 src/modules/moduleList.f95 delete mode 100644 src/modules/moduleOutput.f95 delete mode 100644 src/modules/moduleParallel.f95 delete mode 100644 src/modules/moduleRefParam.f95 delete mode 100644 src/modules/moduleSolver.f95 delete mode 100644 src/modules/moduleSpecies.f95 delete mode 100644 src/modules/moduleTable.f95 diff --git a/makefile b/makefile index d0bf091..c220ac8 100644 --- a/makefile +++ b/makefile @@ -1,21 +1,26 @@ -# compiler -FC := gfortran - +# set folders TOPDIR = $(PWD)# top directory MODDIR := $(TOPDIR)/mod# module folder OBJDIR := $(TOPDIR)/obj# object folder SRCDIR := $(TOPDIR)/src# source folder -JSONDIR := $(TOPDIR)/json-fortran-8.2.0/build +# compiler +# gfortran: +FC := gfortran +JSONDIR := $(TOPDIR)/json-fortran-8.2.0/build-gfortran +# ifort: +# FC := ifort +# JSONDIR := $(TOPDIR)/json-fortran-8.2.0/build-ifort JSONLIB := $(JSONDIR)/lib/libjsonfortran.a JSONINC := $(JSONDIR)/include INCDIR := INCDIR += $(JSONINC) -VAR := "" - -# compile flags -FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -g +# compiler flags +# gfortran: +FCFLAGS := -fopenmp -Ofast -J $(MODDIR) -I $(INCDIR) -Wall -march=native -g +# ifort: +# FCFLAGS := -qopenmp -Ofast -xHost -module $(MODDIR) -I $(INCDIR) -warn all -parallel -free -g #Output file OUTPUT = fpakc diff --git a/runs/1D_Cathode/input.json b/runs/1D_Cathode/input.json index fe9e9d5..330eefe 100644 --- a/runs/1D_Cathode/input.json +++ b/runs/1D_Cathode/input.json @@ -29,7 +29,8 @@ {"name": "Infinite", "type": "dirichlet", "potential": 100.0, "physicalSurface": 2} ], "inject": [ - {"name": "Cathode Electron", "species": "Electron", "flow": 9.0e-5, "units": "A", "v": 27500.0, "T": [2500.0, 2500.0, 2500.0], "n": [ 1, 0, 0], "physicalSurface": 1} + {"name": "Cathode Electron", "species": "Electron", "flow": 9.0e-5, "units": "A", "v": 27500.0, "T": [2500.0, 2500.0, 2500.0], + "velDist": ["Delta", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1} ], "case": { "tau": [1.0e-11], diff --git a/runs/ALPHIE_Grid/input.json b/runs/ALPHIE_Grid/input.json index 1e1e360..010d329 100644 --- a/runs/ALPHIE_Grid/input.json +++ b/runs/ALPHIE_Grid/input.json @@ -25,14 +25,16 @@ {"name": "Axis", "type": "axis", "physicalSurface": 6} ], "boundaryEM": [ - {"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}, {"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4}, {"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5} ], "inject": [ - {"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0], "n": [ 1, 0, 0], "physicalSurface": 1}, - {"name": "Ionization Electron", "species": "Electron", "flow": 27.0e-6, "units": "A", "v": 87000.0, "T": [ 500.0, 500.0, 500.0], "n": [ 1, 0, 0], "physicalSurface": 1}, - {"name": "Cathode Electron", "species": "Electron", "flow": 9.0e-5, "units": "A", "v": 87000.0, "T": [2500.0, 2500.0, 2500.0], "n": [-1, 0, 0], "physicalSurface": 2} + {"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}, + {"name": "Ionization Electron", "species": "Electron", "flow": 27.0e-6, "units": "A", "v": 87000.0, "T": [ 500.0, 500.0, 500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}, + {"name": "Cathode Electron", "species": "Electron", "flow": 9.0e-5, "units": "A", "v": 87000.0, "T": [2500.0, 2500.0, 2500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2} ], "reference": { "density": 1.0e16, diff --git a/runs/ALPHIE_Grid/inputDiffTau.json b/runs/ALPHIE_Grid/inputDiffTau.json new file mode 100644 index 0000000..e8fda33 --- /dev/null +++ b/runs/ALPHIE_Grid/inputDiffTau.json @@ -0,0 +1,55 @@ +{ + "output": { + "path": "./runs/ALPHIE_Grid/", + "triggerOutput": 500, + "triggerCPUTime": 1, + "cpuTime": true, + "numColl": false, + "EMField": true + }, + "geometry": { + "type": "2DCyl", + "meshType": "gmsh", + "meshFile": "mesh.msh" + }, + "species": [ + {"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1}, + {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1} + ], + "boundary": [ + {"name": "Ionization Chanber", "type": "absorption", "physicalSurface": 1}, + {"name": "Vacuum Chamber", "type": "absorption", "physicalSurface": 2}, + {"name": "Exterior", "type": "reflection", "physicalSurface": 3}, + {"name": "Grid Extraction", "type": "absorption", "physicalSurface": 4}, + {"name": "Grid Acceleration", "type": "absorption", "physicalSurface": 5}, + {"name": "Axis", "type": "axis", "physicalSurface": 6} + ], + "boundaryEM": [ + {"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4}, + {"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5} + ], + "inject": [ + {"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}, + {"name": "Ionization Electron", "species": "Electron", "flow": 27.0e-6, "units": "A", "v": 87000.0, "T": [ 500.0, 500.0, 500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1}, + {"name": "Cathode Electron", "species": "Electron", "flow": 9.0e-5, "units": "A", "v": 87000.0, "T": [2500.0, 2500.0, 2500.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2} + ], + "reference": { + "density": 1.0e16, + "mass": 9.109e-31, + "temperature": 2500.0 + }, + "case": { + "tau": [1.0e-9, 1.0e-11], + "time": 2.0e-6, + "pusher": ["2DCylCharged", "2DCylCharged"], + "EMSolver": "Electrostatic" + }, + "parallel": { + "OpenMP":{ + "nThreads": 8 + } + } +} diff --git a/runs/cylFlow/input.json b/runs/cylFlow/input.json index 3dd4b77..c9db450 100644 --- a/runs/cylFlow/input.json +++ b/runs/cylFlow/input.json @@ -21,7 +21,8 @@ {"name": "Axis", "type": "axis", "physicalSurface": 5} ], "inject": [ - {"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} + {"name": "Nozzle", "species": "Argon", "flow": 10.0, "units": "sccm", "v": 300.0, "T": [300.0, 300.0, 300.0], + "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [1, 0, 0], "physicalSurface": 1} ], "reference": { "density": 1.0e20, diff --git a/src/fpakc.f95 b/src/fpakc.f95 deleted file mode 100644 index a24b492..0000000 --- a/src/fpakc.f95 +++ /dev/null @@ -1,103 +0,0 @@ -! FPAKC main program -PROGRAM fpakc - USE moduleInput - USE moduleErrors - USE OMP_LIB - USE moduleInject - USE moduleSolver - USE moduleOutput - USE moduleCompTime - USE moduleMesh - IMPLICIT NONE - - ! 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", "fpakc") - - !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 DEFAULT(SHARED) - DO t = 1, tmax - !Insert new particles and push them - !$OMP SINGLE - tStep = omp_get_wtime() - - !Checks if a species needs to me moved in this iteration - CALL solver%updatePushSpecies(t) - tPush = omp_get_wtime() - !$OMP END SINGLE - - !Injects new particles - CALL doInjects() - - !Push old particles - CALL doPushes() - - !$OMP SINGLE - tPush = omp_get_wtime() - tPush - - !Collisions - tColl = omp_get_wtime() - !$OMP END SINGLE - - CALL doCollisions() - - !$OMP SINGLE - tColl = omp_get_wtime() - tColl - - !Reset particles - tReset = omp_get_wtime() - !$OMP END SINGLE - - CALL doReset() - - !$OMP SINGLE - tReset = omp_get_wtime() - tReset - - !Weight - tWeight = omp_get_wtime() - !$OMP END SINGLE - - CALL doScatter() - - !$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) - !$OMP END SINGLE - - END DO - !$OMP END PARALLEL - -END PROGRAM fpakc - diff --git a/src/makefile b/src/makefile index 6f00e3a..42dc995 100644 --- a/src/makefile +++ b/src/makefile @@ -9,8 +9,8 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolv all: $(OUTPUT) -$(OUTPUT): modules.o $(OUTPUT).f95 - $(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f95 +$(OUTPUT): modules.o $(OUTPUT).f90 + $(FC) $(FCFLAGS) -o $(OBJDIR)/$(OUTPUT).o -c $(OUTPUT).f90 $(FC) $(FCFLAGS) -o $(TOPDIR)/$(OUTPUT) $(OBJECTS) $(OBJDIR)/$(OUTPUT).o $(JSONLIB) -L/usr/local/lib -llapack -lopenblas modules.o: diff --git a/src/modules/makefile b/src/modules/makefile index ea6a1d6..6c8895b 100644 --- a/src/modules/makefile +++ b/src/modules/makefile @@ -9,35 +9,35 @@ 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.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleCollisions.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.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -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.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleInject.o: moduleSpecies.o moduleSolver.o moduleInject.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleRefParam.o: moduleConstParam.o moduleRefParam.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleRefParam.o: moduleConstParam.o moduleRefParam.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleSpecies.o: moduleErrors.o moduleCaseParam.o moduleSpecies.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleSpecies.o: moduleErrors.o moduleCaseParam.o moduleSpecies.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleSolver.o: mesh.o moduleList.o moduleEM.o moduleSpecies.o moduleRefParam.o moduleOutput.o moduleSolver.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleSolver.o: mesh.o moduleList.o moduleEM.o moduleSpecies.o moduleRefParam.o moduleOutput.o moduleSolver.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleTable.o: moduleErrors.o moduleTable.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleTable.o: moduleErrors.o moduleTable.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleEM.o: mesh.o moduleSpecies.o moduleEM.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleEM.o: mesh.o moduleSpecies.o moduleEM.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -%.o: %.f95 +%.o: %.f90 $(FC) $(FCFLAGS) -c $< -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/1DCart/makefile b/src/modules/mesh/1DCart/makefile index 5218e9d..b5b9ff2 100644 --- a/src/modules/mesh/1DCart/makefile +++ b/src/modules/mesh/1DCart/makefile @@ -1,11 +1,11 @@ all: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.o -moduleMesh1D.o: moduleMesh1D.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleMesh1D.o: moduleMesh1D.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMesh1DBoundary.o: moduleMesh1D.o moduleMesh1DBoundary.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleMesh1DBoundary.o: moduleMesh1D.o moduleMesh1DBoundary.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMesh1DRead.o: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleMesh1DRead.o: moduleMesh1D.o moduleMesh1DBoundary.o moduleMesh1DRead.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/1DCart/moduleMesh1D.f95 b/src/modules/mesh/1DCart/moduleMesh1D.f95 deleted file mode 100644 index 2f320ab..0000000 --- a/src/modules/mesh/1DCart/moduleMesh1D.f95 +++ /dev/null @@ -1,489 +0,0 @@ -!moduleMesh1D: 1D cartesian module -! x == x -! y == unused -! z == unused -MODULE moduleMesh1D - USE moduleMesh - IMPLICIT NONE - - TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1D - !Element coordinates - REAL(8):: x = 0.D0 - CONTAINS - PROCEDURE, PASS:: init => initNode1D - PROCEDURE, PASS:: getCoordinates => getCoord1D - - END TYPE meshNode1D - - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdge1D - !Element coordinates - REAL(8):: x = 0.D0 - !Connectivity to nodes - CLASS(meshNode), POINTER:: n1 => NULL() - CONTAINS - PROCEDURE, PASS:: init => initEdge1D - PROCEDURE, PASS:: getNodes => getNodes1D - PROCEDURE, PASS:: randPos => randPos1D - - END TYPE meshEdge1D - - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1D - CONTAINS - PROCEDURE, PASS:: detJac => detJ1D - PROCEDURE, PASS:: invJac => invJ1D - PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshVol1D - - ABSTRACT INTERFACE - PURE FUNCTION fPsi_interface(xi) RESULT(fPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - END FUNCTION fPsi_interface - - PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - END FUNCTION dPsi_interface - - PURE SUBROUTINE partialDer_interface(self, dPsi, dx) - IMPORT meshVol1D - CLASS(meshVol1D), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1):: dx - - END SUBROUTINE partialDer_interface - - END INTERFACE - - TYPE, PUBLIC, EXTENDS(meshVol1D):: meshVol1DSegm - !Element coordinates - REAL(8):: x(1:2) - !Connectivity to nodes - CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() - !Connectivity to adjacent elements - CLASS(*), POINTER:: e1 => NULL(), e2 => NULL() - REAL(8):: arNodes(1:2) - CONTAINS - PROCEDURE, PASS:: init => initVol1DSegm - PROCEDURE, PASS:: area => areaSegm - PROCEDURE, NOPASS:: fPsi => fPsiSegm - PROCEDURE, NOPASS:: dPsi => dPsiSegm - PROCEDURE, PASS:: partialDer => partialDerSegm - PROCEDURE, PASS:: elemK => elemKSegm - PROCEDURE, PASS:: elemF => elemFSegm - PROCEDURE, NOPASS:: weight => weightSegm - PROCEDURE, NOPASS:: inside => insideSegm - PROCEDURE, PASS:: scatter => scatterSegm - PROCEDURE, PASS:: gatherEF => gatherEFSegm - PROCEDURE, PASS:: getNodes => getNodesSegm - PROCEDURE, PASS:: phy2log => phy2logSegm - PROCEDURE, PASS:: nextElement => nextElementSegm - PROCEDURE, PASS:: resetOutput => resetOutputSegm - - END TYPE meshVol1DSegm - - CONTAINS - !NODE FUNCTIONS - !Init node element - SUBROUTINE initNode1D(self, n, r) - USE moduleSpecies - USE moduleRefParam - IMPLICIT NONE - - CLASS(meshNode1D), INTENT(out):: self - INTEGER, INTENT(in):: n - REAL(8), INTENT(in):: r(1:3) - - self%n = n - self%x = r(1)/L_ref - !Node volume, to be determined in mesh - self%v = 0.D0 - - !Allocates output - ALLOCATE(self%output(1:nSpecies)) - - END SUBROUTINE initNode1D - - PURE FUNCTION getCoord1D(self) RESULT(r) - IMPLICIT NONE - - CLASS(meshNode1D), INTENT(in):: self - REAL(8):: r(1:3) - - r = (/ self%x, 0.D0, 0.D0 /) - - END FUNCTION getCoord1D - - !EDGE FUNCTIONS - !Inits edge element - SUBROUTINE initEdge1D(self, n, p, bt, physicalSurface) - IMPLICIT NONE - - CLASS(meshEdge1D), INTENT(out):: self - INTEGER, INTENT(in):: n - INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt - INTEGER, INTENT(in):: physicalSurface - REAL(8), DIMENSION(1:3):: r1 - - self%n = n - self%n1 => mesh%nodes(p(1))%obj - !Get element coordinates - r1 = self%n1%getCoordinates() - - self%x = r1(1) - - self%normal = (/ 1.D0, 0.D0, 0.D0 /) - - !Boundary index - self%bt = bt - !Physical Surface - self%physicalSurface = physicalSurface - - END SUBROUTINE initEdge1D - - !Get nodes from edge - PURE FUNCTION getNodes1D(self) RESULT(n) - IMPLICIT NONE - - CLASS(meshEdge1D), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - ALLOCATE(n(1)) - n = (/ self%n1%n /) - - END FUNCTION getNodes1D - - !Calculates a 'random' position in edge - FUNCTION randPos1D(self) RESULT(r) - CLASS(meshEdge1D), INTENT(in):: self - REAL(8):: r(1:3) - - r = (/ self%x, 0.D0, 0.D0 /) - - END FUNCTION randPos1D - - !VOLUME FUNCTIONS - !SEGMENT FUNCTIONS - !Init segment element - SUBROUTINE initVol1DSegm(self, n, p) - USE moduleRefParam - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(out):: self - INTEGER, INTENT(in):: n - INTEGER, INTENT(in):: p(:) - REAL(8), DIMENSION(1:3):: r1, r2 - - self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - !Get element coordinates - r1 = self%n1%getCoordinates() - r2 = self%n2%getCoordinates() - self%x = (/ r1(1), r2(1) /) - - !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - - self%sigmaVrelMax = sigma_ref/L_ref**2 - - CALL OMP_INIT_LOCK(self%lock) - - END SUBROUTINE initVol1DSegm - - !Computes element area - PURE SUBROUTINE areaSegm(self) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(inout):: self - REAL(8):: l !element length - REAL(8):: fPsi(1:2) - REAL(8):: detJ - REAL(8):: Xii(1:3) - - self%volume = 0.D0 - self%arNodes = 0.D0 - !1 point Gauss integral - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) - l = 2.D0*detJ - self%volume = l - self%arNodes = fPsi*l - - END SUBROUTINE areaSegm - - !Computes element functions at point xii - PURE FUNCTION fPsiSegm(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:2)) - - fPsi(1) = 1.D0 - xi(1) - fPsi(2) = 1.D0 + xi(1) - fPsi = fPsi * 5.D-1 - - END FUNCTION fPsiSegm - - !Computes element derivative shape function at Xii - PURE FUNCTION dPsiSegm(xi) RESULT(dPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - ALLOCATE(dPsi(1:1, 1:2)) - - dPsi(1, 1) = -5.D-1 - dPsi(1, 2) = 5.D-1 - - END FUNCTION dPsiSegm - - !Computes partial derivatives of coordinates - PURE SUBROUTINE partialDerSegm(self, dPsi, dx) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1):: dx - - dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) - - END SUBROUTINE partialDerSegm - - !Computes local stiffness matrix - PURE FUNCTION elemKSegm(self) RESULT(ke) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - REAL(8):: ke(1:2,1:2) - REAL(8):: Xii(1:3) - REAL(8):: dPsi(1:1, 1:2) - REAL(8):: invJ - - ke = 0.D0 - Xii = (/ 0.D0, 0.D0, 0.D0 /) - dPsi = self%dPsi(Xii) - invJ = self%invJac(Xii, dPsi) - ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) - ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) - ke = 2.D0*ke*invJ - - END FUNCTION elemKSegm - - PURE FUNCTION elemFSegm(self, source) RESULT(localF) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: fPsi(1:2) - REAL(8):: detJ - REAL(8):: Xii(1:3) - - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) - ALLOCATE(localF(1:2)) - localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ - - END FUNCTION elemFSegm - - PURE FUNCTION weightSegm(xi) RESULT(w) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: w(1:3) - - w = fPsiSegm(xi) - - END FUNCTION weightSegm - - PURE FUNCTION insideSegm(xi) RESULT(ins) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - LOGICAL:: ins - - ins = xi(1) >=-1.D0 .AND. & - xi(1) <= 1.D0 - - END FUNCTION insideSegm - - SUBROUTINE scatterSegm(self, part) - USE moduleOutput - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - CLASS(particle), INTENT(in):: part - TYPE(outputNode), POINTER:: vertex - REAL(8):: w_p(1:2) - REAL(8):: tensorS(1:3,1:3) - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - 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%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 - - END SUBROUTINE scatterSegm - - !Gathers EF at position Xii - PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1, 1:2) - REAL(8):: phi(1:2) - REAL(8):: EF(1:3) - REAL(8):: invJ - - phi = (/ self%n1%emData%phi, & - self%n2%emData%phi /) - - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - EF(1) = -DOT_PRODUCT(dPsi(1, :), phi)*invJ - EF(2) = 0.D0 - EF(3) = 0.D0 - - END FUNCTION gatherEFSegm - - !Get nodes from 1D volume - PURE FUNCTION getNodesSegm(self) RESULT(n) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - ALLOCATE(n(1:2)) - n = (/ self%n1%n, self%n2%n /) - - END FUNCTION getNodesSegm - - PURE FUNCTION phy2logSegm(self, r) RESULT(xN) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - - xN = 0.D0 - xN(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0 - - END FUNCTION phy2logSegm - - !Get next element for a logical position xi - SUBROUTINE nextElementSegm(self, xi, nextElement) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(*), POINTER, INTENT(out):: nextElement - - NULLIFY(nextElement) - IF (xi(1) < -1.D0) THEN - nextElement => self%e2 - - ELSEIF (xi(1) > 1.D0) THEN - nextElement => self%e1 - - END IF - - END SUBROUTINE nextElementSegm - - !Reset the output of nodes in element - PURE SUBROUTINE resetOutputSegm(self) - USE moduleSpecies - USE moduleOutput - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(inout):: self - INTEGER:: k - - DO k = 1, nSpecies - self%n1%output(k)%den = 0.D0 - self%n1%output(k)%mom = 0.D0 - self%n1%output(k)%tensorS = 0.D0 - - self%n2%output(k)%den = 0.D0 - self%n2%output(k)%mom = 0.D0 - self%n2%output(k)%tensorS = 0.D0 - - END DO - - END SUBROUTINE resetOutputSegm - - - !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS - !Computes the element Jacobian determinant - PURE FUNCTION detJ1D(self, xi, dPsi_in) RESULT(dJ) - IMPLICIT NONE - - CLASS(meshVol1D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) - REAL(8):: dJ - REAL(8):: dx(1) - - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx) - dJ = dx(1) - - END FUNCTION detJ1D - - !Computes the invers Jacobian - PURE FUNCTION invJ1D(self, xi, dPsi_in) RESULT(invJ) - IMPLICIT NONE - - CLASS(meshVol1D), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) - REAL(8):: dx(1) - REAL(8):: invJ - - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dx) - invJ = 1.D0/dx(1) - - END FUNCTION invJ1D - - -END MODULE moduleMesh1D - diff --git a/src/modules/mesh/1DCart/moduleMesh1DBoundary.f95 b/src/modules/mesh/1DCart/moduleMesh1DBoundary.f95 deleted file mode 100644 index bd4e6d6..0000000 --- a/src/modules/mesh/1DCart/moduleMesh1DBoundary.f95 +++ /dev/null @@ -1,40 +0,0 @@ -MODULE moduleMesh1DBoundary - USE moduleMesh1D - - TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DRef - CONTAINS - PROCEDURE, PASS:: fBoundary => reflection - - END TYPE meshEdge1DRef - - TYPE, PUBLIC, EXTENDS(meshEdge1D):: meshEdge1DAbs - CONTAINS - PROCEDURE, PASS:: fBoundary => absorption - - END TYPE meshEdge1DAbs - - CONTAINS - SUBROUTINE reflection(self, part) - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshEdge1DRef), INTENT(inout):: self - CLASS(particle), INTENT(inout):: part - - part%v(1) = -part%v(1) - part%r(1) = 2.D0*self%x - part%r(1) - - END SUBROUTINE reflection - - SUBROUTINE absorption(self, part) - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshEdge1DAbs), INTENT(inout):: self - CLASS(particle), INTENT(inout):: part - - part%n_in = .FALSE. - - END SUBROUTINE absorption - -END MODULE moduleMesh1DBoundary diff --git a/src/modules/mesh/1DCart/moduleMesh1DRead.f95 b/src/modules/mesh/1DCart/moduleMesh1DRead.f95 deleted file mode 100644 index 1c52795..0000000 --- a/src/modules/mesh/1DCart/moduleMesh1DRead.f95 +++ /dev/null @@ -1,268 +0,0 @@ -MODULE moduleMesh1DRead - USE moduleMesh - USE moduleMesh1D - USE moduleMesh1DBoundary - - !TODO: make this abstract to allow different mesh formats - TYPE, EXTENDS(meshGeneric):: mesh1DGeneric - CONTAINS - PROCEDURE, PASS:: init => init1DMesh - PROCEDURE, PASS:: readMesh => readMesh1D - - END TYPE - - INTERFACE connected - MODULE PROCEDURE connectedVolVol, connectedVolEdge - - END INTERFACE connected - - CONTAINS - !Init 1D mesh - SUBROUTINE init1DMesh(self, meshFormat) - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - CLASS(mesh1DGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - SELECT CASE(meshFormat) - CASE ("gmsh") - self%printOutput => printOutputGmsh - self%printColl => printCollGmsh - self%printEM => printEMGmsh - - CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "init1D") - - END SELECT - - END SUBROUTINE init1DMesh - - !Reads 1D mesh - SUBROUTINE readMesh1D(self, filename) - USE moduleBoundary - IMPLICIT NONE - - CLASS(mesh1DGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: x - INTEGER:: p(1:2) - INTEGER:: e, et, n, eTemp, elemType, bt - INTEGER:: totalNumElem - INTEGER:: boundaryType - - !Open file mesh - OPEN(10, FILE=TRIM(filename)) - !Skip header - READ(10, *) - READ(10, *) - READ(10, *) - READ(10, *) - !Read number of nodes - READ(10, *) self%numNodes - !Allocate required matrices and vectors - ALLOCATE(self%nodes(1:self%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 coordinates. Only relevant for x - DO e = 1, self%numNodes - READ(10, *) n, x - ALLOCATE(meshNode1D:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/ x, 0.D0, 0.D0 /)) - - END DO - !Skips comments - READ(10, *) - READ(10, *) - !Reads the total number of elements (edges+vol) - READ(10, *) totalNumElem - self%numEdges = 0 - DO e = 1, totalNumElem - READ(10, *) eTemp, elemType - IF (elemType == 15) THEN !15 is physical node in GMSH - self%numEdges = e - - END IF - - END DO - - !Substract the number of edges to the total number of elements - !to obtain the number of volume elements - self%numVols = totalNumelem - self%numEdges - !Allocates arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beginning of reading elements - DO e = 1, totalNumelem - BACKSPACE(10) - - END DO - - !Reads edges - DO e = 1, self%numEdges - READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1) - !Associate boundary condition - bt = getBoundaryId(boundaryType) - SELECT CASE(boundary(bt)%obj%boundaryType) - CASE ('reflection') - ALLOCATE(meshEdge1DRef:: self%edges(e)%obj) - - CASE ('absorption') - ALLOCATE(meshEdge1DAbs:: self%edges(e)%obj) - - END SELECT - - CALL self%edges(e)%obj%init(n, p(1:1), bt, boundaryType) - - END DO - - !Read and initialize volumes - DO e = 1, self%numVols - READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:2) - ALLOCATE(meshVol1DSegm:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:2)) - - END DO - - CLOSE(10) - - !Build connectivity between elements - DO e = 1, self%numVols - !Connectivity between volumes - DO et = 1, self%numVols - IF (e /= et) THEN - CALL connected(self%vols(e)%obj, self%vols(et)%obj) - - END IF - - END DO - - !Connectivity betwen vols and edges - DO et = 1, self%numEdges - CALL connected(self%vols(e)%obj, self%edges(et)%obj) - - END DO - - !Constructs the global K matrix - CALL constructGlobalK(self%K, self%vols(e)%obj) - - END DO - - END SUBROUTINE readMesh1D - - SUBROUTINE connectedVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVol1DSegm) - SELECT TYPE(elemB) - TYPE IS(meshVol1DSegm) - CALL connectedSegmSegm(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolVol - - SUBROUTINE connectedSegmSegm(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA - CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemB - - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - END IF - - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n1%n == elemB%n2%n) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END SUBROUTINE connectedSegmSegm - - SUBROUTINE connectedVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS (meshVol1DSegm) - SELECT TYPE(elemB) - CLASS IS(meshEdge1D) - CALL connectedSegmEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolEdge - - SUBROUTINE connectedSegmEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol1DSegm), INTENT(inout), TARGET:: elemA - CLASS(meshEdge1D), INTENT(inout), TARGET:: elemB - - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n2%n == elemB%n1%n) THEN - - elemA%e1 => elemB - elemB%e2 => elemA - - END IF - - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n1%n == elemB%n1%n) THEN - - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END SUBROUTINE connectedSegmEdge - - SUBROUTINE constructGlobalK(K, elem) - IMPLICIT NONE - - REAL(8), INTENT(inout):: K(1:,1:) - CLASS(meshVol), INTENT(in):: elem - REAL(8):: localK(1:2,1:2) - INTEGER:: i, j - INTEGER:: n(1:2) - - SELECT TYPE(elem) - TYPE IS(meshVol1DSegm) - localK = elem%elemK() - n = (/ elem%n1%n, elem%n2%n /) - - CLASS DEFAULT - n = 0 - localK = 0.D0 - - END SELECT - - DO i = 1, 2 - DO j = 1, 2 - K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j) - END DO - END DO - - END SUBROUTINE constructGlobalK - -END MODULE moduleMesh1DRead diff --git a/src/modules/mesh/Cyl/makefile b/src/modules/mesh/Cyl/makefile deleted file mode 100644 index 7a60575..0000000 --- a/src/modules/mesh/Cyl/makefile +++ /dev/null @@ -1,11 +0,0 @@ -all : moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.o - -moduleMeshCyl.o: moduleMeshCyl.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ - -moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ - -moduleMeshCylRead.o: moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ - diff --git a/src/modules/mesh/Cyl/moduleMeshCyl.f95 b/src/modules/mesh/Cyl/moduleMeshCyl.f95 deleted file mode 100644 index 2b1362c..0000000 --- a/src/modules/mesh/Cyl/moduleMeshCyl.f95 +++ /dev/null @@ -1,1022 +0,0 @@ -!moduleMeshCyl: 2D axial symmetric extension of generic mesh from GMSH format. -! x == z -! y == r -! z == theta (unused) -MODULE moduleMeshCyl - USE moduleMesh - IMPLICIT NONE - - !Values for Gauss integral - REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) - REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) - - REAL(8), PARAMETER:: xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /) - REAL(8), PARAMETER:: xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /) - REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /) - - TYPE, PUBLIC, EXTENDS(meshNode):: meshNodeCyl - !Element coordinates - REAL(8):: r = 0.D0, z = 0.D0 - CONTAINS - PROCEDURE, PASS:: init => initNodeCyl - PROCEDURE, PASS:: getCoordinates => getCoordCyl - - END TYPE meshNodeCyl - - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshEdge):: meshEdgeCyl - !Element coordinates - REAL(8):: r(1:2) = 0.D0, z(1:2) = 0.D0 - !Connectivity to nodes - CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() - CONTAINS - PROCEDURE, PASS:: init => initEdgeCyl - PROCEDURE, PASS:: getNodes => getNodesCyl - PROCEDURE, PASS:: randPos => randPosCyl - - END TYPE meshEdgeCyl - - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl - CONTAINS - PROCEDURE, PASS:: detJac => detJCyl - PROCEDURE, PASS:: invJac => invJCyl - PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi - PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshVolCyl - - ABSTRACT INTERFACE - PURE FUNCTION fPsi_interface(xi) RESULT(fPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - END FUNCTION fPsi_interface - - PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - END FUNCTION dPsi_interface - - PURE SUBROUTINE partialDer_interface(self, dPsi, dz, dr) - IMPORT meshVolCyl - CLASS(meshVolCyl), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr - - END SUBROUTINE partialDer_interface - - END INTERFACE - - TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylQuad - !Element coordinates - REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 - !Connectivity to nodes - 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() - REAL(8):: arNodes(1:4) = 0.D0 - - CONTAINS - PROCEDURE, PASS:: init => initVolQuadCyl - PROCEDURE, PASS:: area => areaQuad - PROCEDURE, NOPASS:: fPsi => fPsiQuad - PROCEDURE, NOPASS:: dPsi => dPsiQuad - PROCEDURE, NOPASS:: dPsiXi1 => dPsiQuadXi1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiQuadXi2 - PROCEDURE, PASS:: partialDer => partialDerQuad - PROCEDURE, PASS:: elemK => elemKQuad - PROCEDURE, PASS:: elemF => elemFQuad - PROCEDURE, NOPASS:: weight => weightQuad - PROCEDURE, NOPASS:: inside => insideQuad - PROCEDURE, PASS:: scatter => scatterQuad - PROCEDURE, PASS:: gatherEF => gatherEFQuad - PROCEDURE, PASS:: getNodes => getNodesQuad - PROCEDURE, PASS:: phy2log => phy2logQuad - PROCEDURE, PASS:: nextElement => nextElementQuad - PROCEDURE, PASS:: resetOutput => resetOutputQuad - - END TYPE meshVolCylQuad - - TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylTria - !Element coordinates - REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0 - !Connectivity to nodes - CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() - !Connectivity to adjacent elements - CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL() - REAL(8):: arNodes(1:3) = 0.D0 - !Derivatives in z,r real space - REAL(8), DIMENSION(1:3):: dPsiZ, dPsiR - - CONTAINS - PROCEDURE, PASS:: init => initVolTriaCyl - PROCEDURE, PASS:: area => areaTria - PROCEDURE, NOPASS:: fPsi => fPsiTria - PROCEDURE, NOPASS:: dPsi => dPsiTria - PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2 - PROCEDURE, PASS:: partialDer => partialDerTria - PROCEDURE, PASS:: elemK => elemKTria - PROCEDURE, PASS:: elemF => elemFTria - PROCEDURE, NOPASS:: weight => weightTria - PROCEDURE, NOPASS:: inside => insideTria - PROCEDURE, PASS:: scatter => scatterTria - PROCEDURE, PASS:: gatherEF => gatherEFTria - PROCEDURE, PASS:: getNodes => getNodesTria - PROCEDURE, PASS:: phy2log => phy2logTria - PROCEDURE, PASS:: nextElement => nextElementTria - PROCEDURE, PASS:: resetOutput => resetOutputTria - - END TYPE meshVolCylTria - - CONTAINS - - !NODE FUNCTIONS - !Inits node element - SUBROUTINE initNodeCyl(self, n, r) - USE moduleSpecies - USE moduleRefParam - IMPLICIT NONE - - CLASS(meshNodeCyl), INTENT(out):: self - INTEGER, INTENT(in):: n - REAL(8), INTENT(in):: r(1:3) - - self%n = n - self%z = r(2)/L_ref - self%r = r(1)/L_ref - !Node volume, to be determined in mesh - self%v = 0.D0 - - !Allocates output: - ALLOCATE(self%output(1:nSpecies)) - - END SUBROUTINE initNodeCyl - - PURE FUNCTION getCoordCyl(self) RESULT(r) - IMPLICIT NONE - - CLASS(meshNodeCyl), INTENT(in):: self - REAL(8):: r(1:3) - - r = (/self%z, self%r, 0.D0/) - - END FUNCTION getCoordCyl - - !EDGE FUNCTIONS - !Inits edge element - SUBROUTINE initEdgeCyl(self, n, p, bt, physicalSurface) - IMPLICIT NONE - - CLASS(meshEdgeCyl), INTENT(out):: self - INTEGER, INTENT(in):: n - INTEGER, INTENT(in):: p(:) - INTEGER, INTENT(in):: bt - INTEGER, INTENT(in):: physicalSurface - REAL(8), DIMENSION(1:3):: r1, r2 - - self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - !Get element coordinates - r1 = self%n1%getCoordinates() - r2 = self%n2%getCoordinates() - self%z = (/r1(1), r2(1)/) - self%r = (/r1(2), r2(2)/) - !Normal vector - self%normal = (/ self%r(2)-self%r(1), & - self%z(2)-self%z(1), & - 0.D0 /) - !Boundary index - self%bt = bt - !Phyiscal Surface - self%physicalSurface = physicalSurface - - END SUBROUTINE initEdgeCyl - - !Get nodes from edge - PURE FUNCTION getNodesCyl(self) RESULT(n) - IMPLICIT NONE - - CLASS(meshEdgeCyl), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - ALLOCATE(n(1:2)) - n = (/self%n1%n, self%n2%n /) - - END FUNCTION getNodesCyl - - !Calculates a random position in edge - 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 - r(3) = 0.D0 - - END FUNCTION randPosCyl - - !VOLUME FUNCTIONS - !QUAD FUNCTIONS - !Inits quadrilateral element - SUBROUTINE initVolQuadCyl(self, n, p) - USE moduleRefParam - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(out):: self - INTEGER, INTENT(in):: n - INTEGER, INTENT(in):: p(:) - REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 - - self%n = n - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - self%n3 => mesh%nodes(p(3))%obj - self%n4 => mesh%nodes(p(4))%obj - !Get element coordinates - r1 = self%n1%getCoordinates() - r2 = self%n2%getCoordinates() - r3 = self%n3%getCoordinates() - r4 = self%n4%getCoordinates() - self%z = (/r1(1), r2(1), r3(1), r4(1)/) - self%r = (/r1(2), r2(2), r3(2), r4(2)/) - - !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) - self%n4%v = self%n4%v + self%arNodes(4) - - self%sigmaVrelMax = sigma_ref/L_ref**2 - - CALL OMP_INIT_LOCK(self%lock) - - END SUBROUTINE initVolQuadCyl - - !Computes element area - PURE SUBROUTINE areaQuad(self) - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(inout):: self - REAL(8):: r, xi(1:3) - REAL(8):: detJ - REAL(8):: fPsi(1:4) - - self%volume = 0.D0 - self%arNodes = 0.D0 - !2D 1 point Gauss Quad Integral - xi = 0.D0 - detJ = self%detJac(xi)*8.D0*PI !4*2*pi - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - self%volume = r*detJ - self%arNodes = fPsi*r*detJ - - END SUBROUTINE areaQuad - - !Computes element functions in point xi - PURE FUNCTION fPsiQuad(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8),INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:4)) - - fPsi(1) = (1.D0-xi(1))*(1.D0-xi(2)) - fPsi(2) = (1.D0+xi(1))*(1.D0-xi(2)) - fPsi(3) = (1.D0+xi(1))*(1.D0+xi(2)) - fPsi(4) = (1.D0-xi(1))*(1.D0+xi(2)) - fPsi = fPsi*0.25D0 - - END FUNCTION fPsiQuad - - !Derivative element function at coordinates xi - PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - ALLOCATE(dPsi(1:2,1:4)) - - dPsi(1,:) = dPsiQuadXi1(xi(2)) - dPsi(2,:) = dPsiQuadXi2(xi(1)) - - END FUNCTION dPsiQuad - - !Derivative element function (xi1) - PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1) - IMPLICIT NONE - - REAL(8),INTENT(in):: xi2 - REAL(8):: dPsiXi1(1:4) - - dPsiXi1(1) = -(1.D0-xi2) - dPsiXi1(2) = (1.D0-xi2) - dPsiXi1(3) = (1.D0+xi2) - dPsiXi1(4) = -(1.D0+xi2) - dPsiXi1 = dPsiXi1*0.25D0 - - END FUNCTION dPsiQuadXi1 - - !Derivative element function (xi2) - PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2) - IMPLICIT NONE - - REAL(8),INTENT(in):: xi1 - REAL(8):: dPsiXi2(1:4) - - dPsiXi2(1) = -(1.D0-xi1) - dPsiXi2(2) = -(1.D0+xi1) - dPsiXi2(3) = (1.D0+xi1) - dPsiXi2(4) = (1.D0-xi1) - dPsiXi2 = dPsiXi2*0.25D0 - - END FUNCTION dPsiQuadXi2 - - PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr - - 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) - - END SUBROUTINE partialDerQuad - - !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self) RESULT(ke) - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8):: r, xi(1:3) - REAL(8):: fPsi(1:4), dPsi(1:2,1:4) - REAL(8):: ke(1:4,1:4) - REAL(8):: invJ(1:2,1:2), detJ - INTEGER:: l, m - - ke=0.D0 - xi=0.D0 - !Start 2D Gauss Quad Integral - DO l=1, 3 - xi(2) = corQuad(l) - dPsi(1,:) = self%dPsiXi1(xi(2)) - DO m = 1, 3 - xi(1) = corQuad(m) - dPsi(2,:) = self%dPsiXi2(xi(1)) - fPsi = self%fPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - r = DOT_PRODUCT(fPsi,self%r) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wQuad(l)*wQuad(m)/detJ - - END DO - END DO - ke = ke*2.D0*PI - - END FUNCTION elemKQuad - - !Computes the local source vector for a force f - PURE FUNCTION elemFQuad(self, source) RESULT(localF) - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: r, xi(1:3) - REAL(8):: fPsi(1:4) - REAL(8):: detJ, f - INTEGER:: l, m - - ALLOCATE(localF(1:4)) - localF = 0.D0 - xi = 0.D0 - DO l=1, 3 - xi(1) = corQuad(l) - DO m = 1, 3 - xi(2) = corQuad(m) - detJ = self%detJac(xi) - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - f = DOT_PRODUCT(fPsi,source) - localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ - - END DO - END DO - localF = localF*2.D0*PI - - END FUNCTION elemFQuad - - !Computes weights in the element nodes - PURE FUNCTION weightQuad(xi) RESULT(w) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: w(1:4) - - w = fPsiQuad(xi) - - END FUNCTION weightQuad - - !Checks if a particle is inside a quad element - PURE FUNCTION insideQuad(xi) RESULT(ins) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - LOGICAL:: ins - - ins = (xi(1) >= -1.D0 .AND. xi(1) <= 1.D0) .AND. & - (xi(2) >= -1.D0 .AND. xi(2) <= 1.D0) - - END FUNCTION insideQuad - - !Scatter properties of particle into element nodes - SUBROUTINE scatterQuad(self, part) - USE moduleOutput - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - CLASS(particle), INTENT(in):: part - TYPE(outputNode), POINTER:: vertex - REAL(8):: w_p(1:4) - REAL(8):: tensorS(1:3,1:3) - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - 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%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%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%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 - - END SUBROUTINE scatterQuad - - !Gathers the electric field at position xi - PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1:2,1:4) - REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates - REAL(8):: invJ(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 = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - dPsiR = MATMUL(invJ, dPsi)/detJ - EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) - EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) - EF(3) = 0.D0 - - END FUNCTION gatherEFQuad - - !Gets nodes from quadrilateral element - 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 - - !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logQuad(self,r) RESULT(xN) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2) - REAL(8):: dPsi(1:2,1:4), fPsi(1:4) - REAL(8):: conv - - !Iterative newton method to transform coordinates - conv=1.D0 - xO=0.D0 - - DO WHILE(conv>1.D-4) - dPsi = self%dPsi(xO) - invJ = self%invJac(xO, dPsi) - fPsi = self%fPsi(xO) - f(1) = DOT_PRODUCT(fPsi,self%z)-r(1) - f(2) = DOT_PRODUCT(fPsi,self%r)-r(2) - detJ = self%detJac(xO,dPsi) - xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ - conv=MAXVAL(DABS(xN-xO),1) - xO=xN - - END DO - - END FUNCTION phy2logQuad - - !Gets the next element for a logical position xi - SUBROUTINE nextElementQuad(self, xi, nextElement) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(*), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:4) - INTEGER:: nextInt - - xiArray = (/ -xi(2), xi(1), xi(2), -xi(1) /) - nextInt = MAXLOC(xiArray,1) - !Selects the higher value of directions and searches in that direction - NULLIFY(nextElement) - SELECT CASE (nextInt) - CASE (1) - nextElement => self%e1 - CASE (2) - nextElement => self%e2 - CASE (3) - nextElement => self%e3 - CASE (4) - nextElement => self%e4 - END SELECT - - END SUBROUTINE nextElementQuad - - !Reset the output of nodes in quad element - PURE SUBROUTINE resetOutputQuad(self) - USE moduleSpecies - USE moduleOutput - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(inout):: self - INTEGER:: k - - DO k = 1, nSpecies - self%n1%output(k)%den = 0.D0 - self%n1%output(k)%mom = 0.D0 - self%n1%output(k)%tensorS = 0.D0 - - self%n2%output(k)%den = 0.D0 - self%n2%output(k)%mom = 0.D0 - self%n2%output(k)%tensorS = 0.D0 - - self%n3%output(k)%den = 0.D0 - self%n3%output(k)%mom = 0.D0 - self%n3%output(k)%tensorS = 0.D0 - - self%n4%output(k)%den = 0.D0 - self%n4%output(k)%mom = 0.D0 - self%n4%output(k)%tensorS = 0.D0 - - END DO - - END SUBROUTINE resetOutputQuad - - !TRIA ELEMENT - !Init tria element - SUBROUTINE initVolTriaCyl(self, n, p) - USE moduleRefParam - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(out):: self - INTEGER, INTENT(in):: n - INTEGER, INTENT(in):: p(:) - REAL(8), DIMENSION(1:3):: r1, r2, r3 - REAL(8):: A - - !Assign node index - self%n = n - - !Assign nodes to element - self%n1 => mesh%nodes(p(1))%obj - self%n2 => mesh%nodes(p(2))%obj - self%n3 => mesh%nodes(p(3))%obj - !Get element coordinates - r1 = self%n1%getCoordinates() - r2 = self%n2%getCoordinates() - r3 = self%n3%getCoordinates() - self%z = (/r1(1), r2(1), r3(1)/) - self%r = (/r1(2), r2(2), r3(2)/) - !Assign node volume - CALL self%area() - self%n1%v = self%n1%v + self%arNodes(1) - self%n2%v = self%n2%v + self%arNodes(2) - self%n3%v = self%n3%v + self%arNodes(3) - - !Derivatives in z/r for shape functions (node independent) - !TODO: This is used because invJ.dPsi does not produce the right Electric field - A = self%z(2)*self%r(3) - self%z(3)*self%r(2) + & - self%z(3)*self%r(1) - self%z(1)*self%r(3) + & - self%z(1)*self%r(2) - self%z(2)*self%r(1) - - self%dPsiZ = (/ self%r(2)-self%r(3), & - self%r(3)-self%r(1), & - self%r(1)-self%r(2) /)/A - - self%dPsiR = (/ self%z(3)-self%z(2), & - self%z(1)-self%z(3), & - self%z(2)-self%z(1) /)/A - - self%sigmaVrelMax = sigma_ref/L_ref**2 - - CALL OMP_INIT_LOCK(self%lock) - - END SUBROUTINE initVolTriaCyl - - !Calculates area for triangular element - PURE SUBROUTINE areaTria(self) - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(inout):: self - REAL(8):: r, xi(1:3) - REAL(8):: detJ - REAL(8):: fPsi(1:3) - - self%volume = 0.D0 - self%arNodes = 0.D0 - !2D 1 point Gauss Quad Integral - xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) - detJ = self%detJac(xi)*PI !2PI*1/2 - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - self%volume = r*detJ - self%arNodes = fPsi*r*detJ - - END SUBROUTINE areaTria - - !Shape functions for triangular element - PURE FUNCTION fPsiTria(xi) RESULT(fPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:3)) - - fPsi(1) = 1.D0 - xi(1) - xi(2) - fPsi(2) = xi(1) - fPsi(3) = xi(2) - - END FUNCTION fPsiTria - - !Derivative element function at coordinates xi - PURE FUNCTION dPsiTria(xi) RESULT(dPsi) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: dPsi(:,:) - - ALLOCATE(dPsi(1:2,1:3)) - - dPsi(1,:) = dPsiTriaXi1(xi(2)) - dPsi(2,:) = dPsiTriaXi2(xi(1)) - - END FUNCTION dPsiTria - - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi2 - REAL(8):: dPsiXi1(1:3) - - dPsiXi1(1) = -1.D0 - dPsiXi1(2) = 1.D0 - dPsiXi1(3) = 0.D0 - - END FUNCTION dPsiTriaXi1 - - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi1 - REAL(8):: dPsiXi2(1:3) - - dPsiXi2(1) = -1.D0 - dPsiXi2(2) = 0.D0 - dPsiXi2(3) = 1.D0 - - END FUNCTION dPsiTriaXi2 - - PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr) - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:,1:) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr - - 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) - - END SUBROUTINE partialDerTria - - !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self) RESULT(ke) - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - REAL(8):: r, xi(1:3) - REAL(8):: fPsi(1:3), dPsi(1:2,1:3) - REAL(8):: ke(1:3,1:3) - REAL(8):: invJ(1:2,1:2), detJ - INTEGER:: l - - ke=0.D0 - xi=0.D0 - !Start 2D Gauss Quad Integral - DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ - - END DO - ke = ke*2.D0*PI - - END FUNCTION elemKTria - - !Computes element local source vector - PURE FUNCTION elemFTria(self, source) RESULT(localF) - USE moduleConstParam - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - REAL(8), INTENT(in):: source(1:) - REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: fPsi(1:3) - REAL(8):: r, xi(1:3) - REAL(8):: detJ, f - INTEGER:: l - - ALLOCATE(localF(1:3)) - localF = 0.D0 - xi = 0.D0 - !Start 2D Gauss Quad Integral - DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - detJ = self%detJac(xi) - fPsi = self%fPsi(xi) - r = DOT_PRODUCT(fPsi,self%r) - f = DOT_PRODUCT(fPsi,source) - localF = localF + r*f*fPsi*wTria(l)*detJ - - END DO - localF = localF*2.D0*PI - - END FUNCTION elemFTria - - !Computes weights in the element nodes - PURE FUNCTION weightTria(xi) RESULT(w) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: w(:) - - w = fPsiTria(xi) - - END FUNCTION weightTria - - PURE FUNCTION insideTria(xi) RESULT(ins) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - LOGICAL:: ins - - ins = xi(1) >= 0.D0 .AND. & - xi(2) >= 0.D0 .AND. & - 1.D0 - xi(1) - xi(2) >= 0.D0 - - END FUNCTION insideTria - - !Scatter properties of particles into element - SUBROUTINE scatterTria(self, part) - USE moduleOutput - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - CLASS(particle), INTENT(in):: part - TYPE(outputNode), POINTER:: vertex - REAL(8):: w_p(1:3) - REAL(8):: tensorS(1:3,1:3) - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - 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%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%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 - - END SUBROUTINE scatterTria - - !Gathers the electric field at position xi - PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: phi(1:3) - REAL(8):: EF(1:3) - - phi = (/self%n1%emData%phi, & - self%n2%emData%phi, & - self%n3%emData%phi/) - - EF(1) = -DOT_PRODUCT(self%dPsiZ(:), phi) - EF(2) = -DOT_PRODUCT(self%dPsiR(:), phi) - EF(3) = 0.D0 - - END FUNCTION gatherEFTria - - !Gets node indexes from triangular element - PURE FUNCTION getNodesTria(self) RESULT(n) - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) - - ALLOCATE(n(1:3)) - n = (/self%n1%n, self%n2%n, self%n3%n /) - - END FUNCTION getNodesTria - - !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logTria(self,r) RESULT(xi) - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - REAL(8), INTENT(in):: r(1:3) - REAL(8):: xi(1:3) - REAL(8):: invJ(1:2,1:2), detJ - REAL(8):: deltaR(1:2) - REAL(8):: dPsi(1:2,1:3) - - !Direct method to convert coordinates - xi = 0.D0 !Irrelevant, required for input - deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /) - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - detJ = self%detJac(xi, dPsi) - xi(1:2) = MATMUL(invJ,deltaR)/detJ - - END FUNCTION phy2logTria - - SUBROUTINE nextElementTria(self, xi, nextElement) - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(*), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:3) - INTEGER:: nextInt - - xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /) - nextInt = MINLOC(xiArray,1) - NULLIFY(nextElement) - SELECT CASE (nextInt) - CASE (1) - nextElement => self%e1 - CASE (2) - nextElement => self%e2 - CASE (3) - nextElement => self%e3 - END SELECT - - END SUBROUTINE nextElementTria - - !Reset the output of nodes in tria element - PURE SUBROUTINE resetOutputTria(self) - USE moduleSpecies - USE moduleOutput - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(inout):: self - INTEGER:: k - - DO k = 1, nSpecies - self%n1%output(k)%den = 0.D0 - self%n1%output(k)%mom = 0.D0 - self%n1%output(k)%tensorS = 0.D0 - - self%n2%output(k)%den = 0.D0 - self%n2%output(k)%mom = 0.D0 - self%n2%output(k)%tensorS = 0.D0 - - self%n3%output(k)%den = 0.D0 - self%n3%output(k)%mom = 0.D0 - self%n3%output(k)%tensorS = 0.D0 - - END DO - - END SUBROUTINE resetOutputTria - - - !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS - !Computes element Jacobian determinant - PURE FUNCTION detJCyl(self, xi, dPsi_in) RESULT(dJ) - IMPLICIT NONE - - CLASS(meshVolCyl), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) - REAL(8):: dJ - REAL(8):: dz(1:2), dr(1:2) - - IF(PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dz, dr) - dJ = dz(1)*dr(2)-dz(2)*dr(1) - - END FUNCTION detJCyl - - !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJCyl(self,xi,dPsi_in) RESULT(invJ) - IMPLICIT NONE - - CLASS(meshVolCyl), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) - REAL(8), ALLOCATABLE:: dPsi(:,:) - REAL(8):: dz(1:2), dr(1:2) - REAL(8):: invJ(1:2,1:2) - - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in - - ELSE - dPsi = self%dPsi(xi) - - END IF - - CALL self%partialDer(dPsi, dz, dr) - invJ(1,:) = (/ dr(2), -dz(2) /) - invJ(2,:) = (/ -dr(1), dz(1) /) - - END FUNCTION invJCyl - -END MODULE moduleMeshCyl diff --git a/src/modules/mesh/Cyl/moduleMeshCylBoundary.f95 b/src/modules/mesh/Cyl/moduleMeshCylBoundary.f95 deleted file mode 100644 index 7c29115..0000000 --- a/src/modules/mesh/Cyl/moduleMeshCylBoundary.f95 +++ /dev/null @@ -1,80 +0,0 @@ -!moduleMeshCylBoundary: Edge elements for Cylindrical mesh. -MODULE moduleMeshCylBoundary - USE moduleMeshCyl - - TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylRef - CONTAINS - PROCEDURE, PASS:: fBoundary => reflection - - END TYPE meshEdgeCylRef - - TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAbs - CONTAINS - PROCEDURE, PASS:: fBoundary => absorption - - END TYPE meshEdgeCylAbs - - TYPE, PUBLIC, EXTENDS(meshEdgeCyl):: meshEdgeCylAxis - CONTAINS - PROCEDURE, PASS:: fBoundary => symmetryAxis - - END TYPE meshEdgeCylAxis - - CONTAINS - SUBROUTINE reflection(self, part) - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshEdgeCylRef), INTENT(inout):: self - CLASS(particle), INTENT(inout):: part - REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2) - - edgeNorm = DSQRT((self%r(2)-self%r(1))**2 + (self%z(2)-self%z(1))**2) - cosT = (self%z(2)-self%z(1))/edgeNorm - sinT = DSQRT(1-cosT**2) - - rp(1) = part%r(1) - self%z(1); - rp(2) = part%r(2) - self%r(1); - - rpp(1) = cosT*rp(1) - sinT*rp(2) - rpp(2) = sinT*rp(1) + cosT*rp(2) - rpp(2) = -rpp(2) - - vpp(1) = cosT*part%v(1) - sinT*part%v(2) - vpp(2) = sinT*part%v(1) + cosT*part%v(2) - vpp(2) = -vpp(2) - - part%r(1) = cosT*rpp(1) + sinT*rpp(2) + self%z(1); - part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + self%r(1); - part%v(1) = cosT*vpp(1) + sinT*vpp(2) - part%v(2) = -sinT*vpp(1) + cosT*vpp(2) - - part%n_in = .TRUE. - - END SUBROUTINE reflection - - !Absoption in a surface - SUBROUTINE absorption(self, part) - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshEdgeCylAbs), INTENT(inout):: self - CLASS(particle), INTENT(inout):: part - - - !TODO: Add scatter to mesh nodes - part%n_in = .FALSE. - - END SUBROUTINE absorption - - SUBROUTINE symmetryAxis(self, part) - USE moduleSpecies - IMPLICIT NONE - - CLASS(meshEdgeCylAxis), INTENT(inout):: self - CLASS(particle), INTENT(inout):: part - - END SUBROUTINE symmetryAxis - - -END MODULE moduleMeshCylBoundary diff --git a/src/modules/mesh/Cyl/moduleMeshCylRead.f95 b/src/modules/mesh/Cyl/moduleMeshCylRead.f95 deleted file mode 100644 index 2da8c31..0000000 --- a/src/modules/mesh/Cyl/moduleMeshCylRead.f95 +++ /dev/null @@ -1,610 +0,0 @@ -MODULE moduleMeshCylRead - USE moduleMesh - USE moduleMeshCyl - USE moduleMeshCylBoundary - - !TODO: make this abstract to allow different mesh formats - TYPE, EXTENDS(meshGeneric):: meshCylGeneric - CONTAINS - PROCEDURE, PASS:: init => initCylMesh - PROCEDURE, PASS:: readMesh => readMeshCylGmsh - - END TYPE - - INTERFACE connected - MODULE PROCEDURE connectedVolVol, connectedVolEdge - - END INTERFACE connected - - CONTAINS - !Init mesh - SUBROUTINE initCylMesh(self, meshFormat) - USE moduleMesh - USE moduleErrors - IMPLICIT NONE - - CLASS(meshCylGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - SELECT CASE(meshFormat) - CASE ("gmsh") - self%printOutput => printOutputGmsh - self%printColl => printCollGmsh - self%printEM => printEMGmsh - - CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "initCylMesh") - - END SELECT - - END SUBROUTINE initCylMesh - - !Read mesh from gmsh file - SUBROUTINE readMeshCylGmsh(self, filename) - USE moduleBoundary - IMPLICIT NONE - - CLASS(meshCylGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - REAL(8):: r, z - INTEGER:: p(1:4) - INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0 - INTEGER:: totalNumElem - INTEGER:: boundaryType - - !Read msh - OPEN(10, FILE=TRIM(filename)) - !Skip header - READ(10, *) - READ(10, *) - READ(10, *) - READ(10, *) - !Read number of nodes - READ(10, *) self%numNodes - !Allocate required matrices and vectors - ALLOCATE(self%nodes(1:self%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 - ALLOCATE(meshNodeCyl:: self%nodes(n)%obj) - CALL self%nodes(n)%obj%init(n, (/r, z, 0.D0 /)) - - END DO - !Skips comments - READ(10, *) - READ(10, *) - !Reads Totalnumber of elements - READ(10, *) TotalnumElem - !counts edges and volume elements - self%numEdges = 0 - DO e=1, TotalnumElem - READ(10,*) eTemp, elemType - IF (elemType==1) THEN - self%numEdges=e - END IF - END DO - !Substract the number of edges to the total number of elements - !to obtain the number of volume elements - self%numVols = TotalnumElem - self%numEdges - !Allocates arrays - ALLOCATE(self%edges(1:self%numEdges)) - ALLOCATE(self%vols(1:self%numVols)) - - !Go back to the beggining to read elements - DO e=1, totalNumElem - BACKSPACE(10) - END DO - - !Reads edges - DO e=1, self%numEdges - READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2) - !Associate boundary condition procedure. - bt = getBoundaryId(boundaryType) - SELECT CASE(boundary(bt)%obj%boundaryType) - CASE ('reflection') - ALLOCATE(meshEdgeCylRef:: self%edges(e)%obj) - - CASE ('absorption') - ALLOCATE(meshEdgeCylAbs:: self%edges(e)%obj) - - CASE ('axis') - ALLOCATE(meshEdgeCylAxis:: self%edges(e)%obj) - - END SELECT - - CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) - - END DO - - !Read and initialize volumes - DO e=1, self%numVols - READ(10,*) n, elemType - BACKSPACE(10) - - SELECT CASE(elemType) - CASE (2) - !Triangular element - READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) - ALLOCATE(meshVolCylTria:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3)) - - CASE (3) - !Quadrilateral element - READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4) - ALLOCATE(meshVolCylQuad:: self%vols(e)%obj) - CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) - - END SELECT - - END DO - - CLOSE(10) - - !Build connectivity between elements - DO e = 1, self%numVols - !Connectivity between volumes - DO et = 1, self%numVols - IF (e /= et) THEN - CALL connected(self%vols(e)%obj, self%vols(et)%obj) - - END IF - END DO - - !Connectivity between vols and edges - DO et = 1, self%numEdges - CALL connected(self%vols(e)%obj, self%edges(et)%obj) - - END DO - - !Constructs the global K matrix - CALL constructGlobalK(self%K, self%vols(e)%obj) - - END DO - - END SUBROUTINE readMeshCylGmsh - - !Selects type of elements to build connection - SUBROUTINE connectedVolVol(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshVol), INTENT(inout):: elemB - - SELECT TYPE(elemA) - TYPE IS(meshVolCylQuad) - !Element A is a quadrilateral - SELECT TYPE(elemB) - TYPE IS(meshVolCylQuad) - !Element B is a quadrilateral - CALL connectedQuadQuad(elemA, elemB) - - TYPE IS(meshVolCylTria) - !Element B is a triangle - CALL connectedQuadTria(elemA, elemB) - - END SELECT - - TYPE IS(meshVolCylTria) - !Element A is a Triangle - SELECT TYPE(elemB) - TYPE IS(meshVolCylQuad) - !Element B is a quadrilateral - CALL connectedQuadTria(elemB, elemA) - - TYPE IS(meshVolCylTria) - !Element B is a triangle - CALL connectedTriaTria(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolVol - - SUBROUTINE connectedVolEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: elemA - CLASS(meshEdge), INTENT(inout):: elemB - - SELECT TYPE(elemB) - CLASS IS(meshEdgeCyl) - SELECT TYPE(elemA) - TYPE IS(meshVolCylQuad) - !Element A is a quadrilateral - CALL connectedQuadEdge(elemA, elemB) - - TYPE IS(meshVolCylTria) - !Element A is a triangle - CALL connectedTriaEdge(elemA, elemB) - - END SELECT - - END SELECT - - END SUBROUTINE connectedVolEdge - - SUBROUTINE connectedQuadQuad(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1) .AND. & - elemA%n1%n == elemB%n4%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2) .AND. & - elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n4%n) THEN - elemA%e2 => elemB - elemB%e4 => elemA - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3) .AND. & - elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4) .AND. & - elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - END IF - - END SUBROUTINE connectedQuadQuad - - SUBROUTINE connectedQuadTria(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) THEN - elemA%e2 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n3%n) THEN - elemA%e3 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n3%n == elemB%n3%n .AND. & - elemA%n4%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (elemA%n4%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n3%n) THEN - elemA%e4 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n4%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n4%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e4 => elemB - elemB%e1 => elemA - - END IF - - END IF - - END SUBROUTINE connectedQuadTria - - SUBROUTINE connectedTriaTria(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n3%n) THEN - elemA%e1 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n1%n == elemB%n3%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n3%n) THEN - elemA%e2 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n2%n == elemB%n3%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - END IF - - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n3%n) THEN - elemA%e3 => elemB - elemB%e3 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - ELSEIF (elemA%n3%n == elemB%n3%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - END IF - - - END IF - - END SUBROUTINE connectedTriaTria - - SUBROUTINE connectedQuadEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n4%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n4%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 4 - IF (.NOT. ASSOCIATED(elemA%e4)) THEN - IF (elemA%n4%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e4 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n4%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e4 => elemB - elemB%e1 => elemA - - END IF - - END IF - - END SUBROUTINE connectedQuadEdge - - SUBROUTINE connectedTriaEdge(elemA, elemB) - IMPLICIT NONE - - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA - CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB - - !Check direction 1 - IF (.NOT. ASSOCIATED(elemA%e1)) THEN - IF (elemA%n1%n == elemB%n1%n .AND. & - elemA%n2%n == elemB%n2%n) THEN - elemA%e1 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n1%n == elemB%n2%n .AND. & - elemA%n2%n == elemB%n1%n) THEN - elemA%e1 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 2 - IF (.NOT. ASSOCIATED(elemA%e2)) THEN - IF (elemA%n2%n == elemB%n1%n .AND. & - elemA%n3%n == elemB%n2%n) THEN - elemA%e2 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n2%n == elemB%n2%n .AND. & - elemA%n3%n == elemB%n1%n) THEN - elemA%e2 => elemB - elemB%e1 => elemA - - END IF - - END IF - - !Check direction 3 - IF (.NOT. ASSOCIATED(elemA%e3)) THEN - IF (elemA%n3%n == elemB%n1%n .AND. & - elemA%n1%n == elemB%n2%n) THEN - elemA%e3 => elemB - elemB%e2 => elemA - - ELSEIF (elemA%n3%n == elemB%n2%n .AND. & - elemA%n1%n == elemB%n1%n) THEN - elemA%e3 => elemB - elemB%e1 => elemA - - END IF - - END IF - - END SUBROUTINE connectedTriaEdge - - 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 /) - - TYPE IS(meshVolCylTria) - nNodes = 3 - ALLOCATE(localK(1:nNodes,1:nNodes)) - localK = elem%elemK() - ALLOCATE(n(1:nNodes)) - n = (/ elem%n1%n, elem%n2%n, elem%n3%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 - -END MODULE moduleMeshCylRead diff --git a/src/modules/mesh/makefile b/src/modules/mesh/makefile index 19e3a13..67a7629 100644 --- a/src/modules/mesh/makefile +++ b/src/modules/mesh/makefile @@ -1,11 +1,11 @@ -all: moduleMesh.o Cyl.o 1DCart.o +all: moduleMesh.o 2DCyl.o 1DCart.o -Cyl.o: - $(MAKE) -C Cyl all +2DCyl.o: + $(MAKE) -C 2DCyl all 1DCart.o: $(MAKE) -C 1DCart all -moduleMesh.o: moduleMesh.f95 - $(FC) $(FCFLAGS) -c $(subst .o,.f95,$@) -o $(OBJDIR)/$@ +moduleMesh.o: moduleMesh.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 285c0a1..c858724 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -240,7 +240,6 @@ MODULE moduleMesh PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() - CONTAINS PROCEDURE(initMesh_interface), DEFERRED, PASS:: init PROCEDURE(readMesh_interface), DEFERRED, PASS:: readMesh diff --git a/src/modules/mesh/moduleMesh.f95 b/src/modules/mesh/moduleMesh.f95 deleted file mode 100644 index 9e562f4..0000000 --- a/src/modules/mesh/moduleMesh.f95 +++ /dev/null @@ -1,606 +0,0 @@ -!moduleMesh: General module for Finite Element mesh -MODULE moduleMesh - USE moduleList - USE moduleOutput - IMPLICIT NONE - - !Parent of Node element - TYPE, PUBLIC, ABSTRACT:: meshNode - !Node index - INTEGER:: n = 0 - !Node volume - 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 - - END TYPE meshNode - - ABSTRACT INTERFACE - !Interface of init a node (3D generic coordinates) - SUBROUTINE initNode_interface(self, n, r) - IMPORT:: meshNode - CLASS(meshNode), INTENT(out):: self - INTEGER, INTENT(in):: n - REAL(8), INTENT(in):: r(1:3) - - END SUBROUTINE initNode_interface - - !Interface to get coordinates from node - PURE FUNCTION getCoord_interface(self) RESULT(r) - IMPORT:: meshNode - CLASS(meshNode), INTENT(in):: self - REAL(8):: r(1:3) - - END FUNCTION getCoord_interface - - END INTERFACE - - !Containers for nodes in the mesh - TYPE:: meshNodeCont - CLASS(meshNode), ALLOCATABLE:: obj - - END TYPE meshNodeCont - - !Parent of Edge element - TYPE, PUBLIC, ABSTRACT:: meshEdge - !Element index - INTEGER:: n = 0 - !Connectivity to vols - CLASS(meshVol), POINTER:: e1 => NULL(), e2 => NULL() - !Normal vector - REAL(8):: normal(1:3) - !Physical surface in mesh - INTEGER:: physicalSurface - !id for boundary condition - INTEGER:: bt = 0 - CONTAINS - 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(:) - INTEGER, INTENT(in):: bt - INTEGER, INTENT(in):: physicalSurface - - END SUBROUTINE initEdge_interface - - SUBROUTINE boundary_interface(self, part) - USE moduleSpecies - - IMPORT:: meshEdge - CLASS (meshEdge), INTENT(inout):: self - CLASS (particle), INTENT(inout):: part - - END SUBROUTINE - - 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):: r(1:3) - - END FUNCTION randPos_interface - - END INTERFACE - - !Containers for edges in the mesh - TYPE:: meshEdgeCont - CLASS(meshEdge), ALLOCATABLE:: obj - - END TYPE meshEdgeCont - - !Parent of Volume element - TYPE, PUBLIC, ABSTRACT:: meshVol - !Volume index - INTEGER:: n = 0 - !Maximum collision rate - REAL(8):: sigmaVrelMax = 1.D-15 - !Volume - REAL(8):: volume = 0.D0 - !List of particles inside the volume - TYPE(listNode):: listPart_in - !Lock indicator for listPart_in - INTEGER(KIND=OMP_LOCK_KIND):: lock - !Number of collisions per volume - INTEGER:: nColl = 0 - !Total weight of particles inside cell - REAL(8):: totalWeight = 0.D0 - 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, PASS:: findCell - PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log - PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside - PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement - PROCEDURE, PASS:: collision - PROCEDURE(resetOutput_interface), DEFERRED, PASS:: resetOutput - - END TYPE meshVol - - ABSTRACT INTERFACE - SUBROUTINE initVol_interface(self, n, p) - IMPORT:: meshVol - CLASS(meshVol), INTENT(out):: self - INTEGER, INTENT(in):: n - INTEGER, INTENT(in):: p(:) - - 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 - 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 nextElement_interface(self, xi, nextElement) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - CLASS(*), POINTER, INTENT(out):: nextElement - - END SUBROUTINE nextElement_interface - - PURE FUNCTION phy2log_interface(self,r) RESULT(xN) - IMPORT:: meshVol - CLASS(meshVol), INTENT(in):: self - REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - - END FUNCTION phy2log_interface - - PURE FUNCTION inside_interface(xi) RESULT(ins) - IMPORT:: meshVol - REAL(8), INTENT(in):: xi(1:3) - LOGICAL:: ins - - END FUNCTION inside_interface - - SUBROUTINE collision_interface(self) - IMPORT:: meshVol - CLASS(meshVol), INTENT(inout):: self - - END SUBROUTINE collision_interface - - PURE SUBROUTINE resetOutput_interface(self) - IMPORT:: meshVol - CLASS(meshVol), INTENT(inout):: self - - END SUBROUTINE resetOutput_interface - - END INTERFACE - - !Containers for volumes in the mesh - TYPE:: meshVolCont - CLASS(meshVol), ALLOCATABLE:: obj - - END TYPE meshVolCont - - !Abstract type of mesh - TYPE, PUBLIC, ABSTRACT:: meshGeneric - INTEGER:: numEdges, numNodes, numVols - !Array of nodes - TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) - !Array of boundary elements - TYPE(meshEdgeCont), ALLOCATABLE:: edges(:) - !Array of volume elements - TYPE(meshVolCont), ALLOCATABLE:: vols(:) - !Global stiffness matrix - REAL(8), ALLOCATABLE, DIMENSION(:,:):: K - !Permutation matrix for P L U factorization - INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV - PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() - PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() - PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() - - CONTAINS - PROCEDURE(initMesh_interface), DEFERRED, PASS:: init - PROCEDURE(readMesh_interface), DEFERRED, PASS:: readMesh - - END TYPE meshGeneric - - ABSTRACT INTERFACE - !Inits the mesh - SUBROUTINE initMesh_interface(self, meshFormat) - IMPORT meshGeneric - - CLASS(meshGeneric), INTENT(out):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat - - END SUBROUTINE initMesh_interface - - !Reads the mesh from a file - SUBROUTINE readMesh_interface(self, filename) - IMPORT meshGeneric - - CLASS(meshGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: filename - - END SUBROUTINE readMesh_interface - - !Prints Species data - SUBROUTINE printOutput_interface(self, t) - IMPORT meshGeneric - - CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - - END SUBROUTINE printOutput_interface - - !Prints number of collisions - SUBROUTINE printColl_interface(self, t) - IMPORT meshGeneric - - CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - - 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 - CLASS(meshGeneric), ALLOCATABLE, TARGET:: mesh - - CONTAINS - !Find next cell for particle - RECURSIVE SUBROUTINE findCell(self, part, oldCell) - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol), INTENT(inout):: self - CLASS(meshVol), OPTIONAL, INTENT(in):: oldCell - CLASS(particle), INTENT(inout), TARGET:: part - REAL(8):: xi(1:3) - CLASS(*), POINTER:: nextElement - - xi = self%phy2log(part%r) - !Checks if particle is inside 'self' cell - IF (self%inside(xi)) THEN - part%vol = self%n - part%xi = xi - part%n_in = .TRUE. - !Assign particle to listPart_in - CALL OMP_SET_LOCK(self%lock) - CALL self%listPart_in%add(part) - self%totalWeight = self%totalWeight + part%weight - CALL OMP_UNSET_LOCK(self%lock) - - ELSE - !If not, searches for a neighbour and repeats the process. - CALL self%nextElement(xi, nextElement) - !Defines the next step - SELECT TYPE(nextElement) - CLASS IS(meshVol) - !Particle moved to new cell, repeat find procedure - CALL nextElement%findCell(part, self) - - CLASS IS (meshEdge) - !Particle encountered an edge, execute boundary - CALL nextElement%fBoundary(part) - !If particle is still inside the domain, call findCell - IF (part%n_in) THEN - IF(PRESENT(oldCell)) THEN - CALL self%findCell(part, oldCell) - - ELSE - CALL self%findCell(part) - - END IF - END IF - - CLASS DEFAULT - WRITE(*,*) "ERROR, CHECK findCell" - - END SELECT - END IF - - END SUBROUTINE findCell - - !Computes collisions in element - SUBROUTINE collision(self) - USE moduleCollisions - USE moduleSpecies - USE moduleList - use moduleRefParam - 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 - TYPE(particle), POINTER:: part_i, part_j - INTEGER:: n !collision - INTEGER:: ij, k - REAL(8):: sigmaVrelMaxNew - TYPE(pointerArray), ALLOCATABLE:: partTemp(:) - - self%nColl = 0 - nPart = self%listPart_in%amount - IF (nPart > 1) THEN - pMax = self%totalWeight*self%sigmaVrelMax*tauMin/self%volume - self%nColl = INT(REAL(nPart)*pMax*0.5D0) - - !Converts the list of particles to an array for easy access - IF (self%nColl > 0) THEN - partTemp = self%listPart_in%convert2Array() - - END IF - - DO n = 1, self%nColl - !Select random numbers - rnd = 1 + FLOOR(nPart*RAND()) - part_i => partTemp(rnd)%part - rnd = 1 + FLOOR(nPart*RAND()) - part_j => partTemp(rnd)%part - 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) - - END DO - - !Update maximum cross section*v_rel per each collision - IF (sigmaVrelMaxNew > self%sigmaVrelMax) THEN - self%sigmaVrelMax = sigmaVrelMaxNew - - END IF - - END DO - - END IF - - self%totalWeight = 0.D0 - - !Reset output in nodes - CALL self%resetOutput() - - !Erase the list of particles inside the cell - CALL self%listPart_in%erase() - - END SUBROUTINE collision - - SUBROUTINE printOutputGmsh(self, t) - USE moduleRefParam - USE moduleSpecies - USE moduleOutput - IMPLICIT NONE - - CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - INTEGER:: n, i - TYPE(outputFormat):: output(1:self%numNodes) - REAL(8):: time - CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations - - time = DBLE(t)*tauMin*ti_ref - - DO i = 1, nSpecies - WRITE(tstring, '(I6.6)') t - fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh' - WRITE(*, "(6X,A15,A)") "Creating file: ", fileName - OPEN (60, file = path // folder // '/' // fileName) - WRITE(60, "(A)") '$MeshFormat' - WRITE(60, "(A)") '2.2 0 8' - WRITE(60, "(A)") '$EndMeshFormat' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Density (m^-3)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj) - WRITE(60, "(I6,ES20.6E3)") n, output(n)%density - END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Velocity (m/s)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 3 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity - END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Pressure (Pa)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure - END DO - WRITE(60, "(A)") '$EndNodeData' - WRITE(60, "(A)") '$NodeData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Temperature (K)"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numNodes - DO n=1, self%numNodes - WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature - END DO - WRITE(60, "(A)") '$EndNodeData' - CLOSE (60) - - END DO - - END SUBROUTINE printOutputGmsh - - SUBROUTINE printCollGmsh(self, t) - USE moduleRefParam - USE moduleCaseParam - USE moduleOutput - IMPLICIT NONE - - CLASS(meshGeneric), INTENT(in):: self - INTEGER, INTENT(in):: t - INTEGER:: n - REAL(8):: time - CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations - - - IF (collOutput) THEN - time = DBLE(t)*tauMin*ti_ref - WRITE(tstring, '(I6.6)') t - - fileName='OUTPUT_' // tstring// '_Collisions.msh' - WRITE(*, "(6X,A15,A)") "Creating file: ", fileName - OPEN (60, file = path // folder // '/' // fileName) - WRITE(60, "(A)") '$MeshFormat' - WRITE(60, "(A)") '2.2 0 8' - WRITE(60, "(A)") '$EndMeshFormat' - WRITE(60, "(A)") '$ElementData' - WRITE(60, "(A)") '1' - WRITE(60, "(A)") '"Collisions"' - WRITE(60, *) 1 - WRITE(60, *) time - WRITE(60, *) 3 - WRITE(60, *) t - WRITE(60, *) 1 - WRITE(60, *) self%numVols - DO n=1, self%numVols - WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl - END DO - WRITE(60, "(A)") '$EndElementData' - - CLOSE(60) - - END IF - - END SUBROUTINE printCollGmsh - - SUBROUTINE printEMGmsh(self, t) - USE moduleRefParam - USE moduleCaseParam - USE moduleOutput - IMPLICIT NONE - - CLASS(meshGeneric), 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 - REAL(8):: xi(1:3) - - xi = (/ 0.D0, 0.D0, 0.D0 /) - - IF (emOutput) THEN - time = DBLE(t)*tauMin*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(xi)*EF_ref - END DO - WRITE(20, "(A)") '$EndElementData' - CLOSE(20) - - END IF - - END SUBROUTINE printEMGmsh - -END MODULE moduleMesh diff --git a/src/modules/moduleBoundary.f95 b/src/modules/moduleBoundary.f95 deleted file mode 100644 index 71caa16..0000000 --- a/src/modules/moduleBoundary.f95 +++ /dev/null @@ -1,36 +0,0 @@ -MODULE moduleBoundary - - TYPE, PUBLIC:: boundaryGeneric - INTEGER:: id = 0 - CHARACTER(:), ALLOCATABLE:: name - INTEGER:: physicalSurface = 0 - CHARACTER(:), ALLOCATABLE:: boundaryType - - END TYPE boundaryGeneric - - TYPE:: boundaryCont - CLASS(boundaryGeneric), ALLOCATABLE:: obj - - END TYPE boundaryCont - - - INTEGER:: nBoundary = 0 - TYPE(boundaryCont), ALLOCATABLE:: boundary(:) - - CONTAINS - FUNCTION getBoundaryId(physicalSurface) RESULT(id) - IMPLICIT NONE - - INTEGER:: physicalSurface - INTEGER:: id - INTEGER:: i - - id = 0 - DO i = 1, nBoundary - IF (physicalSurface == boundary(i)%obj%physicalSurface) id = boundary(i)%obj%id - - END DO - - END FUNCTION getBoundaryId - -END MODULE moduleBoundary diff --git a/src/modules/moduleCaseParam.f95 b/src/modules/moduleCaseParam.f95 deleted file mode 100644 index 5f1023a..0000000 --- a/src/modules/moduleCaseParam.f95 +++ /dev/null @@ -1,9 +0,0 @@ -!Problems of the case -MODULE moduleCaseParam - !Maximum number of iterations and number of species - INTEGER:: tmax - REAL(8), ALLOCATABLE:: tau(:) - REAL(8):: tauMin - -END MODULE moduleCaseParam - diff --git a/src/modules/moduleCollisions.f95 b/src/modules/moduleCollisions.f95 deleted file mode 100644 index 59fcbbe..0000000 --- a/src/modules/moduleCollisions.f95 +++ /dev/null @@ -1,171 +0,0 @@ -MODULE moduleCollisions - USE moduleTable - - TYPE, ABSTRACT:: collisionBinary - REAL(8):: rMass !reduced mass - TYPE(table1D):: crossSec !cross section of collision - CONTAINS - PROCEDURE(initBinary_interface), PASS, DEFERRED:: init - PROCEDURE(collideBinary_interface), PASS, DEFERRED:: collide - - END TYPE collisionBinary - - ABSTRACT INTERFACE - SUBROUTINE initBinary_interface(self, crossSectionFilename, mass_i, mass_j) - IMPORT:: collisionBinary - CLASS(collisionBinary), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename - REAL(8), INTENT(in):: mass_i, mass_j - - END SUBROUTINE - - SUBROUTINE collideBinary_interface(self, sigmaVRelMax, sigmaVrelMaxNew, part_i, part_j) - USE moduleSpecies - IMPORT:: collisionBinary - - CLASS(collisionBinary), INTENT(in):: self - REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(inout):: sigmaVrelMaxNew - TYPE(particle), INTENT(inout):: part_i, part_j - - END SUBROUTINE - - END INTERFACE - - !Container for binary collisions - TYPE:: collisionCont - CLASS(collisionBinary), ALLOCATABLE:: obj - - END TYPE collisionCont - - TYPE, EXTENDS(collisionBinary):: collisionBinaryElastic - !Weight distribution for Maxwellian function - REAL(8):: w_i = (1.D0+DSQRT(3.D0))/2.D0 - REAL(8):: w_j = (DSQRT(3.D0)-1.D0)/2.D0 - CONTAINS - PROCEDURE, PASS:: init => initBinaryElastic - PROCEDURE, PASS:: collide => collideBinaryElastic - - END TYPE collisionBinaryElastic - - !Type for interaction matrix - TYPE:: interactionsBinary - INTEGER:: amount - TYPE(collisionCont), ALLOCATABLE:: collisions(:) - CONTAINS - PROCEDURE, PASS:: init => initInteractionBinary - - END TYPE interactionsBinary - - !Collision 'Matrix'. A symmetric 2D matrix put into a 1D array to save memory - TYPE(interactionsBinary), ALLOCATABLE:: interactionMatrix(:) - !Folder for collision cross section tables - CHARACTER(:), ALLOCATABLE:: pathCollisions - - CONTAINS - SUBROUTINE initInteractionMatrix(interactionMatrix) - USE moduleSpecies - IMPLICIT NONE - - TYPE(interactionsBinary), INTENT(inout), ALLOCATABLE:: interactionMatrix(:) - INTEGER:: nInteractions - - nInteractions = (nSpecies*(nSpecies+1))/2 - ALLOCATE(interactionMatrix(1:nInteractions)) - - END SUBROUTINE initInteractionMatrix - - FUNCTION interactionIndex(i,j) RESULT(k) - - INTEGER:: i, j - INTEGER:: p - INTEGER:: k - - k = i + j - p = (k + ABS(i - j))/2 - k = k + (p*(p-3))/2 - - END FUNCTION interactionIndex - - SUBROUTINE initInteractionBinary(self, amount) - IMPLICIT NONE - - CLASS(interactionsBinary), INTENT(inout):: self - INTEGER, INTENT(in):: amount - INTEGER:: k - - self%amount = amount - - ALLOCATE(self%collisions(1:self%amount)) - DO k= 1, self%amount - !TODO: make type dependent - ALLOCATE(collisionBinaryElastic:: self%collisions(k)%obj) - - END DO - - END SUBROUTINE initInteractionBinary - - SUBROUTINE initBinaryElastic(self, crossSectionFilename, mass_i, mass_j) - USE moduleTable - USE moduleRefParam - USE moduleConstParam - IMPLICIT NONE - - CLASS(collisionBinaryElastic), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSectionFilename - REAL(8), INTENT(in):: mass_i, mass_j - - !Reads data from file - CALL self%crossSec%init(crossSectionFilename) - - !Convert to no-dimensional units - CALL self%crossSec%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2) - - !Calculates reduced mass - self%rMass = (mass_i*mass_j)/(mass_i+mass_j) - - END SUBROUTINE - - SUBROUTINE collideBinaryElastic(self, sigmaVrelMax, sigmaVrelMaxNew, & - part_i, part_j) - USE moduleConstParam - USE moduleSpecies - USE moduleTable - IMPLICIT NONE - - CLASS(collisionBinaryElastic), INTENT(in):: self - REAL(8), INTENT(in):: sigmaVrelMax - REAL(8), INTENT(inout):: sigmaVrelMaxNew - TYPE(particle), INTENT(inout):: part_i, part_j - REAL(8):: sigmaVrel - REAL(8):: vRel !relative velocity - REAL(8):: eRel !relative energy - 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 - - !eRel (in units of [m][L]^2[s]^-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 - IF (sigmaVrelMaxNew/sigmaVrelMax > RAND()) 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 - alpha = PI*RAND() - part_i%v(1) = v_i*DCOS(alpha) - part_i%v(2) = v_i*DSIN(alpha) - alpha = PI*RAND() - part_j%v(1) = v_j*DCOS(alpha) - part_j%v(2) = v_j*DSIN(alpha) - - END IF - - - END SUBROUTINE - -END MODULE moduleCollisions diff --git a/src/modules/moduleCompTime.f95 b/src/modules/moduleCompTime.f95 deleted file mode 100644 index 28273ce..0000000 --- a/src/modules/moduleCompTime.f95 +++ /dev/null @@ -1,11 +0,0 @@ -!Information to calculate computation time -MODULE moduleCompTime - REAL(8):: tStep=0.D0 - REAL(8):: tPush=0.D0 - 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 deleted file mode 100644 index 634f845..0000000 --- a/src/modules/moduleConstParam.f95 +++ /dev/null @@ -1,15 +0,0 @@ -!Physical and mathematical constants -MODULE moduleConstParam - IMPLICIT NONE - - PUBLIC - - 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:: qe = 1.60217662D-19 !Elementary charge - REAL(8), PARAMETER:: kb = 1.38064852D-23 !Boltzmann constants SI - 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 deleted file mode 100644 index 876093e..0000000 --- a/src/modules/moduleEM.f95 +++ /dev/null @@ -1,128 +0,0 @@ -!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%xi - - elField = mesh%vols(part%vol)%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 DO REDUCTION(+:vectorF) - 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 - - !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) - DEALLOCATE(nodes, rho) - - END DO - !$OMP END DO - - !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 deleted file mode 100644 index 0a95db8..0000000 --- a/src/modules/moduleErrors.f95 +++ /dev/null @@ -1,43 +0,0 @@ -!moduleErrors: Manages the different type of errors the program can produce. -! By errors we understand critical errors (that stop the program), -! warnings (that only output a message with a WARNING tag), -! o verbose outputs that can be used for debugging porpouse. -MODULE moduleErrors - CONTAINS - SUBROUTINE criticalError(msg, pgr) - IMPLICIT NONE - - CHARACTER(*), INTENT(in):: msg, pgr - CHARACTER(:), ALLOCATABLE:: errorMsg - - errorMsg = "CRITICAL error in " // pgr // " with message:" // NEW_LINE('A') // msg - - ERROR STOP errorMsg - - END SUBROUTINE criticalError - - SUBROUTINE warningError(msg) - IMPLICIT NONE - - CHARACTER(*), INTENT(in):: msg - CHARACTER(:), ALLOCATABLE:: errorMsg - - errorMsg = "WARNING: " // msg - - WRITE (*, '(A)') errorMsg - - END SUBROUTINE warningError - - SUBROUTINE verboseError(msg) - IMPLICIT NONE - - CHARACTER(*), INTENT(in):: msg - CHARACTER(:), ALLOCATABLE:: errorMsg - - errorMsg = msg - - WRITE (*, '(A)') errorMsg - - END SUBROUTINE verboseError - -END MODULE moduleErrors diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index eb3c737..47edbf0 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -1,21 +1,63 @@ !injection of particles MODULE moduleInject + !Generic type for velocity distribution function + TYPE, ABSTRACT:: velDistGeneric + CONTAINS + !Returns random velocity from distribution function + PROCEDURE(randomVel_interface), DEFERRED, PASS:: randomVel + + END TYPE velDistGeneric + + ABSTRACT INTERFACE + FUNCTION randomVel_interface(self) RESULT(v) + IMPORT velDistGeneric + + CLASS(velDistGeneric), INTENT(in):: self + REAL(8):: v + + END FUNCTION randomVel_interface + + END INTERFACE + + !Container for velocity distributions + TYPE:: velDistCont + CLASS(velDistGeneric), ALLOCATABLE:: obj + + END TYPE velDistCont + + !Maxwellian distribution function + TYPE, EXTENDS(velDistGeneric):: velDistMaxwellian + REAL(8):: v !Velocity + REAL(8):: vTh !Thermal Velocity + CONTAINS + PROCEDURE, PASS:: randomVel => randomVelMaxwellian + + END TYPE velDistMaxwellian + + !Dirac's delta distribution function + TYPE, EXTENDS(velDistGeneric):: velDistDelta + REAL(8):: v !Velocity + CONTAINS + PROCEDURE, PASS:: randomVel => randomVelDelta + + END TYPE velDistDelta + + !Generic injection of particles TYPE:: injectGeneric INTEGER:: id CHARACTER(:), ALLOCATABLE:: name REAL(8):: vMod !Velocity (module) REAL(8):: T(1:3) !Temperature - REAL(8):: v(1:3) !Velocity(vector) - 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:: sp !Species of injection INTEGER:: nEdges INTEGER, ALLOCATABLE:: edges(:) !Array with edges + TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction CONTAINS PROCEDURE, PASS:: init => initInject - PROCEDURE, PASS:: addParticles => addParticlesMaxwellian + PROCEDURE, PASS:: addParticles END TYPE injectGeneric @@ -64,9 +106,6 @@ MODULE moduleInject self%nParticles = self%nParticles * solver%pusher(sp)%every self%sp = sp - self%v = self%vMod * self%n - self%vTh = DSQRT(self%T/species(self%sp)%obj%m) - !Gets the edge elements from which particles are injected !TODO: Improve this A LOT DO e = 1, mesh%numEdges @@ -88,7 +127,7 @@ MODULE moduleInject END DO ! self%sumWeight = SUM(self%weight) - END SUBROUTINE + END SUBROUTINE initInject !Injection of particles SUBROUTINE doInjects() @@ -114,12 +153,35 @@ MODULE moduleInject END SUBROUTINE doInjects - !Random velocity from maxwellian distribution - REAL(8) FUNCTION vBC(u, vth) + SUBROUTINE initVelDistMaxwellian(velDist, v, T, m) + IMPLICIT NONE + + CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist + REAL(8), INTENT(in):: v, T, m + + velDist = velDistMaxwellian(v = v, vTh = DSQRT(T/m)) + + END SUBROUTINE initVelDistMaxwellian + + SUBROUTINE initVelDistDelta(velDist, v) + IMPLICIT NONE + + CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist + REAL(8), INTENT(in):: v + + velDist = velDistDelta(v = v) + + END SUBROUTINE initVelDistDelta + + !Random velocity from Maxwellian distribution + FUNCTION randomVelMaxwellian(self) RESULT (v) USE moduleConstParam, ONLY: PI - REAL(8), INTENT (in):: u, vth + IMPLICIT NONE + + CLASS(velDistMaxwellian), INTENT(in):: self + REAL(8):: v REAL(8):: x, y - vBC=0.D0 + v = 0.D0 x = 0.D0 DO WHILE (x == 0.D0) @@ -127,11 +189,23 @@ MODULE moduleInject END DO CALL RANDOM_NUMBER(y) - vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) + v = self%v + self%vTh*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) - END FUNCTION vBC + END FUNCTION randomVelMaxwellian - SUBROUTINE addParticlesMaxwellian(self) + !Random velocity from Dirac's delta distribution + PURE FUNCTION randomVelDelta(self) RESULT(v) + IMPLICIT NONE + + CLASS(velDistDelta), INTENT(in):: self + REAL(8):: v + + v = self%v + + END FUNCTION randomVelDelta + + !Add particles for the injection + SUBROUTINE addParticles(self) USE moduleSpecies USE moduleSolver USE moduleMesh @@ -185,9 +259,9 @@ MODULE moduleInject END IF - partInj(n)%v = (/ vBC(self%v(1), self%vTh(1)), & - vBC(self%v(2), self%vTh(2)), & - vBC(self%v(3), self%vTh(3)) /) + partInj(n)%v = (/ self%v(1)%obj%randomVel(), & + self%v(2)%obj%randomVel(), & + self%v(3)%obj%randomVel() /) !Push new particle CALL solver%pusher(self%sp)%pushParticle(partInj(n)) @@ -197,6 +271,6 @@ MODULE moduleInject END DO !$OMP END DO - END SUBROUTINE addParticlesMaxwellian + END SUBROUTINE addParticles END MODULE moduleInject diff --git a/src/modules/moduleInject.f95 b/src/modules/moduleInject.f95 deleted file mode 100644 index c524257..0000000 --- a/src/modules/moduleInject.f95 +++ /dev/null @@ -1,202 +0,0 @@ -!injection of particles -MODULE moduleInject - - TYPE:: injectGeneric - INTEGER:: id - CHARACTER(:), ALLOCATABLE:: name - REAL(8):: vMod !Velocity (module) - REAL(8):: T(1:3) !Temperature - REAL(8):: v(1:3) !Velocity(vector) - 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:: sp !Species of injection - INTEGER:: nEdges - INTEGER, ALLOCATABLE:: edges(:) !Array with edges - ! REAL(8), ALLOCATABLE:: weight(:) !weight of cells for injection - ! REAL(8):: sumWeight - CONTAINS - PROCEDURE, PASS:: init => initInject - PROCEDURE, PASS:: addParticles => addParticlesMaxwellian - - END TYPE injectGeneric - - INTEGER:: nInject - TYPE(injectGeneric), ALLOCATABLE:: inject(:) - - CONTAINS - !Initialize an injection of particles - SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface) - USE moduleMesh - USE moduleRefParam - USE moduleConstParam - USE moduleSpecies - USE moduleSolver - USE moduleErrors - IMPLICIT NONE - - CLASS(injectGeneric), INTENT(inout):: self - INTEGER, INTENT(in):: i - REAL(8), INTENT(in):: v, n(1:3), T(1:3) - INTEGER, INTENT(in):: sp, physicalSurface - REAL(8), INTENT(in):: flow - CHARACTER(:), ALLOCATABLE, INTENT(in):: units - INTEGER:: e, et - INTEGER:: phSurface(1:mesh%numEdges) - - self%id = i - self%vMod = v/v_ref - self%n = n - self%T = T/T_ref - SELECT CASE(units) - CASE ("sccm") - !Standard cubic centimeter per minute - self%nParticles = INT(flow*sccm2atomPerS*tauMin*ti_ref/species(sp)%obj%weight) - - CASE ("A") - !Input current in Ampers - self%nParticles = INT(flow*tauMin*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%nParticles = self%nParticles * solver%pusher(sp)%every - self%sp = sp - - self%v = self%vMod * self%n - self%vTh = DSQRT(self%T/species(self%sp)%obj%m) - - !Gets the edge elements from which particles are injected - !TODO: Improve this A LOT - DO e = 1, mesh%numEdges - phSurface(e) = mesh%edges(e)%obj%physicalSurface - - END DO - - self%nEdges = COUNT(phSurface == physicalSurface) - ALLOCATE(inject(i)%edges(1:self%nEdges)) - ! ALLOCATE(inject(i)%weight(1:self%nEdges)) - et = 0 - DO e=1, mesh%numEdges - IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN - et = et + 1 - self%edges(et) = mesh%edges(e)%obj%n - - END IF - - END DO - ! self%sumWeight = SUM(self%weight) - - END SUBROUTINE - - !Injection of particles - SUBROUTINE doInjects() - USE moduleSpecies - USE moduleSolver - IMPLICIT NONE - - INTEGER:: i - - !$OMP SINGLE - nPartInj = 0 - DO i = 1, nInject - IF (solver%pusher(inject(i)%sp)%pushSpecies) nPartInj = nPartInj + inject(i)%nParticles - - END DO - IF (ALLOCATED(partInj)) DEALLOCATE(partInj) - ALLOCATE(partInj(1:nPartInj)) - !$OMP END SINGLE - - DO i=1, nInject - IF (solver%pusher(inject(i)%sp)%pushSpecies) CALL inject(i)%addParticles() - END DO - - END SUBROUTINE doInjects - - !Random velocity from maxwellian distribution - REAL(8) FUNCTION vBC(u, vth) - USE moduleConstParam, ONLY: PI - REAL(8), INTENT (in):: u, vth - REAL(8):: x, y - vBC=0.D0 - x = 0.D0 - - DO WHILE (x == 0.D0) - x=RAND() - END DO - y=RAND() - - vBC = u + vth*DSQRT(-2.D0*DLOG(x))*DCOS(2.D0*PI*y) - - END FUNCTION vBC - - SUBROUTINE addParticlesMaxwellian(self) - USE moduleSpecies - USE moduleSolver - USE moduleMesh - IMPLICIT NONE - - CLASS(injectGeneric), INTENT(in):: self - INTEGER:: randomX - INTEGER:: i!, j - INTEGER, SAVE:: nMin, nMax !Min and Max index in partInj array - INTEGER:: n - CLASS(meshEdge), POINTER:: randomEdge - - !Insert particles - !$OMP SINGLE - nMin = 0 - DO i = 1, self%id - 1 - IF (solver%pusher(inject(i)%sp)%pushSpecies) nMin = nMin + inject(i)%nParticles - - END DO - nMin = nMin + 1 - nMax = nMin + self%nParticles - 1 - !Assign particle type - partInj(nMin:nMax)%sp = self%sp - !Assign weight to particle. - 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%sp)%obj) - TYPE IS(speciesCharged) - partInj(nMin:nMax)%qm = sp%qm - - END SELECT - !$OMP END SINGLE - - !$OMP DO - DO n = nMin, nMax - randomX = INT(DBLE(self%nEdges-1)*RAND()) + 1 - - randomEdge => mesh%edges(self%edges(randomX))%obj - !Random position in edge - partInj(n)%r = randomEdge%randPos() - !Volume associated to the edge: - IF (DOT_PRODUCT(self%n, randomEdge%normal) >= 0.D0) THEN - partInj(n)%vol = randomEdge%e1%n - - ELSE - partInj(n)%vol = randomEdge%e2%n - - END IF - - partInj(n)%v = (/ vBC(self%v(1), self%vTh(1)), & - vBC(self%v(2), self%vTh(2)), & - vBC(self%v(3), self%vTh(3)) /) - - !Push new particle - CALL solver%pusher(self%sp)%pushParticle(partInj(n)) - !Assign cell to new particle - CALL solver%updateParticleCell(partInj(n)) - - END DO - !$OMP END DO - - END SUBROUTINE addParticlesMaxwellian - -END MODULE moduleInject diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index d31b5cf..e084352 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -520,6 +520,8 @@ MODULE moduleInput CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface) + CALL readVelDistr(config, inject(i), object) + END DO !Allocate array for injected particles @@ -530,6 +532,52 @@ MODULE moduleInput END SUBROUTINE readInject + !Reads the velocity distribution functions for each inject + SUBROUTINE readVelDistr(config, inj, object) + USE moduleErrors + USE moduleInject + USE moduleSpecies + USE json_module + IMPLICIT NONE + + TYPE(json_file), INTENT(inout):: config + TYPE(injectGeneric), INTENT(inout):: inj + CHARACTER(:), ALLOCATABLE, INTENT(in):: object + INTEGER:: i + CHARACTER(2):: istring + CHARACTER(:), ALLOCATABLE:: fvType + LOGICAL:: found + REAL(8):: v, T, m + + !Reads species mass + m = species(inj%sp)%obj%m + !Reads distribution functions for velocity + DO i = 1, 3 + WRITE(istring, '(i2)') i + CALL config%get(object // '.velDist('// TRIM(istring) //')', fvType, found) + IF (.NOT. found) CALL criticalError("No velocity distribution in direction " // istring // & + " found for " // object, 'readVelDistr') + + SELECT CASE(fvType) + CASE ("Maxwellian") + v = inj%vMod*inj%n(i) + T = inj%T(i) + CALL initVelDistMaxwellian(inj%v(i)%obj, v, t, m) + + CASE ("Delta") + v = inj%vMod*inj%n(i) + CALL initVelDistDelta(inj%v(i)%obj, v) + + CASE DEFAULT + CALL criticalError("No velocity distribution type " // fvType // " defined", 'readVelDistr') + + END SELECT + + END DO + + END SUBROUTINE readVelDistr + + !Reads configuration for parallel run SUBROUTINE readParallel(config) USE moduleParallel USE moduleErrors diff --git a/src/modules/moduleInput.f95 b/src/modules/moduleInput.f95 deleted file mode 100644 index 25b00d3..0000000 --- a/src/modules/moduleInput.f95 +++ /dev/null @@ -1,559 +0,0 @@ -! moduleInput: Reads JSON configuration file -MODULE moduleInput - USE json_module - IMPLICIT NONE - - CONTAINS - !Main routine to read the input JSON file - SUBROUTINE readConfig(inputFile) - USE json_module - USE moduleErrors - USE moduleBoundary - USE moduleInject - IMPLICIT NONE - - CHARACTER(:), ALLOCATABLE, INTENT(in):: inputFile - TYPE(json_file):: config - - !Initialize the json file variable - CALL config%initialize() - - !Loads the config file - CALL verboseError('Loading input file...') - CALL config%load(filename = inputFile) - - !Reads reference parameters - CALL readReference(config) - - !Reads output parameters - CALL readOutput(config) - - !Read species - CALL readSpecies(config) - - !Read interactions between species - CALL readInteractions(config) - - !Read boundaries - CALL readBoundary(config) - - !Read Geometry - CALL verboseError('Reading Geometry...') - CALL readGeometry(config) - - !Reads case parameters - CALL verboseError('Reading Case Parameters...') - CALL readCase(config) - - !Read injection of particles - CALL verboseError('Reading Interactions between species...') - CALL readInject(config) - - !Read parallel parameters - CALL verboseError('Reading Parallel configuration...') - CALL readParallel(config) - - END SUBROUTINE readConfig - - !Reads the reference parameters - SUBROUTINE readReference(config) - USE moduleRefParam - USE moduleConstParam - USE moduleErrors - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - LOGICAL:: found, found_r - CHARACTER(:), ALLOCATABLE:: object - - object = 'reference' - !Mandatory parameters that define the case and computes derived parameters - CALL config%get(object // '.density', n_ref, found) - IF (.NOT. found) CALL criticalError('Reference density not found','readReference') - - CALL config%get(object // '.mass', m_ref, found) - IF (.NOT. found) CALL criticalError('Reference mass not found','readReference') - - 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_r) - - !Derived parameters - v_ref = DSQRT(kb*T_ref/m_ref) !reference velocity - !TODO: Make this solver dependent - IF (found_r) THEN - sigma_ref = PI*(r_ref+r_ref)**2 !reference cross section - L_ref = 1.D0/(sigma_ref*n_ref) !mean free path - - 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 - - !Reads the specific case parameters - SUBROUTINE readCase(config) - USE moduleRefParam - USE moduleErrors - USE moduleCaseParam - USE moduleSolver - USE moduleSpecies - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: object - REAL(8):: time !simulation time in [t] - CHARACTER(:), ALLOCATABLE:: pusherType, EMType, NAType - INTEGER:: nTau, nSolver - INTEGER:: i - CHARACTER(2):: iString - - object = 'case' - - !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(tau(1:nSpecies)) - DO i = 1, nTau - WRITE(iString, '(I2)') i - CALL config%get(object // '.tau(' // TRIM(iString) // ')', tau(i), found) - - END DO - IF (nTau < nSpecies) THEN - CALL warningError('Using minimum time step for some species') - tau(nTau+1:nSpecies) = MINVAL(tau(1:nTau)) - - END IF - tauMin = MINVAL(tau) - - !Gets the simulation time - CALL config%get(object // '.time', time, found) - IF (.NOT. found) CALL criticalError('Required parameter time not found','readCase') - !Convert simulation time to number of iterations - tmax = INT(time/tauMin) - - !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, tauMin, tau(i)) - - END DO - - !Gets the solver for the electromagnetic field - CALL config%get(object // '.EMSolver', EMType, found) - CALL solver%initEM(EMType) - SELECT CASE(EMType) - CASE("Electrostatic") - CALL readEMBoundary(config) - - END SELECT - - !Gest the non-analogue scheme - CALL config%get(object // '.NAScheme', NAType, found) - CALL solver%initNA(NAType) - - - !Makes tau non-dimensional - tau = tau / ti_ref - tauMin = tauMin / ti_ref - - END SUBROUTINE readCase - - !Reads configuration for the output files - SUBROUTINE readOutput(config) - USE moduleErrors - USE moduleOutput - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: object - CHARACTER(8) :: date_now='' - CHARACTER(10) :: time_now='' - - object = 'output' - CALL config%get(object // '.path', path, found) - CALL config%get(object // '.triggerOutput', triggerOutput, found) - IF (.NOT. found) THEN - triggerOutput = 100 - CALL warningError('Using default trigger for output file of 100 iterations') - - END IF - - !Creates output folder - !TODO: Add option for custon name output_folder - CALL DATE_AND_TIME(date_now, time_now) - folder = date_now(1:4) // '-' // date_now(5:6) // '-' // date_now(7:8) // '_' & - // time_now(1:2) // '.' // time_now(3:4) // '.' // time_now(5:6) - CALL SYSTEM('mkdir ' // path // folder ) - - CALL config%get(object // '.cpuTime', timeOutput, found) - CALL config%get(object // '.numColl', collOutput, found) - CALL config%get(object // '.EMField', emOutput, found) - - CALL config%get(object // '.triggerCPUTime', triggerCPUTime, found) - IF (.NOT. found) THEN - triggerCPUTime = triggerOutput - IF (timeOutput) CALL warningError('No triggerCPUTime found, using same vale as triggerOutput') - - END IF - - END SUBROUTINE readOutput - - !Reads information about the case species - SUBROUTINE readSpecies(config) - USE moduleSpecies - USE moduleErrors - USE moduleRefParam - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - CHARACTER(2):: iString - CHARACTER(:), ALLOCATABLE:: object - CHARACTER(:), ALLOCATABLE:: speciesType - REAL(8):: mass, charge - LOGICAL:: found - INTEGER:: i - - !Gets the number of species - CALL config%info('species', found, n_children = nSpecies) - !Zero species means critical error - IF (nSpecies == 0) CALL criticalError("No species found", "configRead") - - ALLOCATE(species(1:nSpecies)) - - !Reads information of individual species - DO i = 1, nSpecies - 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()) - - 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%sp = i - species(i)%obj%m = mass - - END DO - - !Set number of particles to 0 for init state - !TODO: In a future, this should include the particles from init states - nPartOld = 0 - - !Initialize the lock for the non-analogue (NA) list of particles - CALL OMP_INIT_LOCK(lockNAScheme) - - END SUBROUTINE readSpecies - - !Reads information about interactions between species - SUBROUTINE readInteractions(config) - USE moduleSpecies - USE moduleCollisions - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - CHARACTER(2):: iString, kString - CHARACTER(:), ALLOCATABLE:: object - 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, found) - - 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, found) - pt_i = speciesName2Index(species_i) - CALL config%get(object // '.species_j', species_j, found) - pt_j = speciesName2Index(species_j) - 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, found) - crossSecFilePath = pathCollisions // crossSecFile - CALL interactionMatrix(ij)%collisions(k)%obj%init(crossSecFilePath, species(pt_i)%obj%m, species(pt_j)%obj%m) - - END DO - - END DO - - END SUBROUTINE readInteractions - - !Reads boundary conditions for the mesh - SUBROUTINE readBoundary(config) - USE moduleBoundary - USE moduleErrors - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - CHARACTER(2):: istring - CHARACTER(:), ALLOCATABLE:: object - LOGICAL:: found - INTEGER:: i - - CALL config%info('boundary', found, n_children = nBoundary) - ALLOCATE(boundary(1:nBoundary)) - DO i = 1, nBoundary - WRITE(istring, '(i2)') i - object = 'boundary(' // trim(istring) // ')' - - ALLOCATE(boundaryGeneric:: boundary(i)%obj) - - CALL config%get(object // '.type', boundary(i)%obj%boundaryType, found) - CALL config%get(object // '.physicalSurface', boundary(i)%obj%physicalSurface, found) - boundary(i)%obj%id = i - - END DO - - END SUBROUTINE readBoundary - - !Read the geometry (mesh) for the case - SUBROUTINE readGeometry(config) - USE moduleMesh - USE moduleMeshCylRead, ONLY: meshCylGeneric - USE moduleMesh1DRead, ONLY: mesh1DGeneric - USE moduleErrors - USE moduleOutput - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: geometryType, meshFormat, meshFile - CHARACTER(:), ALLOCATABLE:: fullPath - - !Selects the type of geometry. - CALL config%get('geometry.type', geometryType, found) - SELECT CASE(geometryType) - CASE ("2DCyl") - !Creates a 2D cylindrical mesh - ALLOCATE(meshCylGeneric:: mesh) - - CASE ("1DCart") - !Creates a 1D cartesian mesh - ALLOCATE(mesh1DGeneric:: mesh) - - CASE DEFAULT - CALL criticalError("Geometry type " // geometryType // " not supported.", "readGeometry") - - END SELECT - - !Gets the type of mesh - CALL config%get('geometry.meshType', meshFormat, found) - CALL mesh%init(meshFormat) - !Reads the mesh - CALL config%get('geometry.meshFile', meshFile, found) - fullpath = path // meshFile - CALL mesh%readMesh(fullPath) - - 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 - USE moduleErrors - USE moduleInject - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - INTEGER:: i - CHARACTER(2):: istring - CHARACTER(:), ALLOCATABLE:: object - LOGICAL:: found - CHARACTER(:), ALLOCATABLE:: speciesName - CHARACTER(:), ALLOCATABLE:: name - REAL(8):: v - REAL(8), ALLOCATABLE:: T(:), normal(:) - REAL(8):: flow - CHARACTER(:), ALLOCATABLE:: units - INTEGER:: physicalSurface - INTEGER:: sp - - CALL config%info('inject', found, n_children = nInject) - ALLOCATE(inject(1:nInject)) - nPartInj = 0 - DO i = 1, nInject - WRITE(istring, '(i2)') i - object = 'inject(' // trim(istring) // ')' - - !Find species - CALL config%get(object // '.species', speciesName, found) - sp = speciesName2Index(speciesName) - - CALL config%get(object // '.name', name, found) - CALL config%get(object // '.v', v, found) - 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, units, sp, physicalSurface) - - END DO - - !Allocate array for injected particles - IF (nPartInj > 0) THEN - ALLOCATE(partInj(1:nPartInj)) - - END IF - - END SUBROUTINE readInject - - SUBROUTINE readParallel(config) - USE moduleParallel - USE moduleErrors - USE json_module - IMPLICIT NONE - - TYPE(json_file), INTENT(inout):: config - CHARACTER(:), ALLOCATABLE:: object - LOGICAL:: found - - !Reads OpenMP parameters - object = 'parallel.OpenMP' - - CALL config%get(object // '.nThreads', openMP%nThreads, found) - - IF (.NOT. found) THEN - openMP%nthreads = 8 - CALL warningError('No OpenMP configuration detected, using 8 threads as default.') - - END IF - - CALL initParallel() - - END SUBROUTINE readParallel - - -END MODULE moduleInput diff --git a/src/modules/moduleList.f95 b/src/modules/moduleList.f95 deleted file mode 100644 index bea9ba0..0000000 --- a/src/modules/moduleList.f95 +++ /dev/null @@ -1,91 +0,0 @@ -!Linked list of particles -MODULE moduleList - USE moduleSpecies - IMPLICIT NONE - - TYPE lNode - TYPE(particle), POINTER:: part => NULL() - TYPE(lNode), POINTER:: next => NULL() - - END TYPE lNode - - TYPE listNode - INTEGER:: amount = 0!TODO: Make private - TYPE(lNode),POINTER:: head => NULL() - TYPE(lNode),POINTER:: tail => NULL() - CONTAINS - PROCEDURE,PASS:: add => addToList - PROCEDURE,PASS:: convert2Array - PROCEDURE,PASS:: erase => eraseList - - END TYPE listNode - - TYPE(listNode):: partNAScheme !Particles comming from the nonAnalogue scheme - - TYPE pointerArray - TYPE(particle), POINTER:: part - - END TYPE - - CONTAINS - !Adds element to list - SUBROUTINE addToList(self,part) - USE moduleSpecies - CLASS(listNode), INTENT(inout):: self - TYPE(particle),INTENT(in), TARGET:: part - TYPE(lNode),POINTER:: temp - - ALLOCATE(temp) - temp%part => part - temp%next => NULL() - self%amount = self%amount + 1 - IF (.NOT. ASSOCIATED(self%head)) THEN - !First element - self%head => temp - self%tail => temp - ELSE - !Append element - self%tail%next => temp - self%tail => temp - END IF - - END SUBROUTINE addToList - - !converts list to array - FUNCTION convert2Array(self) RESULT(partArray) - IMPLICIT NONE - - CLASS(listNode), INTENT(in):: self - TYPE(pointerArray), ALLOCATABLE:: partArray(:) - TYPE(lNode), POINTER:: tempNode - INTEGER:: n - - ALLOCATE(partArray(1:self%amount)) - tempNode => self%head - DO n=1, self%amount - !Point element in array to element in list - partArray(n)%part => tempNode%part - !Go to next element - tempNode => tempNode%next - - END DO - - END FUNCTION convert2Array - - !Erase list - SUBROUTINE eraseList(self) - CLASS(listNode):: self - TYPE(lNode),POINTER:: current, next - - current => self%head - DO WHILE (ASSOCIATED(current)) - next => current%next - DEALLOCATE(current) - current => next - END DO - IF (ASSOCIATED(self%head)) NULLIFY(self%head) - IF (ASSOCIATED(self%tail)) NULLIFY(self%tail) - self%amount = 0 - END SUBROUTINE eraseList - -END MODULE moduleList diff --git a/src/modules/moduleOutput.f95 b/src/modules/moduleOutput.f95 deleted file mode 100644 index 33d9a94..0000000 --- a/src/modules/moduleOutput.f95 +++ /dev/null @@ -1,131 +0,0 @@ -!Contains information about output -MODULE moduleOutput - IMPLICIT NONE - !Output for each node - TYPE outputNode - REAL(8):: den, mom(1:3), tensorS(1:3,1:3) - - 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 - - END TYPE - - CHARACTER(:), ALLOCATABLE:: path - CHARACTER(:), ALLOCATABLE:: folder - INTEGER:: triggerOutput, counterOutput = 0 - INTEGER:: triggerCPUTime, counterCPUTime = 0 - LOGICAL:: timeOutput = .FALSE. - LOGICAL:: collOutput = .FALSE. - LOGICAL:: emOutput = .FALSE. - - CONTAINS - FUNCTION outerProduct(a,b) RESULT(s) - IMPLICIT NONE - - REAL(8), DIMENSION(1:3):: a, b - REAL(8):: s(1:3,1:3) - - s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3) - - END FUNCTION outerProduct - - FUNCTION tensorTrace(a) RESULT(t) - IMPLICIT NONE - - REAL(8), DIMENSION(1:3,1:3):: a - REAL(8):: t - - t = 0.D0 - t = a(1,1)+a(2,2)+a(3,3) - - END FUNCTION tensorTrace - - SUBROUTINE calculateOutput(rawValues, formatValues, nodeVol, speciesIn) - USE moduleConstParam - USE moduleRefParam - USE moduleSpecies - IMPLICIT NONE - - TYPE(outputNode), INTENT(in):: rawValues - TYPE(outputFormat), INTENT(out):: formatValues - REAL(8), INTENT(in):: nodeVol - CLASS(speciesGeneric), INTENT(in):: speciesIn - REAL(8), DIMENSION(1:3,1:3):: tensorTemp - REAL(8), DIMENSION(1:3):: tempVel - REAL(8):: tempVol - - !Resets the node outputs - formatValues%density = 0.D0 - formatValues%velocity = 0.D0 - formatValues%pressure = 0.D0 - formatValues%temperature = 0.D0 - tempVol = 1.D0/(nodeVol*Vol_ref) - IF (rawValues%den > 0.D0) THEN - tempVel = rawValues%mom(:)/rawValues%den - IF ((tempVel(1) - 1.D0) .EQ. tempVel(1)) THEN - PRINT *, rawValues%mom - END IF - tensorTemp = (rawValues%tensorS(:,:) - rawValues%den*outerProduct(tempVel,tempVel)) - formatValues%density = rawValues%den*tempVol - formatValues%velocity(:) = tempVel - IF (tensorTrace(tensorTemp) > 0.D0) THEN - formatValues%pressure = speciesIn%m*tensorTrace(tensorTemp)*tempVol/3.D0 - formatValues%temperature = formatValues%pressure/(formatValues%density*kb) - - END IF - END IF - - formatValues%velocity = formatValues%velocity*v_ref - formatValues%pressure = formatValues%pressure*m_ref*v_ref**2 - formatValues%temperature = formatValues%temperature*m_ref*v_ref**2 - - END SUBROUTINE calculateOutput - - SUBROUTINE printTime(t, first) - USE moduleSpecies - USE moduleCompTime - IMPLICIT NONE - - INTEGER, INTENT(in):: t - LOGICAL, INTENT(in), OPTIONAL:: first - CHARACTER(:), ALLOCATABLE:: fileName - - fileName = 'cpuTime.dat' - - IF (timeOutput) THEN - IF (PRESENT(first)) THEN - IF (first) THEN - OPEN(20, file = path // folder // '/' // fileName, action = 'write') - WRITE(20, "(A1, 8X, A1, 9X, A1, 5(A20))") "#","t","n","total","push","reset","collision","weighting" - WRITE(*, "(6X,A15,A)") "Creating file: ", fileName - - ELSE - - END IF - OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - - ELSE - OPEN(20, file = path // folder // '/' // fileName, position = 'append', action = 'write') - - END IF - - WRITE (20, "(I10, I10, 5(ES20.6E3))") t, nPartOld, tStep, tPush, tReset, tColl, tWeight - - CLOSE(20) - - END IF - - END SUBROUTINE printTime - -END MODULE moduleOutput - diff --git a/src/modules/moduleParallel.f95 b/src/modules/moduleParallel.f95 deleted file mode 100644 index 564724c..0000000 --- a/src/modules/moduleParallel.f95 +++ /dev/null @@ -1,20 +0,0 @@ -MODULE moduleParallel - IMPLICIT NONE - - TYPE openMP_type - INTEGER:: nThreads - - END TYPE - - TYPE(openMP_type):: openMP - - CONTAINS - SUBROUTINE initParallel() - IMPLICIT NONE - - !Starts threads for OpenMP parallelization - CALL OMP_SET_NUM_THREADS(openMP%nThreads) - - END SUBROUTINE initParallel - -END MODULE moduleParallel diff --git a/src/modules/moduleRefParam.f95 b/src/modules/moduleRefParam.f95 deleted file mode 100644 index 9d2b8d9..0000000 --- a/src/modules/moduleRefParam.f95 +++ /dev/null @@ -1,9 +0,0 @@ -!Reference parameters -MODULE moduleRefParam - !Parameters that define the problem (inputs) - 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, EF_ref, Volt_ref - -END MODULE moduleRefParam - diff --git a/src/modules/moduleSolver.f95 b/src/modules/moduleSolver.f95 deleted file mode 100644 index 60c1bd3..0000000 --- a/src/modules/moduleSolver.f95 +++ /dev/null @@ -1,584 +0,0 @@ -MODULE moduleSolver - - !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, 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 - - 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 - - !Apply nonAnalogue scheme to a particle - SUBROUTINE nonAnalogue_interface(part, volOld, volNew) - USE moduleSpecies - USE moduleMesh - IMPLICIT NONE - - TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER, INTENT(in):: volOld - CLASS(meshVol), POINTER, INTENT(inout):: volNew - - END SUBROUTINE nonAnalogue_interface - - END INTERFACE - - 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('1DCartCharged') - self%pushParticle => push1DCharged - - CASE DEFAULT - CALL criticalError('Solver ' // pusherType // ' not found','readCase') - - END SELECT - - self%pushSpecies = .FALSE. - self%every = INT(tauSp/tau) - - 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 - - END SELECT - - END SUBROUTINE initEM - - !Initialize the non-analogue scheme - SUBROUTINE initNA(self, NAType) - USE moduleMesh - IMPLICIT NONE - - CLASS(solverGeneric), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE:: NAType - - SELECT CASE(NAType) - CASE ('Volume') - self%nonAnalogue => volumeNAScheme - - 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 - !Checks if the species sp is update this iteration - IF (solver%pusher(sp)%pushSpecies) THEN - !Push particle - CALL solver%pusher(sp)%pushParticle(partOld(n)) - !Find cell in wich particle reside - CALL solver%updateParticleCell(partOld(n)) - - END IF - - END DO - !$OMP END DO - - END SUBROUTINE doPushes - - !Push one particle. Boris pusher for 2D Cyl Neutral particle - PURE SUBROUTINE pushCylNeutral(part) - USE moduleSpecies - IMPLICIT NONE - - TYPE(particle), INTENT(inout):: part - TYPE(particle):: part_temp - REAL(8):: tauSp - REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha - REAL(8):: v_p_oh_star(2:3) - - part_temp = part - !Time step for the species - tauSp = tau(part_temp%sp) - !z - part_temp%v(1) = part%v(1) - part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp - !r,theta - v_p_oh_star(2) = part%v(2) - x_new = part%r(2) + v_p_oh_star(2)*tauSp - v_p_oh_star(3) = part%v(3) - y_new = v_p_oh_star(3)*tauSp - 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 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):: tauSp - REAL(8):: qmEFt(1:3)!charge*tauSp*EF/mass - - part_temp = part - !Time step for the species - tauSp = tau(part_temp%sp) - !Get electric field at particle position - qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp - !z - part_temp%v(1) = part%v(1) + qmEFt(1) - part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp - !r,theta - v_p_oh_star(2) = part%v(2) + qmEFt(2) - x_new = part%r(2) + v_p_oh_star(2)*tauSp - v_p_oh_star(3) = part%v(3) + qmEFt(3) - y_new = v_p_oh_star(3)*tauSp - 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 - - !Push charged particles in 1D cartesian coordinates - PURE SUBROUTINE push1DCharged(part) - USE moduleSPecies - USE moduleEM - IMPLICIT NONE - - TYPE(particle), INTENT(inout):: part - TYPE(particle):: part_temp - REAL(8):: tauSp - REAL(8):: qmEFt(1:3) - - part_temp = part - !Time step for particle species - tauSp = tau(part_temp%sp) - !Get the electric field at particle position - qmEFt = part_temp%qm*gatherElecField(part_temp)*tauSp - - !x - part_temp%v(1) = part%v(1) + qmEFt(1) - part_temp%r(1) = part%r(1) + part_temp%v(1)*tauSp - - part_temp%n_in = .FALSE. - - part = part_temp - - END SUBROUTINE push1DCharged - - !Do the collisions in all the cells - SUBROUTINE doCollisions() - USE moduleMesh - IMPLICIT NONE - - INTEGER:: e - - !$OMP DO SCHEDULE(DYNAMIC) - DO e=1, mesh%numVols - CALL mesh%vols(e)%obj%collision() - END DO - !$OMP END DO - - END SUBROUTINE doCollisions - - SUBROUTINE doReset() - USE moduleSpecies - USE moduleList - IMPLICIT NONE - - INTEGER:: nn, n - INTEGER, SAVE:: nPartNew - INTEGER, SAVE:: nInjIn, nOldIn, nNAScheme - TYPE(particle), ALLOCATABLE, SAVE:: partTemp(:) - TYPE(lNode), POINTER:: partCurr, partNext - - !$OMP SECTIONS - !$OMP SECTION - nInjIn = 0 - IF (ALLOCATED(partInj)) THEN - nInjIn = COUNT(partInj%n_in) - - END IF - !$OMP SECTION - nOldIn = 0 - IF (ALLOCATED(partOld)) THEN - nOldIn = COUNT(partOld%n_in) - - END IF - !$OMP SECTION - nNAScheme = partNAScheme%amount - !$OMP END SECTIONS - - !$OMP BARRIER - - !$OMP SINGLE - CALL MOVE_ALLOC(partOld, partTemp) - nPartNew = nInjIn + nOldIn + nNAScheme - ALLOCATE(partOld(1:nPartNew)) - !$OMP END SINGLE - - !$OMP SECTIONS - !$OMP SECTION - nn = 0 - DO n = 1, nPartInj - IF (partInj(n)%n_in) THEN - nn = nn + 1 - partOld(nn) = partInj(n) - - END IF - - END DO - - !$OMP SECTION - nn = nInjIn - DO n = 1, nPartOld - IF (partTemp(n)%n_in) THEN - nn = nn + 1 - partOld(nn) = partTemp(n) - - END IF - - END DO - !$OMP SECTION - nn = nInjIn + nOldIn - partCurr => partNAScheme%head - DO n = 1, nNAScheme - partNext => partCurr%next - partOld(nn+n) = partCurr%part - DEALLOCATE(partCurr) - partCurr => partNext - - END DO - IF (ASSOCIATED(partNAScheme%head)) NULLIFY(partNAScheme%head) - IF (ASSOCIATED(partNAScheme%tail)) NULLIFY(partNAScheme%tail) - partNAScheme%amount = 0 - - !$OMP END SECTIONS - - !$OMP SINGLE - nPartOld = nPartNew - !$OMP END SINGLE - - END SUBROUTINE doReset - - !Scatter particles in the grid - SUBROUTINE doScatter - USE moduleSpecies - USE moduleMesh - IMPLICIT NONE - - INTEGER:: n - - !Loops over the particles to scatter them - !$OMP DO - DO n=1, nPartOld - CALL mesh%vols(partOld(n)%vol)%obj%scatter(partOld(n)) - - END DO - !$OMP END DO - - END SUBROUTINE doScatter - - SUBROUTINE doEMField() - IMPLICIT NONE - - IF (ASSOCIATED(solver%solveEM)) THEN - CALL solver%solveEM() - - END IF - - 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 - - !Modify particle weight as a function of cell volume and splits particle - SUBROUTINE volumeNAScheme(part, volOld, volNew) - USE moduleSpecies - USE moduleMesh - IMPLICIT NONE - - TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER, INTENT(in):: volOld - CLASS(meshVol), POINTER, INTENT(inout):: volNew - REAL(8):: fractionVolume, fractionWeight - INTEGER:: nSplit - - !If particle has change cell, call nonAnalogue scheme - IF (volOld%n /= volNew%n) THEN - fractionVolume = volOld%volume/volNew%volume - - part%weight = part%weight * fractionVolume - - fractionWeight = part%weight / species(part%sp)%obj%weight - - IF (fractionWeight >= 2.D0) THEN - nSplit = FLOOR(fractionWeight) - CALL splitParticle(part, nSplit, volNew) - - ELSEIF (part%weight < 1.D0) THEN - !Particle has lost statistical meaning and will be terminated - part%n_in = .FALSE. - - END IF - - END IF - - END SUBROUTINE volumeNAScheme - - !Subroutine to split the particle 'part' into a number 'nSplit' of particles. - !'nSplit-1' particles are added to the partNAScheme list - SUBROUTINE splitParticle(part, nSplit, vol) - USE moduleSpecies - USE moduleList - USE moduleMesh - USE OMP_LIB - IMPLICIT NONE - - TYPE(particle), INTENT(inout):: part - INTEGER, INTENT(in):: nSplit - CLASS(meshVol), INTENT(inout):: vol - REAL(8):: newWeight - TYPE(particle), POINTER:: newPart - INTEGER:: p - - newWeight = part%weight / nSplit - - !Assign new weight to original particle - part%weight = newWeight - - !Add new particles to list of NA particles - DO p = 2, nSplit - !Allocate the pointer for the new particles - ALLOCATE(newPart) - !Copy data from original particle - newPart = part - CALL OMP_SET_LOCK(lockNAScheme) - CALL partNAScheme%add(newPart) - CALL OMP_UNSET_LOCK(lockNASCheme) - !Add particle to cell list - CALL OMP_SET_lock(vol%lock) - CALL vol%listPart_in%add(newPart) - CALL OMP_UNSET_lock(vol%lock) - - END DO - - END SUBROUTINE splitParticle - - SUBROUTINE updateParticleCell(self, part) - USE moduleSpecies - USE moduleMesh - IMPLICIT NONE - - CLASS(solverGeneric), INTENT(in):: self - TYPE(particle), INTENT(inout):: part - CLASS(meshVol), POINTER:: volOld, volNew - - volOld => mesh%vols(part%vol)%obj - CALL volOld%findCell(part) - volNew => mesh%vols(part%vol)%obj - !Call the NA shcme - IF (ASSOCIATED(self%nonAnalogue)) THEN - CALL self%nonAnalogue(part, volOld, volNew) - - 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 - USE moduleSpecies - USE moduleCompTime - IMPLICIT NONE - - INTEGER, INTENT(in):: t - - counterOutput = counterOutput + 1 - IF (counterOutput >= triggerOutput .OR. & - t == tmax .OR. t == 0) THEN - - !Resets output counter - counterOutput=0 - - CALL mesh%printOutput(t) - CALL mesh%printColl(t) - 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" - - ELSE - WRITE(*, "(5X,A21,F8.1,A2)") "iteration time: ", 1.D3*tStep, "ms" - - END IF - - IF (nPartOld > 0) THEN - WRITE(*, "(5X,A21,F8.1,A2)") "avg t/particle: ", 1.D9*tStep/DBLE(nPartOld), "ns" - - END IF - WRITE(*,*) - - END IF - - counterCPUTime = counterCPUTime + 1 - IF (counterCPUTime >= triggerCPUTime .OR. & - t == tmax .OR. t == 0) THEN - - !Reset CPU Time counter - counterCPUTime = 0 - - CALL printTime(t, t == 0) - - END IF - - END SUBROUTINE doOutput - - -END MODULE moduleSolver - diff --git a/src/modules/moduleSpecies.f95 b/src/modules/moduleSpecies.f95 deleted file mode 100644 index da5128c..0000000 --- a/src/modules/moduleSpecies.f95 +++ /dev/null @@ -1,71 +0,0 @@ -!Contains the information about species (particles) -MODULE moduleSpecies - USE moduleCaseParam - USE OMP_LIB - IMPLICIT NONE - - TYPE, ABSTRACT:: speciesGeneric - CHARACTER(:), ALLOCATABLE:: name - REAL(8):: m=0.D0, weight=0.D0 - INTEGER:: sp=0 - END TYPE speciesGeneric - - TYPE, EXTENDS(speciesGeneric):: speciesNeutral - - END TYPE speciesNeutral - - TYPE, EXTENDS(speciesGeneric):: speciesCharged - REAL(8):: q=0.D0, qm=0.D0 - - END TYPE speciesCharged - - TYPE:: speciesCont - CLASS(speciesGeneric), ALLOCATABLE:: obj - - END TYPE - - INTEGER:: nSpecies - TYPE(speciesCont), ALLOCATABLE:: species(:) - - TYPE particle - REAL(8):: r(1:3) !Position - REAL(8):: v(1:3) !Velocity - 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 - REAL(8):: qm = 0.D0 !charge over mass - - END TYPE particle - - !Number of old particles - INTEGER:: nPartOld - INTEGER:: nPartInj - !Arrays that contain the particles - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partOld !array of particles from previous iteration - TYPE(particle), ALLOCATABLE, DIMENSION(:), TARGET:: partInj !array of inject particles - INTEGER(KIND=OMP_LOCK_KIND):: lockNAScheme !Lock for the NA list of particles - - CONTAINS - FUNCTION speciesName2Index(speciesName) RESULT(sp) - USE moduleErrors - IMPLICIT NONE - - CHARACTER(:), ALLOCATABLE:: speciesName - INTEGER:: sp - INTEGER:: n - - sp = 0 - DO n = 1, nSpecies - IF (speciesName == species(n)%obj%name) THEN - sp = species(n)%obj%sp - EXIT - END IF - END DO - !If no species is found, call a critical error - IF (sp == 0) CALL criticalError('Species ' // speciesName // ' not found.', 'speciesName2Index') - - END FUNCTION speciesName2Index - -END MODULE moduleSpecies diff --git a/src/modules/moduleTable.f95 b/src/modules/moduleTable.f95 deleted file mode 100644 index 1ea89e5..0000000 --- a/src/modules/moduleTable.f95 +++ /dev/null @@ -1,121 +0,0 @@ -MODULE moduleTable - - TYPE:: table1D - REAL(8):: xMin, xMax - REAL(8):: fMin, fMax - REAL(8), ALLOCATABLE, DIMENSION(:):: x, f, k - CONTAINS - PROCEDURE, PASS:: init => initTable1D - PROCEDURE, PASS:: get => getValueTable1D - PROCEDURE, PASS:: convert => convertUnits - - END TYPE table1D - - CONTAINS - SUBROUTINE initTable1D(self, tableFile) - USE moduleErrors - IMPLICIT NONE - - CLASS(table1D), INTENT(inout):: self - CHARACTER(:), ALLOCATABLE, INTENT(IN):: tableFile - CHARACTER(100):: dummy - INTEGER:: amount - INTEGER:: i - INTEGER:: stat - INTEGER:: id = 20 - - OPEN(id, file = tableFile) - amount = 0 - DO - READ(id, '(A)', iostat = stat) dummy - !Skip comment - IF (INDEX(dummy,'#') /= 0) CYCLE - !If EOF or error, exit file - IF (stat /= 0) EXIT - !Add row - amount = amount + 1 - - END DO - - IF (amount == 0) CALL criticalError('Table ' // tableFile // ' is empty', 'initTable1D') - IF (amount == 1) CALL criticalError('Table ' // tableFile // ' has only one row', 'initTable1D') - - !Go bback to initial point - REWIND(id) - - !Allocate table arrays - ALLOCATE(self%x(1:amount), self%f(1:amount), self%k(1:amount)) - self%x = 0.D0 - self%f = 0.D0 - self%k = 0.D0 - - i = 0 - DO - READ(id, '(A)', iostat = stat) dummy - IF (INDEX(dummy,'#') /= 0) CYCLE - IF (stat /= 0) EXIT - !Add data - !TODO: substitute with extracting information from dummy - BACKSPACE(id) - i = i + 1 - READ(id, *) self%x(i), self%f(i) - - END DO - - CLOSE(10) - - self%xMin = self%x(1) - self%xMax = self%x(amount) - self%fMin = self%f(1) - self%fMax = self%f(amount) - - DO i = 1, amount - 1 - self%k(i) = (self%f(i+1) - self%f(i))/(self%x(i+1) - self%x(i)) - - END DO - - END SUBROUTINE initTable1D - - FUNCTION getValueTable1D(self, x) RESULT(f) - IMPLICIT NONE - - CLASS(table1D), INTENT(in):: self - REAL(8):: x - REAL(8):: f - REAL(8):: deltaX - INTEGER:: i - - IF (x <= self%xMin) THEN - f = self%fMin - - ELSEIF (x >= self%xMax) THEN - f = self%fMax - - ELSE - i = MINLOC(x - self%x, 1) - deltaX = x - self%x(i) - IF (deltaX < 0 ) THEN - i = i - 1 - deltaX = x - self%x(i) - - END IF - - f = self%f(i) + self%k(i)*deltaX - - END IF - - END FUNCTION getValueTable1D - - SUBROUTINE convertUnits(self, data_x, data_f) - IMPLICIT NONE - - CLASS(table1D), INTENT(inout):: self - REAL(8):: data_x, data_f - - self%x = self%x * data_x - self%f = self%f * data_f - self%k = self%k * data_f / data_x - - END SUBROUTINE convertUnits - -END MODULE moduleTable