JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Structure for 3D Cartesian Grid created.

Browse source 
Unification of boundary conditions into one file.

Some changes to input file for reference cases. This should have been
done in another branch but I wanto to commit to save progress and I
don't want to deal with tswitching branches right now, I'm very busy
watching Futurama.
This commit is contained in:
Jorge Gonzalez 2021-02-27 16:24:44 +01:00
commit ac2965621a
29 changed files with 1549 additions and 40455 deletions
Download patch file Download diff file Expand all files Collapse all files

1
.gitignore vendored
View file

@ -4,3 +4,4 @@ obj/
doc/user_manual/ doc/user_manual/
doc/coding_style/ doc/coding_style/
json-fortran-8.2.0/ json-fortran-8.2.0/
runs/

1
makefile
View file

@ -39,5 +39,6 @@ src.o:
clean: clean:
rm -f $(OUTPUT) rm -f $(OUTPUT)
rm -f $(MODDIR)/*.mod rm -f $(MODDIR)/*.mod
rm -f $(MODDIR)/*.smod
rm -f $(OBJDIR)/*.o rm -f $(OBJDIR)/*.o

19
runs/1D_Cathode/inputRad.json
View file

@ -19,13 +19,16 @@
"meshFile": "mesh.msh" "meshFile": "mesh.msh"
}, },
"species": [ "species": [
{"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1} {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1},
{"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e1}
], ],
"boundary": [ "boundary": [
{"name": "Cathode", "physicalSurface": 1, "bTypes": [ {"name": "Cathode", "physicalSurface": 1, "bTypes": [
{"type": "absorption"},
{"type": "absorption"} {"type": "absorption"}
]}, ]},
{"name": "Infinite", "physicalSurface": 2, "bTypes": [ {"name": "Infinite", "physicalSurface": 2, "bTypes": [
{"type": "transparent"},
{"type": "transparent"} {"type": "transparent"}
]} ]}
], ],
@ -35,17 +38,21 @@
], ],
"inject": [ "inject": [
{"name": "Cathode Electron", "species": "Electron", "flow": 1.0e-3, "units": "A", "v": 27500.0, "T": [2500.0, 2500.0, 2500.0], {"name": "Cathode Electron", "species": "Electron", "flow": 1.0e-3, "units": "A", "v": 27500.0, "T": [2500.0, 2500.0, 2500.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1} "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1},
{"name": "Plasma Inf Ar+", "species": "Argon+", "flow": 1.0e-6, "units": "A", "v": 500.0, "T": [500.0, 500.0, 500.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2},
{"name": "Plasma Inf e", "species": "Electron", "flow": 1.0e-6, "units": "A", "v": 500.0, "T": [500.0, 500.0, 500.0],
"velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [-1, 0, 0], "physicalSurface": 2}
], ],
"case": { "case": {
"tau": [1.0e-11], "tau": [1.0e-11, 1.0e-11],
"time": 1.0e-7, "time": 1.0e-5,
"pusher": ["1DRadCharged"], "pusher": ["1DRadCharged", "1DRadCharged"],
"EMSolver": "Electrostatic" "EMSolver": "Electrostatic"
}, },
"parallel": { "parallel": {
"OpenMP":{ "OpenMP":{
"nThreads": 1 "nThreads": 24
} }
} }
} }

16
runs/1D_Cathode/mesh.geo Normal file
View file

@ -0,0 +1,16 @@
Lcell = 0.00005;
x0 = 0.001;
xf = x0 + 100.0*Lcell;
Point(1) = {x0, 0, 0, 1};
Point(2) = {xf, 0, 0, 1};
Line(1) = {1, 2};
Physical Point(1) = {1};
Physical Point(2) = {2};
Physical Line(1) = {1};
Transfinite Line {1} = (xf-x0)/Lcell + 1 Using Progression 1;

39806
runs/1D_Cathode/mesh.msh
View file

File diff suppressed because it is too large Load diff

3
runs/ALPHIE_Grid/inputDiffTau.json
View file

@ -45,7 +45,8 @@
], ],
"boundaryEM": [ "boundaryEM": [
{"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4}, {"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
{"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5} {"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
{"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}
], ],
"inject": [ "inject": [
{"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0], {"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0],

9
runs/ALPHIE_Grid/inputSameTau.json
View file

@ -14,8 +14,8 @@
"meshFile": "mesh.msh" "meshFile": "mesh.msh"
}, },
"species": [ "species": [
{"name": "Argon+", "type": "charged", "mass": 6.633e-26, "charge": 1.0, "weight": 1.0e2}, {"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.0e2} {"name": "Electron", "type": "charged", "mass": 9.109e-31, "charge":-1.0, "weight": 1.0e1}
], ],
"boundary": [ "boundary": [
{"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [ {"name": "Ionization Chanber", "physicalSurface": 1, "bTypes": [
@ -45,7 +45,8 @@
], ],
"boundaryEM": [ "boundaryEM": [
{"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4}, {"name": "Extraction Grid", "type": "dirichlet", "potential": -150.0, "physicalSurface": 4},
{"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5} {"name": "Acceleration Grid", "type": "dirichlet", "potential": -600.0, "physicalSurface": 5},
{"name": "Ionization Chamber", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1}
], ],
"inject": [ "inject": [
{"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0], {"name": "Ionization Argon+", "species": "Argon+", "flow": 27.0e-6, "units": "A", "v": 322.0, "T": [ 500.0, 500.0, 500.0],
@ -63,7 +64,7 @@
}, },
"case": { "case": {
"tau": [1.0e-11, 1.0e-11], "tau": [1.0e-11, 1.0e-11],
"time": 2.0e-6, "time": 1.0e-6,
"pusher": ["2DCylCharged", "2DCylCharged"], "pusher": ["2DCylCharged", "2DCylCharged"],
"EMSolver": "Electrostatic" "EMSolver": "Electrostatic"
}, },

37
runs/ALPHIE_Grid/mesh.geo
View file

@ -1,5 +1,3 @@
Lz = 0.0100;
Lr = 0.0006;
zg1 = 0.0025; zg1 = 0.0025;
tg1 = 0.0004; tg1 = 0.0004;
rg1 = 0.0005; rg1 = 0.0005;
@ -7,6 +5,9 @@ dg = 0.0025;
zg2 = zg1+tg1+dg; zg2 = zg1+tg1+dg;
tg2 = tg1; tg2 = tg1;
rg2 = rg1; rg2 = rg1;
zEnd = 0.0042;
Lz = zg2 + tg2 + zEnd;
Lr = rg1 + 0.0001;
Lcell = 0.0001; Lcell = 0.0001;
@ -96,19 +97,19 @@ Transfinite Line {18, 19, 20, 21, 22, 6} = rg1/Lcell + 1 Using Progression 1;
Transfinite Line {17, 15, 13, 11, 9, 7} = (Lr-rg1)/Lcell + 1 Using Progression 1; Transfinite Line {17, 15, 13, 11, 9, 7} = (Lr-rg1)/Lcell + 1 Using Progression 1;
#Transfinite Surface{1}; Transfinite Surface{1};
#Recombine Surface {1}; Recombine Surface {1};
#Transfinite Surface{2}; Transfinite Surface{2};
#Recombine Surface {2}; Recombine Surface {2};
#Transfinite Surface{3}; Transfinite Surface{3};
#Recombine Surface {3}; Recombine Surface {3};
#Transfinite Surface{4}; Transfinite Surface{4};
#Recombine Surface {4}; Recombine Surface {4};
#Transfinite Surface{5}; Transfinite Surface{5};
#Recombine Surface {5}; Recombine Surface {5};
#Transfinite Surface{6}; Transfinite Surface{6};
#Recombine Surface {6}; Recombine Surface {6};
#Transfinite Surface{7}; Transfinite Surface{7};
#Recombine Surface {7}; Recombine Surface {7};
#Transfinite Surface{8}; Transfinite Surface{8};
#Recombine Surface {8}; Recombine Surface {8};

10
src/makefile
View file

@ -1,13 +1,13 @@
OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolver.o \
$(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \ $(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \
$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \ $(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \ $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \
$(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o $(OBJDIR)/moduleMesh2DCylBoundary.o \ $(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o \
$(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o $(OBJDIR)/moduleMesh2DCartBoundary.o \ $(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o \
$(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o $(OBJDIR)/moduleMesh1DCartBoundary.o \ $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o \
$(OBJDIR)/moduleMesh1DRad.o $(OBJDIR)/moduleMesh1DRadRead.o $(OBJDIR)/moduleMesh1DRadBoundary.o $(OBJDIR)/moduleMesh1DRad.o $(OBJDIR)/moduleMesh1DRadRead.o
all: $(OUTPUT) all: $(OUTPUT)

7
src/modules/mesh/1DCart/makefile
View file

@ -1,11 +1,8 @@
all: moduleMesh1DCart.o moduleMesh1DCartBoundary.o moduleMesh1DCartRead.o all: moduleMesh1DCart.o moduleMesh1DCartRead.o
moduleMesh1DCart.o: moduleMesh1DCart.f90 moduleMesh1DCart.o: moduleMesh1DCart.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh1DCartBoundary.o: moduleMesh1DCart.o moduleMesh1DCartBoundary.f90 moduleMesh1DCartRead.o: moduleMesh1DCart.o moduleMesh1DCartRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh1DCartRead.o: moduleMesh1DCart.o moduleMesh1DCartBoundary.o moduleMesh1DCartRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

53
src/modules/mesh/1DCart/moduleMesh1DCart.f90
View file

@ -4,6 +4,7 @@
! z == unused ! z == unused
MODULE moduleMesh1DCart MODULE moduleMesh1DCart
USE moduleMesh USE moduleMesh
USE moduleMeshBoundary
IMPLICIT NONE IMPLICIT NONE
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart
@ -23,50 +24,11 @@ MODULE moduleMesh1DCart
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initEdge1DCart PROCEDURE, PASS:: init => initEdge1DCart
PROCEDURE, PASS:: getNodes => getNodes1DCart PROCEDURE, PASS:: getNodes => getNodes1DCart
PROCEDURE, PASS:: intersection => intersection1DCart
PROCEDURE, PASS:: randPos => randPosEdge PROCEDURE, PASS:: randPos => randPosEdge
END TYPE meshEdge1DCart END TYPE meshEdge1DCart
!Boundary functions defined in the submodule Boundary
INTERFACE
MODULE SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE reflection
MODULE SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE absorption
MODULE SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE transparent
MODULE SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE wallTemperature
END INTERFACE
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DCart
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ1DCart PROCEDURE, PASS:: detJac => detJ1DCart
@ -226,6 +188,17 @@ MODULE moduleMesh1DCart
END FUNCTION getNodes1DCart END FUNCTION getNodes1DCart
PURE FUNCTION intersection1DCart(self, r0, v0) RESULT(r)
IMPLICIT NONE
CLASS(meshEdge1DCart), INTENT(in):: self
REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
REAL(8), DIMENSION(1:3):: r
r = (/ self%x, 0.D0, 0.D0 /)
END FUNCTION intersection1DCart
!Calculates a 'random' position in edge !Calculates a 'random' position in edge
FUNCTION randPosEdge(self) RESULT(r) FUNCTION randPosEdge(self) RESULT(r)
CLASS(meshEdge1DCart), INTENT(in):: self CLASS(meshEdge1DCart), INTENT(in):: self

91
src/modules/mesh/1DCart/moduleMesh1DCartBoundary.f90
View file

@ -1,91 +0,0 @@
SUBMODULE (moduleMesh1DCart) moduleMesh1DCartBoundary
USE moduleMesh1DCart
CONTAINS
SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
SELECT TYPE(edge)
TYPE IS(meshEdge1DCart)
part%v(1) = -part%v(1)
part%r(1) = 2.D0*edge%x - part%r(1)
END SELECT
END SUBROUTINE reflection
SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: rEdge(1) !Position of particle in the edge
REAL(8):: d !Distance from particle to edge
SELECT TYPE(edge)
TYPE IS(meshEdge1DCart)
rEdge(1) = edge%x
d = DABS(part%r(1) - rEdge(1))
IF (d > 0.D0) THEN
part%weight = part%weight / d
part%r(1) = rEdge(1)
END IF
IF (ASSOCIATED(edge%e1)) THEN
CALL edge%e1%scatter(part)
ELSE
CALL edge%e2%scatter(part)
END IF
END SELECT
part%n_in = .FALSE.
END SUBROUTINE absorption
SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
part%n_in = .FALSE.
END SUBROUTINE transparent
SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
USE moduleBoundary
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
!Modifies particle velocity according to wall temperature
SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
TYPE IS(boundaryWallTemperature)
part%v(1) = part%v(1) + bound%vTh*randomMaxwellian()
END SELECT
SELECT TYPE(edge)
TYPE IS(meshEdge1DCart)
part%v(1) = -part%v(1)
part%r(1) = 2.D0*edge%x - part%r(1)
END SELECT
END SUBROUTINE wallTemperature
END SUBMODULE moduleMesh1DCartBoundary

7
src/modules/mesh/1DRad/makefile
View file

@ -1,11 +1,8 @@
all: moduleMesh1DRad.o moduleMesh1DRadBoundary.o moduleMesh1DRadRead.o all: moduleMesh1DRad.o moduleMesh1DRadRead.o
moduleMesh1DRad.o: moduleMesh1DRad.f90 moduleMesh1DRad.o: moduleMesh1DRad.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh1DRadBoundary.o: moduleMesh1DRad.o moduleMesh1DRadBoundary.f90 moduleMesh1DRadRead.o: moduleMesh1DRad.o moduleMesh1DRadRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh1DRadRead.o: moduleMesh1DRad.o moduleMesh1DRadBoundary.o moduleMesh1DRadRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

53
src/modules/mesh/1DRad/moduleMesh1DRad.f90
View file

@ -4,6 +4,7 @@
! z == unused ! z == unused
MODULE moduleMesh1DRad MODULE moduleMesh1DRad
USE moduleMesh USE moduleMesh
USE moduleMeshBoundary
IMPLICIT NONE IMPLICIT NONE
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad
@ -23,50 +24,11 @@ MODULE moduleMesh1DRad
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initEdge1DRad PROCEDURE, PASS:: init => initEdge1DRad
PROCEDURE, PASS:: getNodes => getNodes1DRad PROCEDURE, PASS:: getNodes => getNodes1DRad
PROCEDURE, PASS:: intersection => intersection1DRad
PROCEDURE, PASS:: randPos => randPos1DRad PROCEDURE, PASS:: randPos => randPos1DRad
END TYPE meshEdge1DRad END TYPE meshEdge1DRad
!Boundary functions defined in the submodule Boundary
INTERFACE
MODULE SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE reflection
MODULE SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE absorption
MODULE SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE transparent
MODULE SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE wallTemperature
END INTERFACE
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol1DRad
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ1DRad PROCEDURE, PASS:: detJac => detJ1DRad
@ -227,6 +189,17 @@ MODULE moduleMesh1DRad
END FUNCTION getNodes1DRad END FUNCTION getNodes1DRad
PURE FUNCTION intersection1DRad(self, r0, v0) RESULT(r)
IMPLICIT NONE
CLASS(meshEdge1DRad), INTENT(in):: self
REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
REAL(8), DIMENSION(1:3):: r
r = (/ self%r, 0.D0, 0.D0 /)
END FUNCTION intersection1DRad
!Calculates a 'random' position in edge !Calculates a 'random' position in edge
FUNCTION randPos1DRad(self) RESULT(r) FUNCTION randPos1DRad(self) RESULT(r)
CLASS(meshEdge1DRad), INTENT(in):: self CLASS(meshEdge1DRad), INTENT(in):: self

93
src/modules/mesh/1DRad/moduleMesh1DRadBoundary.f90
View file

@ -1,93 +0,0 @@
SUBMODULE (moduleMesh1DRad) moduleMesh1DRadBoundary
USE moduleMesh1DRad
CONTAINS
SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
SELECT TYPE(edge)
TYPE IS(meshEdge1DRad)
part%v(1) = -part%v(1)
part%r(1) = 2.D0*edge%r - part%r(1)
END SELECT
part%n_in = .TRUE.
END SUBROUTINE reflection
SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: rEdge(1) !Position of particle in the edge
REAL(8):: d !Distance from particle to edge
SELECT TYPE(edge)
TYPE IS(meshEdge1DRad)
rEdge(1) = edge%r
d = DABS(part%r(1) - rEdge(1))
IF (d > 0.D0) THEN
part%weight = part%weight / d
part%r(1) = rEdge(1)
END IF
IF (ASSOCIATED(edge%e1)) THEN
CALL edge%e1%scatter(part)
ELSE
CALL edge%e2%scatter(part)
END IF
END SELECT
part%n_in = .FALSE.
END SUBROUTINE absorption
SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
part%n_in = .FALSE.
END SUBROUTINE transparent
SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
USE moduleBoundary
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
!Modifies particle velocity according to wall temperature
SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
TYPE IS(boundaryWallTemperature)
part%v(1) = part%v(1) + bound%vTh*randomMaxwellian()
END SELECT
SELECT TYPE(edge)
TYPE IS(meshEdge1DRad)
part%v(1) = -part%v(1)
part%r(1) = 2.D0*edge%r - part%r(1)
END SELECT
END SUBROUTINE wallTemperature
END SUBMODULE moduleMesh1DRadBoundary

7
src/modules/mesh/2DCart/makefile
View file

@ -1,11 +1,8 @@
all : moduleMesh2DCart.o moduleMesh2DCartBoundary.o moduleMesh2DCartRead.o all : moduleMesh2DCart.o moduleMesh2DCartRead.o
moduleMesh2DCart.o: moduleMesh2DCart.f90 moduleMesh2DCart.o: moduleMesh2DCart.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh2DCartBoundary.o: moduleMesh2DCart.o moduleMesh2DCartBoundary.f90 moduleMesh2DCartRead.o: moduleMesh2DCart.o moduleMesh2DCartRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh2DCartRead.o: moduleMesh2DCart.o moduleMesh2DCartBoundary.o moduleMesh2DCartRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

62
src/modules/mesh/2DCart/moduleMesh2DCart.f90
View file

@ -4,6 +4,7 @@
! z == unused ! z == unused
MODULE moduleMesh2DCart MODULE moduleMesh2DCart
USE moduleMesh USE moduleMesh
USE moduleMeshBoundary
IMPLICIT NONE IMPLICIT NONE
!Values for Gauss integral !Values for Gauss integral
@ -31,50 +32,11 @@ MODULE moduleMesh2DCart
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initEdge2DCart PROCEDURE, PASS:: init => initEdge2DCart
PROCEDURE, PASS:: getNodes => getNodes2DCart PROCEDURE, PASS:: getNodes => getNodes2DCart
PROCEDURE, PASS:: intersection => intersection2DCartEdge
PROCEDURE, PASS:: randPos => randPosEdge PROCEDURE, PASS:: randPos => randPosEdge
END TYPE meshEdge2DCart END TYPE meshEdge2DCart
!Boundary functions defined in the submodule Boundary
INTERFACE
MODULE SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE reflection
MODULE SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE absorption
MODULE SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE wallTemperature
MODULE SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE transparent
END INTERFACE
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCart TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCart
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ2DCart PROCEDURE, PASS:: detJac => detJ2DCart
@ -295,6 +257,23 @@ MODULE moduleMesh2DCart
END FUNCTION getNodes2DCart END FUNCTION getNodes2DCart
PURE FUNCTION intersection2DCartEdge(self, r0, v0) RESULT(r)
IMPLICIT NONE
CLASS(meshEdge2DCart), INTENT(in):: self
REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
REAL(8), DIMENSION(1:3):: r
REAL(8), DIMENSION(1:3):: rS !base point of surface
REAL(8):: d
rS = (/ self%x(1), self%y(1), 0.D0 /)
d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal)
r = r0 + v0*d
END FUNCTION intersection2DCartEdge
!Calculates a random position in edge !Calculates a random position in edge
FUNCTION randPosEdge(self) RESULT(r) FUNCTION randPosEdge(self) RESULT(r)
USE moduleRandom USE moduleRandom
@ -365,7 +344,7 @@ MODULE moduleMesh2DCart
self%arNodes = 0.D0 self%arNodes = 0.D0
!2D 1 point Gauss Quad Integral !2D 1 point Gauss Quad Integral
xi = 0.D0 xi = 0.D0
detJ = self%detJac(xi)*4.D0 !4*2*pi detJ = self%detJac(xi)*4.D0 !4
fPsi = self%fPsi(xi) fPsi = self%fPsi(xi)
self%volume = detJ self%volume = detJ
self%arNodes = fPsi*detJ self%arNodes = fPsi*detJ
@ -930,7 +909,6 @@ MODULE moduleMesh2DCart
REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates
REAL(8):: invJ(1:2,1:2), detJ REAL(8):: invJ(1:2,1:2), detJ
REAL(8):: phi(1:3) REAL(8):: phi(1:3)
REAL(8):: dummy
REAL(8):: EF(1:3) REAL(8):: EF(1:3)
phi = (/self%n1%emData%phi, & phi = (/self%n1%emData%phi, &

154
src/modules/mesh/2DCart/moduleMesh2DCartBoundary.f90
View file

@ -1,154 +0,0 @@
!moduleMesh2DCartBoundary: Boundary functions for cylindrical coordinates
SUBMODULE (moduleMesh2DCart) moduleMesh2DCartBoundary
USE moduleMesh2DCart
CONTAINS
SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
!TODO: Try to do this without select
SELECT TYPE(edge)
TYPE IS(meshEdge2DCart)
edgeNorm = DSQRT((edge%y(2)-edge%y(1))**2 + (edge%x(2)-edge%x(1))**2)
cosT = (edge%x(2)-edge%x(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - edge%x(1);
rp(2) = part%r(2) - edge%y(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) + edge%x(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%y(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
END SELECT
part%n_in = .TRUE.
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: rEdge(1:2) !Position of particle projected to the edge
REAL(8):: a, b, c
REAL(8):: a2b2
REAL(8):: d !Distance from particle to edge
SELECT TYPE(edge)
TYPE IS(meshEdge2DCart)
a = (edge%x(1) - edge%x(2))
b = (edge%y(1) - edge%y(2))
c = edge%x(1)*edge%y(2) - edge%x(2)*edge%y(1)
a2b2 = a**2 + b**2
rEdge(1) = (b*( b*part%r(1) - a*part%r(2)) - a*c)/a2b2
rEdge(2) = (a*(-b*part%r(1) + a*part%r(2)) - b*c)/a2b2
d = NORM2(rEdge - part%r(1:2))
!Reduce weight of particle by the distance to the edge and move it to the edge
IF (d > 0.D0) THEN
part%weight = part%weight / d
part%r(1:2) = rEdge
END IF
!Scatter particle in associated volume
IF (ASSOCIATED(edge%e1)) THEN
CALL edge%e1%scatter(part)
ELSE
CALL edge%e2%scatter(part)
END IF
END SELECT
!Remove particle from the domain
part%n_in = .FALSE.
END SUBROUTINE absorption
!Transparent boundary condition
SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
!Removes particle from domain
part%n_in = .FALSE.
END SUBROUTINE transparent
!Wall with temperature
SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
USE moduleBoundary
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
INTEGER:: i
!Modifies particle velocity according to wall temperature
SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
TYPE IS(boundaryWallTemperature)
DO i = 1, 3
part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
END DO
END SELECT
!Reflects particle in the edge
SELECT TYPE(edge)
TYPE IS(meshEdge2DCart)
edgeNorm = DSQRT((edge%y(2)-edge%y(1))**2 + (edge%x(2)-edge%x(1))**2)
cosT = (edge%x(2)-edge%x(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - edge%x(1);
rp(2) = part%r(2) - edge%y(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) + edge%x(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%y(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
END SELECT
part%n_in = .TRUE.
END SUBROUTINE wallTemperature
END SUBMODULE moduleMesh2DCartBoundary

2
src/modules/mesh/2DCart/moduleMesh2DCartRead.f90
View file

@ -65,7 +65,7 @@ MODULE moduleMesh2DCartRead
ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes)) ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
self%K = 0.D0 self%K = 0.D0
self%IPIV = 0 self%IPIV = 0
!Read nodes cartesian coordinates (x=x, y=y, z=null) !Read node cartesian coordinates (x=x, y=y, z=null)
DO e=1, self%numNodes DO e=1, self%numNodes
READ(10, *) n, x, y READ(10, *) n, x, y
ALLOCATE(meshNode2DCart:: self%nodes(n)%obj) ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)

7
src/modules/mesh/2DCyl/makefile
View file

@ -1,11 +1,8 @@
all : moduleMesh2DCyl.o moduleMesh2DCylBoundary.o moduleMesh2DCylRead.o all : moduleMesh2DCyl.o moduleMesh2DCylRead.o
moduleMesh2DCyl.o: moduleMesh2DCyl.f90 moduleMesh2DCyl.o: moduleMesh2DCyl.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh2DCylBoundary.o: moduleMesh2DCyl.o moduleMesh2DCylBoundary.f90 moduleMesh2DCylRead.o: moduleMesh2DCyl.o moduleMesh2DCylRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh2DCylRead.o: moduleMesh2DCyl.o moduleMesh2DCylBoundary.o moduleMesh2DCylRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

70
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -4,6 +4,7 @@
! z == theta (unused) ! z == theta (unused)
MODULE moduleMesh2DCyl MODULE moduleMesh2DCyl
USE moduleMesh USE moduleMesh
USE moduleMeshBoundary
IMPLICIT NONE IMPLICIT NONE
!Values for Gauss integral !Values for Gauss integral
@ -31,59 +32,11 @@ MODULE moduleMesh2DCyl
CONTAINS CONTAINS
PROCEDURE, PASS:: init => initEdge2DCyl PROCEDURE, PASS:: init => initEdge2DCyl
PROCEDURE, PASS:: getNodes => getNodes2DCyl PROCEDURE, PASS:: getNodes => getNodes2DCyl
PROCEDURE, PASS:: intersection => intersection2DCylEdge
PROCEDURE, PASS:: randPos => randPosEdge PROCEDURE, PASS:: randPos => randPosEdge
END TYPE meshEdge2DCyl END TYPE meshEdge2DCyl
!Boundary functions defined in the submodule Boundary
INTERFACE
MODULE SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE reflection
MODULE SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE absorption
MODULE SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE wallTemperature
MODULE SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE transparent
MODULE SUBROUTINE symmetryAxis(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE symmetryAxis
END INTERFACE
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCyl TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCyl
CONTAINS CONTAINS
PROCEDURE, PASS:: detJac => detJ2DCyl PROCEDURE, PASS:: detJac => detJ2DCyl
@ -286,6 +239,23 @@ MODULE moduleMesh2DCyl
END FUNCTION getNodes2DCyl END FUNCTION getNodes2DCyl
PURE FUNCTION intersection2DCylEdge(self, r0, v0) RESULT(r)
IMPLICIT NONE
CLASS(meshEdge2DCyl), INTENT(in):: self
REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
REAL(8), DIMENSION(1:3):: r
REAL(8), DIMENSION(1:3):: rS !base point of surface
REAL(8):: d
rS = (/ self%z(1), self%r(1), 0.D0 /)
d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal)
r = r0 + v0*d
END FUNCTION intersection2DCylEdge
!Calculates a random position in edge !Calculates a random position in edge
FUNCTION randPosEdge(self) RESULT(r) FUNCTION randPosEdge(self) RESULT(r)
USE moduleRandom USE moduleRandom
@ -734,6 +704,7 @@ MODULE moduleMesh2DCyl
xii(2) = random( 0.D0, 1.D0) xii(2) = random( 0.D0, 1.D0)
xii(3) = 0.D0 xii(3) = 0.D0
ALLOCATE(fPsi(1:3))
fPsi = self%fPsi(xii) fPsi = self%fPsi(xii)
r(1) = DOT_PRODUCT(fPsi, self%z) r(1) = DOT_PRODUCT(fPsi, self%z)
@ -959,7 +930,6 @@ MODULE moduleMesh2DCyl
REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates
REAL(8):: invJ(1:2,1:2), detJ REAL(8):: invJ(1:2,1:2), detJ
REAL(8):: phi(1:3) REAL(8):: phi(1:3)
REAL(8):: dummy
REAL(8):: EF(1:3) REAL(8):: EF(1:3)
phi = (/self%n1%emData%phi, & phi = (/self%n1%emData%phi, &

164
src/modules/mesh/2DCyl/moduleMesh2DCylBoundary.f90
View file

@ -1,164 +0,0 @@
!moduleMesh2DCylBoundary: Boundary functions for cylindrical coordinates
SUBMODULE (moduleMesh2DCyl) moduleMesh2DCylBoundary
USE moduleMesh2DCyl
CONTAINS
SUBROUTINE reflection(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
!TODO: Try to do this without select
SELECT TYPE(edge)
TYPE IS(meshEdge2DCyl)
edgeNorm = DSQRT((edge%r(2)-edge%r(1))**2 + (edge%z(2)-edge%z(1))**2)
cosT = (edge%z(2)-edge%z(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - edge%z(1);
rp(2) = part%r(2) - edge%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) + edge%z(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%r(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
END SELECT
part%n_in = .TRUE.
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: rEdge(1:2) !Position of particle projected to the edge
REAL(8):: a, b, c
REAL(8):: a2b2
REAL(8):: d !Distance from particle to edge
SELECT TYPE(edge)
TYPE IS(meshEdge2DCyl)
a = (edge%z(1) - edge%z(2))
b = (edge%r(1) - edge%r(2))
c = edge%z(1)*edge%r(2) - edge%z(2)*edge%r(1)
a2b2 = a**2 + b**2
rEdge(1) = (b*( b*part%r(1) - a*part%r(2)) - a*c)/a2b2
rEdge(2) = (a*(-b*part%r(1) + a*part%r(2)) - b*c)/a2b2
d = NORM2(rEdge - part%r(1:2))
!Reduce weight of particle by the distance to the edge and move it to the edge
IF (d > 0.D0) THEN
part%weight = part%weight / d
part%r(1:2) = rEdge
END IF
!Scatter particle in associated volume
IF (ASSOCIATED(edge%e1)) THEN
CALL edge%e1%scatter(part)
ELSE
CALL edge%e2%scatter(part)
END IF
END SELECT
!Remove particle from the domain
part%n_in = .FALSE.
END SUBROUTINE absorption
!Transparent boundary condition
SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
!Removes particle from domain
part%n_in = .FALSE.
END SUBROUTINE transparent
!Wall with temperature
SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
USE moduleBoundary
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: edgeNorm, cosT, sinT, rp(1:2), rpp(1:2), vpp(1:2)
INTEGER:: i
!Modifies particle velocity according to wall temperature
SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
TYPE IS(boundaryWallTemperature)
DO i = 1, 3
part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
END DO
END SELECT
!Reflects particle in the edge
SELECT TYPE(edge)
TYPE IS(meshEdge2DCyl)
edgeNorm = DSQRT((edge%r(2)-edge%r(1))**2 + (edge%z(2)-edge%z(1))**2)
cosT = (edge%z(2)-edge%z(1))/edgeNorm
sinT = DSQRT(1-cosT**2)
rp(1) = part%r(1) - edge%z(1);
rp(2) = part%r(2) - edge%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) + edge%z(1);
part%r(2) = -sinT*rpp(1) + cosT*rpp(2) + edge%r(1);
part%v(1) = cosT*vpp(1) + sinT*vpp(2)
part%v(2) = -sinT*vpp(1) + cosT*vpp(2)
END SELECT
part%n_in = .TRUE.
END SUBROUTINE wallTemperature
!Symmetry axis. Dummy function
SUBROUTINE symmetryAxis(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE symmetryAxis
END SUBMODULE moduleMesh2DCylBoundary

8
src/modules/mesh/3DCart/makefile Normal file
View file

@ -0,0 +1,8 @@
all : moduleMesh3DCart.o moduleMesh3DCartRead.o
moduleMesh3DCart.o: moduleMesh3DCart.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMesh3DCartRead.o: moduleMesh3DCart.o moduleMesh3DCartRead.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

672
src/modules/mesh/3DCart/moduleMesh3DCart.f90 Normal file
View file

@ -0,0 +1,672 @@
!moduleMesh3DCart: 3D Cartesian coordinate system
! x == x
! y == y
! z == z
MODULE moduleMesh3DCart
USE moduleMesh
USE moduleMeshBoundary
IMPLICIT NONE
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode3DCart
!Element coordinates
REAL(8):: x, y, z
CONTAINS
PROCEDURE, PASS:: init => initNode3DCart
PROCEDURE, PASS:: getCoordinates => getCoord3DCart
END TYPE meshNode3DCart
!Triangular surface element
TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge3DCartTria
!Element coordinates
REAL(8):: x(1:3) = 0.D0, y(1:3) = 0.D0, z(1:3) = 0.D0
!Connectivity to nodes
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
CONTAINS
PROCEDURE, PASS:: init => initEdge3DCartTria
PROCEDURE, PASS:: getNodes => getNodes3DCartTria
PROCEDURE, PASS:: intersection => intersection3DCartTria
PROCEDURE, PASS:: randPos => randPosEdgeTria
PROCEDURE, NOPASS:: fPsi => fPsiEdgeTria
END TYPE meshEdge3DCartTria
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol3DCart
CONTAINS
PROCEDURE, PASS:: detJac => detJ3DCart
PROCEDURE, PASS:: invJac => invJ3DCart
PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
END TYPE meshVol3DCart
ABSTRACT INTERFACE
PURE FUNCTION fPsi_interface(xii) RESULT(fPsi)
REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
END FUNCTION fPsi_interface
PURE FUNCTION dPsi_interface(xii) RESULT(dPsi)
REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
END FUNCTION dPsi_interface
PURE SUBROUTINE partialDer_interface(self, dPsi, dx, dy, dz)
IMPORT meshVol3DCart
CLASS(meshVol3DCart), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:,1:)
REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz
END SUBROUTINE partialDer_interface
END INTERFACE
!Tetrahedron volume element
TYPE, PUBLIC, EXTENDS(meshVol3DCart):: meshVol3DCartTetra
!Element Coordinates
REAL(8):: x(1:4) = 0.D0, y(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()
CONTAINS
PROCEDURE, PASS:: init => initVolTetra3DCart
PROCEDURE, PASS:: randPos => randPosVolTetra
PROCEDURE, PASS:: calcVol => volumeTetra
PROCEDURE, NOPASS:: fPsi => fPsiTetra
PROCEDURE, NOPASS:: dPsi => dPsiTetra
PROCEDURE, NOPASS:: dPsiXi1 => dPsiTetraXii1
PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXii2
PROCEDURE, PASS:: partialDer => partialDerTetra
PROCEDURE, PASS:: elemK => elemKTetra
PROCEDURE, PASS:: elemF => elemFTetra
PROCEDURE, NOPASS:: weight => weightTetra
PROCEDURE, NOPASS:: inside => insideTetra
PROCEDURE, PASS:: scatter => scatterTetra
PROCEDURE, PASS:: gatherEF => gatherEFTetra
PROCEDURE, PASS:: getNodes => getNodesTetra
PROCEDURE, PASS:: phy2log => phy2logTetra
PROCEDURE, PASS:: nextElement => nextElementTetra
END TYPE meshVol3DCartTetra
CONTAINS
!NODE FUNCTIONS
!Inits node element
SUBROUTINE initNode3DCart(self, n, r)
USE moduleSpecies
USE moduleRefParam
IMPLICIT NONE
CLASS(meshNode3DCart), INTENT(out):: self
INTEGER, INTENT(in):: n
REAL(8), INTENT(in):: r(1:3)
self%n = n
self%x = r(1)/L_ref
self%y = r(2)/L_ref
self%z = r(3)/L_ref
!Node volume, to be determined in mesh
self%v = 0.D0
!Allocates output:
ALLOCATE(self%output(1:nSpecies))
END SUBROUTINE initNode3DCart
!Get coordinates from node
PURE FUNCTION getCoord3DCart(self) RESULT(r)
IMPLICIT NONE
CLASS(meshNode3DCart), INTENT(in):: self
REAL(8):: r(1:3)
r = (/self%x, self%y, self%z/)
END FUNCTION getCoord3DCart
!SURFACE FUNCTIONS
!Inits surface element
SUBROUTINE initEdge3DCartTria(self, n, p, bt, physicalSurface)
USE moduleSpecies
USE moduleBoundary
USE moduleErrors
IMPLICIT NONE
CLASS(meshEdge3DCartTria), 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, r3
INTEGER:: s
self%n = n
self%n1 => mesh%nodes(p(1))%obj
self%n3 => 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%x = (/r1(1), r2(1), r3(1)/)
self%y = (/r1(2), r2(2), r3(2)/)
self%z = (/r1(3), r2(3), r3(3)/)
!Normal vector
self%normal = (/ (self%y(2)-self%y(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)), &
(self%x(2)-self%x(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%x(3)-self%x(1)), &
(self%x(2)-self%x(1))*(self%y(3)-self%y(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)) /)
!Boundary index
self%boundary => boundary(bt)
ALLOCATE(self%fBoundary(1:nSpecies))
!Assign functions to boundary
DO s = 1, nSpecies
SELECT TYPE(obj => self%boundary%bTypes(s)%obj)
TYPE IS(boundaryAbsorption)
self%fBoundary(s)%apply => absorption
TYPE IS(boundaryReflection)
self%fBoundary(s)%apply => reflection
TYPE IS(boundaryTransparent)
self%fBoundary(s)%apply => transparent
TYPE IS(boundaryWallTemperature)
self%fBoundary(s)%apply => wallTemperature
CLASS DEFAULT
CALL criticalError("Boundary type not defined in this geometry", 'initEdge3DCart')
END SELECT
END DO
!Physical surface
self%physicalSurface = physicalSurface
END SUBROUTINE initEdge3DCartTria
!Get nodes from surface
PURE FUNCTION getNodes3DCartTria(self) RESULT(n)
IMPLICIT NONE
CLASS(meshEdge3DCartTria), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1:3))
n = (/self%n1%n, self%n2%n, self%n3%n/)
END FUNCTION getNodes3DCartTria
PURE FUNCTION intersection3DCartTria(self, r0, v0) RESULT(r)
IMPLICIT NONE
CLASS(meshEdge3DCartTria), INTENT(in):: self
REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
REAL(8), DIMENSION(1:3):: r
REAL(8), DIMENSION(1:3):: rS !base point of surface
REAL(8):: d
rS = (/ self%x(1), self%y(1), self%z(1) /)
d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal)
r = r0 + v0*d
END FUNCTION intersection3DCartTria
!Calculates a random position in the surface
FUNCTION randPosEdgeTria(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge3DCartTria), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8):: fPsi(1:3)
xii = (/random(), random(), 0.D0 /)
fPsi = self%fPsi(xii)
r = (/DOT_PRODUCT(fPsi, self%x), &
DOT_PRODUCT(fPsi, self%y), &
DOT_PRODUCT(fPsi, self%z)/)
END FUNCTION randPosEdgeTria
!Shape functions for triangular surface
PURE FUNCTION fPsiEdgeTria(xii) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:3))
fPsi(1) = 1.D0 - xii(1) - xii(2)
fPsi(2) = xii(1)
fPsi(3) = xii(2)
END FUNCTION fPsiEdgeTria
!VOLUME FUNCTIONS
!TETRA FUNCTIONS
!Inits tetrahedron element
SUBROUTINE initVolTetra3DCart(self, n, p)
USE moduleRefParam
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(out):: self
INTEGER, INTENT(in):: n
INTEGER, INTENT(in):: p(:)
REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 !Positions of each node
REAL(8):: volNodes(1:4) !Volume of each node
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%x = (/r1(1), r2(1), r3(1), r4(1)/)
self%y = (/r1(2), r2(2), r3(2), r4(2)/)
self%z = (/r1(3), r2(3), r3(3), r4(3)/)
!Computes the element volume
CALL self%calcVol()
!Assign proportional volume to each node
!TODO: Review this to apply to all elements in the future
volNodes = self%fPsi((/0.25D0, 0.25D0, 0.25D0/))*self%volume
self%n1%v = self%n1%v + volNodes(1)
self%n2%v = self%n2%v + volNodes(2)
self%n3%v = self%n3%v + volNodes(3)
self%n4%v = self%n4%v + volNodes(4)
self%sigmaVrelMax = sigma_ref/L_ref**2
CALL OMP_INIT_LOCK(self%lock)
END SUBROUTINE initVolTetra3DCart
!Random position in volume tetrahedron
FUNCTION randPosVolTetra(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
xii(1) = random(0.D0, 1.D0)
xii(2) = random(0.D0, 1.D0)
xii(3) = random(0.D0, 1.D0)
ALLOCATE(fPsi(1:4))
fPsi = self%fPsi(xii)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = DOT_PRODUCT(fPsi, self%z)
END FUNCTION randPosVolTetra
!Computes the element volume
PURE SUBROUTINE volumeTetra(self)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(inout):: self
REAL(8):: xii(1:3)
self%volume = 0.D0
xii = (/0.25D0, 0.25D0, 0.25D0/)
self%volume = self%detJac(xii)
END SUBROUTINE volumeTetra
!Computes element functions in point xii
PURE FUNCTION fPsiTetra(xii) RESULT(fPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: fPsi(:)
ALLOCATE(fPsi(1:4))
fPsi(1) = 1.D0 - xii(1) - xii(2) - xii(3)
fPsi(2) = xii(1)
fPsi(3) = xii(2)
fPsi(4) = xii(3)
END FUNCTION fPsiTetra
!Derivative element function at coordinates xii
PURE FUNCTION dPsiTetra(xii) RESULT(dPsi)
IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: dPsi(:,:)
ALLOCATE(dPsi(1:3,1:4))
dPsi(1,:) = dPsiTetraXii1(xii(2), xii(3))
dPsi(2,:) = dPsiTetraXii2(xii(1), xii(3))
dPsi(3,:) = dPsiTetraXii3(xii(1), xii(2))
END FUNCTION dPsiTetra
!Derivative element function respect to xii1
PURE FUNCTION dPsiTetraXii1(xii2, xii3) RESULT(dPsiXii1)
IMPLICIT NONE
REAL(8), INTENT(in):: xii2, xii3
REAL(8):: dPsiXii1(1:4)
dPsiXii1(1) = -1.D0
dPsiXii1(2) = 1.D0
dPsiXii1(3) = 0.D0
dPsiXii1(4) = 0.D0
END FUNCTION dPsiTetraXii1
!Derivative element function respect to xii2
PURE FUNCTION dPsiTetraXii2(xii1, xii3) RESULT(dPsiXii2)
IMPLICIT NONE
REAL(8), INTENT(in):: xii1, xii3
REAL(8):: dPsiXii2(1:4)
dPsiXii2(1) = -1.D0
dPsiXii2(2) = 0.D0
dPsiXii2(3) = 1.D0
dPsiXii2(4) = 0.D0
END FUNCTION dPsiTetraXii2
!Derivative element function respect to xii3
PURE FUNCTION dPsiTetraXii3(xii1, xii2) RESULT(dPsiXii3)
IMPLICIT NONE
REAL(8), INTENT(in):: xii1, xii2
REAL(8):: dPsiXii3(1:4)
dPsiXii3(1) = -1.D0
dPsiXii3(2) = 0.D0
dPsiXii3(3) = 0.D0
dPsiXii3(4) = 1.D0
END FUNCTION dPsiTetraXii3
!Computes the derivatives in global coordinates
PURE SUBROUTINE partialDerTetra(self, dPsi, dx, dy, dz)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: dPsi(1:, 1:)
REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz
dx(1) = DOT_PRODUCT(dPsi(1,:), self%x)
dx(2) = DOT_PRODUCT(dPsi(2,:), self%x)
dx(3) = DOT_PRODUCT(dPsi(3,:), self%x)
dy(1) = DOT_PRODUCT(dPsi(1,:), self%y)
dy(2) = DOT_PRODUCT(dPsi(2,:), self%y)
dy(3) = DOT_PRODUCT(dPsi(3,:), self%y)
dz(1) = DOT_PRODUCT(dPsi(1,:), self%z)
dz(2) = DOT_PRODUCT(dPsi(2,:), self%z)
dz(3) = DOT_PRODUCT(dPsi(3,:), self%z)
END SUBROUTINE partialDerTetra
PURE FUNCTION elemKTetra(self) RESULT(ke)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8):: xii(1:3)
REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
REAL(8):: ke(1:4,1:4)
REAL(8):: invJ(1:3,1:3), detJ
!TODO: One point Gauss integral. Upgrade when possible
ke = 0.D0
xii = (/ 0.25D0, 0.25D0, 0.25D0 /)
dPsi = self%dPsi(xii)
detJ = self%detJac(xii, dPsi)
invJ = self%invJac(xii, dPsi)
fPsi = self%fPsi(xii)
ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ
END FUNCTION elemKTetra
PURE FUNCTION elemFTetra(self, source) RESULT(localF)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: source(1:)
REAL(8), ALLOCATABLE:: localF(:)
REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
REAL(8):: xii(1:3)
REAL(8):: detJ, f
ALLOCATE(localF(1:4))
localF = 0.D0
xii = 0.D0
!TODO: One point Gauss integral. Upgrade when possible
xii = (/ 0.25D0, 0.25D0, 0.25D0 /)
dPsi = self%dPsi(xii)
detJ = self%detJac(xii, dPsi)
fPsi = self%fPsi(xii)
f = DOT_PRODUCT(fPsi, source)
localF = localF + f*fPsi*1.D0*detJ
END FUNCTION elemFTetra
PURE FUNCTION weightTetra(xii) RESULT(w)
IMPLICIT NONE
REAL(8), INTENT(in):: xii(1:3)
REAL(8), ALLOCATABLE:: w(:)
w = fPsiTetra(xii)
END FUNCTION weightTetra
PURE FUNCTION insideTetra(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. &
xi(3) >= 0.D0 .AND. &
1.D0 - xi(1) - xi(2) - xi(3) >= 0.D0
END FUNCTION insideTetra
SUBROUTINE scatterTetra(self, part)
USE moduleOutput
USE moduleSpecies
IMPLICIT NONE
CLASS(meshVol3DCartTetra), 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 scatterTetra
PURE FUNCTION gatherEFTetra(self, xi) RESULT(EF)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
REAL(8):: dPsi(1:3, 1:4)
REAL(8):: dPsiR(1:3, 1:4)
REAL(8):: invJ(1:3, 1:3), 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) = -DOT_PRODUCT(dPsiR(3,:), phi)
END FUNCTION gatherEFTetra
PURE FUNCTION getNodesTetra(self) RESULT(n)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), 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 getNodesTetra
PURE FUNCTION phy2logTetra(self,r) RESULT(xi)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: xi(1:3)
REAL(8):: invJ(1:3, 1:3), detJ
REAL(8):: deltaR(1:3)
REAL(8):: dPsi(1:3, 1:4)
xi = 0.D0
deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /)
dPsi = self%dPsi(xi)
invJ = self%invJac(xi, dPsi)
detJ = self%detJac(xi, dPsi)
xi = MATMUL(invJ, deltaR)/detJ
END FUNCTION phy2logTetra
SUBROUTINE nextElementTetra(self, xi, nextElement)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(in):: self
REAL(8), INTENT(in):: xi(1:3)
CLASS(*), POINTER, INTENT(out):: nextElement
REAL(8):: xiArray(1:4)
INTEGER:: nextInt
!TODO: Review when connectivity
xiArray = (/ xi(3), xi(2), 1.D0 - xi(1) - xi(2) - xi(3), 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
CASE (4)
nextElement => self%e4
END SELECT
END SUBROUTINE nextElementTetra
!COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D
!Computes element Jacobian determinant
PURE FUNCTION detJ3DCart(self, xi, dPsi_in) RESULT(dJ)
IMPLICIT NONE
CLASS(meshVol3DCart), INTENT(in)::self
REAL(8), INTENT(in):: xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:, 1:)
REAL(8):: dJ
REAL(8), ALLOCATABLE:: dPsi(:,:)
REAL(8):: dx(1:3), dy(1:3), dz(1:3)
IF (PRESENT(dPsi_in)) THEN
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx, dy, dz)
dJ = dx(1)*(dy(2)*dz(3) - dy(3)*dz(2)) &
- dx(2)*(dy(1)*dz(3) - dy(3)*dz(1)) &
+ dx(3)*(dy(1)*dz(2) - dy(2)*dz(1))
END FUNCTION detJ3DCart
PURE FUNCTION invJ3DCart(self,xi,dPsi_in) RESULT(invJ)
IMPLICIT NONE
CLASS(meshVol3DCart), 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), DIMENSION(1:3):: dx, dy, dz
REAL(8):: invJ(1:3,1:3)
IF(PRESENT(dPsi_in)) THEN
dPsi=dPsi_in
ELSE
dPsi = self%dPsi(xi)
END IF
CALL self%partialDer(dPsi, dx, dy, dz)
invJ(1,1) = (dy(2)*dz(3) - dy(3)*dz(2))
invJ(1,2) = -(dy(1)*dz(3) - dy(3)*dz(1))
invJ(1,3) = (dy(1)*dz(2) - dy(2)*dz(1))
invJ(2,1) = -(dx(2)*dz(3) - dx(3)*dz(2))
invJ(2,2) = (dx(1)*dz(3) - dx(3)*dz(1))
invJ(2,3) = -(dx(1)*dz(2) - dx(2)*dz(1))
invJ(3,1) = -(dx(2)*dy(3) - dx(3)*dy(2))
invJ(3,2) = (dx(1)*dy(3) - dx(3)*dy(1))
invJ(3,3) = -(dx(1)*dy(2) - dx(2)*dy(1))
END FUNCTION invJ3DCart
END MODULE moduleMesh3DCart

470
src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 Normal file
View file

@ -0,0 +1,470 @@
MODULE moduleMesh3DCartRead
USE moduleMesh
USE moduleMesh3DCart
TYPE, EXTENDS(meshGeneric):: mesh3DCartGeneric
CONTAINS
PROCEDURE, PASS:: init => init3DCartMesh
PROCEDURE, PASS:: readMesh => readMesh3DCartGmsh
END TYPE
INTERFACE connect
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connect
CONTAINS
!Init mesh
SUBROUTINE init3DCartMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(mesh3DCartGeneric), 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.", "init3DCartMesh")
END SELECT
END SUBROUTINE init3DCartMesh
!Read mesh from gmsh file
SUBROUTINE readMesh3DCartGmsh(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(mesh3DCartGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: x, y, z
INTEGER:: p(1:4)
INTEGER:: e = 0, et = 0, n = 0, eTemp = 0, elemType = 0, bt = 0
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Read 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 node cartesian coordinates (x = x, y = y, z = z)
DO e = 1, self%numNodes
READ(10, *) n, x, y, z
ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/x, y, z /))
END DO
!Skip comments
READ(10, *)
READ(10, *)
!Reads total number of elements
READ(10, *) totalNumElem
!conts edges and volume elements
self%numEdges = 0
DO e = 1, totalNumElem
READ(10, *) eTemp, elemType
IF (elemType == 2) 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
!Allocate required arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beggining to read each specific element
DO e = 1, totalNumElem
BACKSPACE(10)
END DO
!Reads surfaces
DO e = 1, self%numEdges
READ(10, *) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE(2)
!Triangular surface
READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1:3)
bt = getBoundaryID(boundaryType)
ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj)
CALL self%edges(e)%obj%init(n, p(1:3), bt, boundaryType)
END SELECT
END DO
!Read and initialize volumes
DO e = 1, self%numVols
READ(10, *) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE(4)
!Tetrahedron element
READ(10, *) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVol3DCartTetra:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
END SELECT
END DO
CLOSE(10)
!Build connectivy 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 surfaces
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 readMesh3DCartGmsh
!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(meshVol3DCartTetra)
!Element A is a tetrahedron
SELECT TYPE(elemB)
TYPE IS(meshVol3DCartTetra)
!Element B is a tetrahedron
CALL connectedTetraTetra(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(meshEdge3DCartTria)
SELECT TYPE(elemA)
TYPE IS(meshVol3DCartTetra)
!Element A is a tetrahedron
CALL connectedTetraEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
SUBROUTINE connectedTetraTetra(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemB
!TODO: Try to find a much clear way to do this
!Check surface 1
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
IF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n1%n)) THEN
elemA%e1 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n2%n)) THEN
elemA%e1 => elemB
elemB%e2 => elemA
ELSEIF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n1%n)) THEN
elemA%e1 => elemB
elemB%e3 => elemA
ELSEIF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n1%n)) THEN
elemA%e1 => elemB
elemB%e4 => elemA
END IF
END IF
!Check surface 2
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
IF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e2 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n2%n)) THEN
elemA%e2 => elemB
elemB%e2 => elemA
ELSEIF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n2%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e2 => elemB
elemB%e3 => elemA
ELSEIF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n2%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n2%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n2%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e2 => elemB
elemB%e4 => elemA
END IF
END IF
!Check surface 3
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF ((elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e3 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n2%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n2%n)) THEN
elemA%e3 => elemB
elemB%e2 => elemA
ELSEIF ((elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n2%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n2%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e3 => elemB
elemB%e3 => elemA
ELSEIF ((elemA%n2%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n2%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n2%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e3 => elemB
elemB%e4 => elemA
END IF
END IF
!Check surface 4
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
IF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e4 => elemB
elemB%e1 => elemA
ELSEIF ((elemA%n1%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n2%n)) THEN
elemA%e4 => elemB
elemB%e2 => elemA
ELSEIF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n2%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n2%n) .OR. &
(elemA%n1%n == elemB%n2%n .AND. &
elemA%n3%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e4 => elemB
elemB%e3 => elemA
ELSEIF ((elemA%n1%n == elemB%n1%n .AND. &
elemA%n3%n == elemB%n3%n .AND. &
elemA%n4%n == elemB%n4%n) .OR. &
(elemA%n1%n == elemB%n4%n .AND. &
elemA%n3%n == elemB%n1%n .AND. &
elemA%n4%n == elemB%n3%n) .OR. &
(elemA%n1%n == elemB%n3%n .AND. &
elemA%n3%n == elemB%n4%n .AND. &
elemA%n4%n == elemB%n1%n)) THEN
elemA%e4 => elemB
elemB%e4 => elemA
END IF
END IF
END SUBROUTINE connectedTetraTetra
SUBROUTINE connectedTetraEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB
END SUBROUTINE connectedTetraEdge
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(:)
END SUBROUTINE constructGlobalK
END MODULE moduleMesh3DCartRead

7
src/modules/mesh/makefile
View file

@ -1,4 +1,7 @@
all: moduleMesh.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o all: moduleMesh.o moduleMeshBoundary.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
3DCart.o:
$(MAKE) -C 3DCart all
2DCyl.o: 2DCyl.o:
$(MAKE) -C 2DCyl all $(MAKE) -C 2DCyl all
@ -15,3 +18,5 @@ all: moduleMesh.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o
moduleMesh.o: moduleMesh.f90 moduleMesh.o: moduleMesh.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90
$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@

18
src/modules/mesh/moduleMesh.f90
View file

@ -70,11 +70,13 @@ MODULE moduleMesh
CONTAINS CONTAINS
PROCEDURE(initEdge_interface), DEFERRED, PASS:: init PROCEDURE(initEdge_interface), DEFERRED, PASS:: init
PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes
PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection
PROCEDURE(randPosEdge_interface), DEFERRED, PASS:: randPos PROCEDURE(randPosEdge_interface), DEFERRED, PASS:: randPos
END TYPE meshEdge END TYPE meshEdge
ABSTRACT INTERFACE ABSTRACT INTERFACE
!Inits the edge parameters
SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface) SUBROUTINE initEdge_interface(self, n, p, bt, physicalSurface)
IMPORT:: meshEdge IMPORT:: meshEdge
@ -86,13 +88,24 @@ MODULE moduleMesh
END SUBROUTINE initEdge_interface END SUBROUTINE initEdge_interface
!Get nodes index from node
PURE FUNCTION getNodesEdge_interface(self) RESULT(n) PURE FUNCTION getNodesEdge_interface(self) RESULT(n)
IMPORT:: meshEdge IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self CLASS(meshEdge), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:) INTEGER, ALLOCATABLE:: n(:)
END FUNCTION END FUNCTION getNodesEdge_interface
!Returns the intersecction between an edge and a line defined by point r0 and vector v0
PURE FUNCTION intersectionEdge_interface(self, r0, v0) RESULT(r)
IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self
REAL(8), INTENT(in), DIMENSION(1:3):: r0, v0
REAL(8):: r(1:3)
END FUNCTION intersectionEdge_interface
!Returns a random position in the edge
FUNCTION randPosEdge_interface(self) RESULT(r) FUNCTION randPosEdge_interface(self) RESULT(r)
IMPORT:: meshEdge IMPORT:: meshEdge
CLASS(meshEdge), INTENT(in):: self CLASS(meshEdge), INTENT(in):: self
@ -391,11 +404,8 @@ MODULE moduleMesh
IMPLICIT NONE IMPLICIT NONE
CLASS(meshVol), INTENT(inout):: self CLASS(meshVol), INTENT(inout):: self
INTEGER:: modCollisions !Remain of current iteration and everyCollisions
INTEGER:: iterToCollisions !Number of iterations from current to next collision
INTEGER:: nPart !Number of particles inside the cell INTEGER:: nPart !Number of particles inside the cell
REAL(8):: pMax !Maximum probability of collision REAL(8):: pMax !Maximum probability of collision
INTEGER:: nCollIter !Number of collisions to be computed in this iteration
INTEGER:: rnd !random index INTEGER:: rnd !random index
TYPE(particle), POINTER:: part_i, part_j TYPE(particle), POINTER:: part_i, part_j
INTEGER:: n !collision INTEGER:: n !collision

121
src/modules/mesh/moduleMeshBoundary.f90 Normal file
View file

@ -0,0 +1,121 @@
!moduleMeshBoundary: Boundary functions
MODULE moduleMeshBoundary
USE moduleMesh
CONTAINS
SUBROUTINE reflection(edge, part)
USE moduleCaseParam
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
!rp = intersection between particle and edge
!rpp = final position of particle
!vpp = final velocity of particle
REAL(8), DIMENSION(1:3):: rp, rpp, vpp
REAL(8):: taup !time step for reflecting process
!Reflect particle velocity
vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal
!Computes the intersection between particle and surface
rp = edge%intersection(part%r, part%v)
!Computes the reflection time step
taup = NORM2(part%r - rp)*tau(part%sp)
!New position of particle
rpp = rp + vpp*taup
!assign new parameters to particle
part%r = rpp
part%v = vpp
part%n_in = .TRUE.
END SUBROUTINE reflection
!Absoption in a surface
SUBROUTINE absorption(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
REAL(8):: rpp(1:3) !Position of particle projected to the edge
REAL(8):: d !Distance from particle to edge
rpp = edge%intersection(part%r, part%v)
d = NORM2(rpp - part%r)
IF (d >= 0.D0) THEN
part%weight = part%weight/d
END IF
!Assign new position to particle
part%r = rpp
!Remove particle from the domain
part%n_in = .FALSE.
!Scatter particle in associated volume
IF (ASSOCIATED(edge%e1)) THEN
CALL edge%e1%scatter(part)
ELSE
CALL edge%e2%scatter(part)
END IF
END SUBROUTINE absorption
!Transparent boundary condition
SUBROUTINE transparent(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
!Removes particle from domain
part%n_in = .FALSE.
END SUBROUTINE transparent
!Wall with temperature
SUBROUTINE wallTemperature(edge, part)
USE moduleSpecies
USE moduleBoundary
USE moduleRandom
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
INTEGER:: i
!Modifies particle velocity according to wall temperature
SELECT TYPE(bound => edge%boundary%bTypes(part%sp)%obj)
TYPE IS(boundaryWallTemperature)
DO i = 1, 3
part%v(i) = part%v(i) + bound%vTh*randomMaxwellian()
END DO
END SELECT
CALL reflection(edge, part)
END SUBROUTINE wallTemperature
!Symmetry axis. Dummy function
SUBROUTINE symmetryAxis(edge, part)
USE moduleSpecies
IMPLICIT NONE
CLASS(meshEdge), INTENT(inout):: edge
CLASS(particle), INTENT(inout):: part
END SUBROUTINE symmetryAxis
END MODULE moduleMeshBoundary

4
src/modules/moduleCollisions.f90
View file

@ -277,7 +277,7 @@ MODULE moduleCollisions
REAL(8), INTENT(in):: sigmaVrelMax REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew REAL(8), INTENT(inout):: sigmaVrelMaxNew
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
TYPE(particle), POINTER:: electron, neutral TYPE(particle), POINTER:: electron => NULL(), neutral => NULL()
TYPE(particle), POINTER:: newElectron TYPE(particle), POINTER:: newElectron
REAL(8):: vRel, eRel REAL(8):: vRel, eRel
REAL(8):: sigmaVrel REAL(8):: sigmaVrel
@ -408,7 +408,7 @@ MODULE moduleCollisions
REAL(8), INTENT(in):: sigmaVrelMax REAL(8), INTENT(in):: sigmaVrelMax
REAL(8), INTENT(inout):: sigmaVrelMaxNew REAL(8), INTENT(inout):: sigmaVrelMaxNew
TYPE(particle), INTENT(inout), TARGET:: part_i, part_j TYPE(particle), INTENT(inout), TARGET:: part_i, part_j
TYPE(particle), POINTER:: electron, ion TYPE(particle), POINTER:: electron => NULL(), ion => NULL()
REAL(8):: vRel, eRel REAL(8):: vRel, eRel
REAL(8):: sigmaVrel REAL(8):: sigmaVrel
REAL(8), DIMENSION(1:3):: vp_i REAL(8), DIMENSION(1:3):: vp_i