diff --git a/doc/user-manual/fpakc_UserManual.pdf b/doc/user-manual/fpakc_UserManual.pdf index e4bd94e..fa6ad94 100644 Binary files a/doc/user-manual/fpakc_UserManual.pdf and b/doc/user-manual/fpakc_UserManual.pdf differ diff --git a/doc/user-manual/fpakc_UserManual.tex b/doc/user-manual/fpakc_UserManual.tex index 567e1f5..efb090e 100644 --- a/doc/user-manual/fpakc_UserManual.tex +++ b/doc/user-manual/fpakc_UserManual.tex @@ -128,15 +128,18 @@ When a 2D cylindrical geometry is used ($z$, $r$), a Boris solver\cite{boris1970relativistic} is used to move particles accounting for the effect of the symmetry axis. This pusher removes the issue with particles going to infinite velocity when $r \rightarrow 0$ by pushing the particles in Cartesian space and then converting it to $r-z$ geometry. Velocity in the $\theta$ direction is updated for collision processes, although the dynamic in the angular direction is assumed as symmetric. - - Cross-sections are read from files. - These files contains two columns: one for the relative energy (in $\unit{eV}$) and another one for the cross-section (in $\unit{m^{-2}}$). - + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - \subsection{1D Cartesian pusher} + \subsection{2D Cartesian pusher} + Moving particles in a simple 2D Cartesian space. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{1D Radial pusher} + Same implementation as 2D cylindrical pusher but direction $z$ is ignored. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + \subsection{1D Cartesian pusher} + Moving particles in a simple 1D Cartesian space. Same implementation as in 2D Cartesian but $y$ direction is ignored. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Find new cell} @@ -338,12 +341,14 @@ make Type of geometry. Current accepted vaules are \begin{itemize} - \item \textbf{2DCyl}: Two-dimensional grid (z-r) with symmetry axis at $r = 0$. + \item \textbf{2DCyl}: Two-dimensional grid ($z \hyphen r$) with symmetry axis at $r = 0$. For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $z$ and $r$ respectively. - \item \textbf{1DCart}: One-dimensional grid (x) in Cartesian coordinates - For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $x$. - \item \textbf{1DRad}: One-dimensional grid (r) in radial coordinates + \item \textbf{2DCart}: Two-dimensional grid ($x \hyphen y$) in Cartesian coordinates.. + For \Gls{gmsh} mesh format, the coordinates $x$ and $y$ correspond to $x$ and $y$ respectively. + \item \textbf{1DRad}: One-dimensional grid ($r$) in radial coordinates For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $r$. + \item \textbf{1DCart}: One-dimensional grid ($x$) in Cartesian coordinates + For \Gls{gmsh} mesh format, the coordinates $x$ corresponds to $x$. \end{itemize} \item \textbf{meshType}: Character. Format of mesh file. @@ -494,10 +499,14 @@ make Array dimension 'number of species'. Indicates the type of pusher used for each species: \begin{itemize} - \item \textbf{2DCylNeutral}: Pushes particles in a 2D z-r space without any external force. - \item \textbf{2DCylCharged}: Pushes particles in a 2D z-r space including the effect of the electrostatic field. - \item \textbf{1DCartCharged}: Pushes particles in a 1D Cartesian space accounting the the electrostatic field. - \item \textbf{1DRadCharged}: Pushes particles in a 1D cylindrical space (r) accounting the the electrostatic field. + \item \textbf{2DCylNeutral}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) without any external force. + \item \textbf{2DCylCharged}: Pushes particles in a 2D cylindrical space ($z \hyphen r$) including the effect of the electrostatic field. + \item \textbf{2DCartNeutral}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) without any external force. + \item \textbf{2DCartCharged}: Pushes particles in a 2D Cartesian space ($x \hyphen y$) including the effect of the electrostatic field. + \item \textbf{1DRadNeutral}: Pushes particles in a 1D cylindrical space ($r$) without any external force. + \item \textbf{1DRadCharged}: Pushes particles in a 1D cylindrical space ($r$) accounting the the electrostatic field. + \item \textbf{1DCartNeutral}: Pushes particles in a 1D Cartesian space ($x$) without any external force. + \item \textbf{1DCartCharged}: Pushes particles in a 1D Cartesian space ($x$) accounting the the electrostatic field. \end{itemize} \item \textbf{WeightingScheme}: Character. Indicates the variable weighting scheme to be used in the simulation. diff --git a/src/makefile b/src/makefile index 3e151a3..004955c 100644 --- a/src/makefile +++ b/src/makefile @@ -4,7 +4,8 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleCompTime.o $(OBJDIR)/moduleSolv $(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \ $(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \ $(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \ - $(OBJDIR)/moduleMeshCyl.o $(OBJDIR)/moduleMeshCylRead.o $(OBJDIR)/moduleMeshCylBoundary.o \ + $(OBJDIR)/moduleMesh2DCyl.o $(OBJDIR)/moduleMesh2DCylRead.o $(OBJDIR)/moduleMesh2DCylBoundary.o \ + $(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o $(OBJDIR)/moduleMesh2DCartBoundary.o \ $(OBJDIR)/moduleMesh1DCart.o $(OBJDIR)/moduleMesh1DCartRead.o $(OBJDIR)/moduleMesh1DCartBoundary.o \ $(OBJDIR)/moduleMesh1DRad.o $(OBJDIR)/moduleMesh1DRadRead.o $(OBJDIR)/moduleMesh1DRadBoundary.o diff --git a/src/modules/mesh/2DCart/makefile b/src/modules/mesh/2DCart/makefile new file mode 100644 index 0000000..312916b --- /dev/null +++ b/src/modules/mesh/2DCart/makefile @@ -0,0 +1,11 @@ +all : moduleMesh2DCart.o moduleMesh2DCartBoundary.o moduleMesh2DCartRead.o + +moduleMesh2DCart.o: moduleMesh2DCart.f90 + $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ + +moduleMesh2DCartBoundary.o: moduleMesh2DCart.o moduleMesh2DCartBoundary.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)/$@ + diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 new file mode 100644 index 0000000..8718804 --- /dev/null +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -0,0 +1,1057 @@ +!moduleMesh2DCart: 2D Cartesian coordinate system +! x == x +! y == y +! z == unused +MODULE moduleMesh2DCart + 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):: meshNode2DCart + !Element coordinates + REAL(8):: x = 0.D0, y = 0.D0 + CONTAINS + PROCEDURE, PASS:: init => initNode2DCart + PROCEDURE, PASS:: getCoordinates => getCoord2DCart + + END TYPE meshNode2DCart + + TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge2DCart + !Element coordinates + REAL(8):: x(1:2) = 0.D0, y(1:2) = 0.D0 + !Connectivity to nodes + CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() + CONTAINS + PROCEDURE, PASS:: init => initEdge2DCart + PROCEDURE, PASS:: getNodes => getNodes2DCart + PROCEDURE, PASS:: randPos => randPosEdge + + 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 + CONTAINS + PROCEDURE, PASS:: detJac => detJ2DCart + PROCEDURE, PASS:: invJac => invJ2DCart + PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi + PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi + PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer + + END TYPE meshVol2DCart + + 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, dy) + IMPORT meshVol2DCart + CLASS(meshVol2DCart), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy + + END SUBROUTINE partialDer_interface + + END INTERFACE + + !Quadrilateral volume element + TYPE, PUBLIC, EXTENDS(meshVol2DCart):: meshVol2DCartQuad + !Element coordinates + REAL(8):: x(1:4) = 0.D0, y(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 => initVolQuad2DCart + PROCEDURE, PASS:: randPos => randPosVolQuad + 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 + + END TYPE meshVol2DCartQuad + + !Triangular volume element + TYPE, PUBLIC, EXTENDS(meshVol2DCart):: meshVol2DCartTria + !Element coordinates + REAL(8):: x(1:3) = 0.D0, y(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 + + CONTAINS + PROCEDURE, PASS:: init => initVolTria2DCart + PROCEDURE, PASS:: randPos => randPosVolTria + 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 + + END TYPE meshVol2DCartTria + + CONTAINS + !NODE FUNCTIONS + !Inits node element + SUBROUTINE initNode2DCart(self, n, r) + USE moduleSpecies + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshNode2DCart), 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 + !Node volume, to be determined in mesh + self%v = 0.D0 + + !Allocates output: + ALLOCATE(self%output(1:nSpecies)) + + END SUBROUTINE initNode2DCart + + !Get coordinates from node + PURE FUNCTION getCoord2DCart(self) RESULT(r) + IMPLICIT NONE + + CLASS(meshNode2DCart), INTENT(in):: self + REAL(8):: r(1:3) + + r = (/self%x, self%y, 0.D0/) + + END FUNCTION getCoord2DCart + + !EDGE FUNCTIONS + !Inits edge element + SUBROUTINE initEdge2DCart(self, n, p, bt, physicalSurface) + USE moduleSpecies + USE moduleBoundary + USE moduleErrors + IMPLICIT NONE + + CLASS(meshEdge2DCart), 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 + INTEGER:: s + + 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)/) + self%y = (/r1(2), r2(2)/) + !Normal vector + self%normal = (/ self%y(2)-self%y(1), & + self%x(2)-self%x(1), & + 0.D0 /) + !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", 'initEdge2DCart') + + END SELECT + + END DO + + !Physical surface + self%physicalSurface = physicalSurface + + END SUBROUTINE initEdge2DCart + + !Random position in quadrilateral volume + FUNCTION randPosVolQuad(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(in):: self + REAL(8):: r(1:3) + REAL(8):: xii(1:3) + REAL(8), ALLOCATABLE:: fPsi(:) + + xii(1) = random(-1.D0, 1.D0) + xii(2) = random(-1.D0, 1.D0) + xii(3) = 0.D0 + + fPsi = self%fPsi(xii) + + r(1) = DOT_PRODUCT(fPsi, self%x) + r(2) = DOT_PRODUCT(fPsi, self%y) + r(3) = 0.D0 + + END FUNCTION randposVolQuad + + !Get nodes from edge + PURE FUNCTION getNodes2DCart(self) RESULT(n) + IMPLICIT NONE + + CLASS(meshEdge2DCart), INTENT(in):: self + INTEGER, ALLOCATABLE:: n(:) + + ALLOCATE(n(1:2)) + n = (/self%n1%n, self%n2%n /) + + END FUNCTION getNodes2DCart + + !Calculates a random position in edge + FUNCTION randPosEdge(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshEdge2DCart), INTENT(in):: self + REAL(8):: rnd + REAL(8):: r(1:3) + REAL(8):: p1(1:2), p2(1:2) + + rnd = random() + + p1 = (/self%x(1), self%y(1) /) + p2 = (/self%x(2), self%y(2) /) + r(1:2) = (1.D0 - rnd)*p1 + rnd*p2 + r(3) = 0.D0 + + END FUNCTION randPosEdge + + !VOLUME FUNCTIONS + !QUAD FUNCTIONS + !Inits quadrilateral element + SUBROUTINE initVolQuad2DCart(self, n, p) + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), 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%x = (/r1(1), r2(1), r3(1), r4(1)/) + self%y = (/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 initVolQuad2DCart + + !Computes element area + PURE SUBROUTINE areaQuad(self) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(inout):: self + REAL(8):: 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)*4.D0 !4*2*pi + fPsi = self%fPsi(xi) + self%volume = detJ + self%arNodes = fPsi*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, dx, dy) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy + + dx(1) = DOT_PRODUCT(dPsi(1,:),self%x) + dx(2) = DOT_PRODUCT(dPsi(2,:),self%x) + dy(1) = DOT_PRODUCT(dPsi(1,:),self%y) + dy(2) = DOT_PRODUCT(dPsi(2,:),self%y) + + END SUBROUTINE partialDerQuad + + !Computes element local stiffness matrix + PURE FUNCTION elemKQuad(self) RESULT(ke) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(in):: self + REAL(8):: 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) + ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ + + END DO + END DO + + END FUNCTION elemKQuad + + !Computes the local source vector for a force f + PURE FUNCTION elemFQuad(self, source) RESULT(localF) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: source(1:) + REAL(8), ALLOCATABLE:: localF(:) + REAL(8):: 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) + f = DOT_PRODUCT(fPsi,source) + localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ + + END DO + END DO + + 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(meshVol2DCartQuad), 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(meshVol2DCartQuad), 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(meshVol2DCartQuad), 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(meshVol2DCartQuad), 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%x)-r(1) + f(2) = DOT_PRODUCT(fPsi,self%y)-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(meshVol2DCartQuad), 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 + + !TRIA ELEMENT + !Init tria element + SUBROUTINE initVolTria2DCart(self, n, p) + USE moduleRefParam + IMPLICIT NONE + + CLASS(meshVol2DCartTria), 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%x = (/r1(1), r2(1), r3(1)/) + self%y = (/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) + + self%sigmaVrelMax = sigma_ref/L_ref**2 + + CALL OMP_INIT_LOCK(self%lock) + + END SUBROUTINE initVolTria2DCart + + !Random position in quadrilateral volume + FUNCTION randPosVolTria(self) RESULT(r) + USE moduleRandom + IMPLICIT NONE + + CLASS(meshVol2DCartTria), 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) = 0.D0 + + fPsi = self%fPsi(xii) + + r(1) = DOT_PRODUCT(fPsi, self%x) + r(2) = DOT_PRODUCT(fPsi, self%y) + r(3) = 0.D0 + + END FUNCTION randposVolTria + + !Calculates area for triangular element + PURE SUBROUTINE areaTria(self) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), INTENT(inout):: self + REAL(8):: 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)/2.D0 + fPsi = self%fPsi(xi) + self%volume = detJ + self%arNodes = fPsi*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, dx, dy) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: dPsi(1:,1:) + REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy + + dx(1) = DOT_PRODUCT(dPsi(1,:),self%x) + dx(2) = DOT_PRODUCT(dPsi(2,:),self%x) + dy(1) = DOT_PRODUCT(dPsi(1,:),self%y) + dy(2) = DOT_PRODUCT(dPsi(2,:),self%y) + + END SUBROUTINE partialDerTria + + !Computes element local stiffness matrix + PURE FUNCTION elemKTria(self) RESULT(ke) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), INTENT(in):: self + REAL(8):: 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) + ke = ke + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ + + END DO + + END FUNCTION elemKTria + + !Computes element local source vector + PURE FUNCTION elemFTria(self, source) RESULT(localF) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: source(1:) + REAL(8), ALLOCATABLE:: localF(:) + REAL(8):: fPsi(1:3) + REAL(8):: 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) + f = DOT_PRODUCT(fPsi,source) + localF = localF + f*fPsi*wTria(l)*detJ + + END DO + + 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(meshVol2DCartTria), 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(meshVol2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: dPsi(1:2,1:3) + REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates + REAL(8):: invJ(1:2,1:2), detJ + REAL(8):: phi(1:3) + REAL(8):: dummy + REAL(8):: EF(1:3) + + phi = (/self%n1%emData%phi, & + self%n2%emData%phi, & + self%n3%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 gatherEFTria + + !Gets node indexes from triangular element + PURE FUNCTION getNodesTria(self) RESULT(n) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), 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(meshVol2DCartTria), 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%x(1), r(2) - self%y(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(meshVol2DCartTria), 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 + + !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D + !Computes element Jacobian determinant + PURE FUNCTION detJ2DCart(self, xi, dPsi_in) RESULT(dJ) + IMPLICIT NONE + + CLASS(meshVol2DCart), 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:2), dy(1:2) + + IF(PRESENT(dPsi_in)) THEN + dPsi = dPsi_in + + ELSE + dPsi = self%dPsi(xi) + + END IF + + CALL self%partialDer(dPsi, dx, dy) + dJ = dx(1)*dy(2)-dx(2)*dy(1) + + END FUNCTION detJ2DCart + + !Computes element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ2DCart(self,xi,dPsi_in) RESULT(invJ) + IMPLICIT NONE + + CLASS(meshVol2DCart), 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:2), dy(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, dx, dy) + invJ(1,:) = (/ dy(2), -dx(2) /) + invJ(2,:) = (/ -dy(1), dx(1) /) + + END FUNCTION invJ2DCart + +END MODULE moduleMesh2DCart diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartBoundary.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartBoundary.f90 new file mode 100644 index 0000000..e4cbe1c --- /dev/null +++ b/src/modules/mesh/2DCart/moduleMesh2DCartBoundary.f90 @@ -0,0 +1,154 @@ +!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 diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 new file mode 100644 index 0000000..db537c6 --- /dev/null +++ b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 @@ -0,0 +1,599 @@ +MODULE moduleMesh2DCartRead + USE moduleMesh + USE moduleMesh2DCart + + TYPE, EXTENDS(meshGeneric):: mesh2DCartGeneric + CONTAINS + PROCEDURE, PASS:: init => init2DCartMesh + PROCEDURE, PASS:: readMesh => readMesh2DCartGmsh + + END TYPE + + INTERFACE connected + MODULE PROCEDURE connectedVolVol, connectedVolEdge + + END INTERFACE connected + + CONTAINS + !Init mesh + SUBROUTINE init2DCartMesh(self, meshFormat) + USE moduleMesh + USE moduleErrors + IMPLICIT NONE + + CLASS(mesh2DCartGeneric), 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.", "init2DCartMesh") + + END SELECT + + END SUBROUTINE init2DCartMesh + + !Read mesh from gmsh file + SUBROUTINE readMesh2DCartGmsh(self, filename) + USE moduleBoundary + IMPLICIT NONE + + CLASS(mesh2DCartGeneric), INTENT(inout):: self + CHARACTER(:), ALLOCATABLE, INTENT(in):: filename + REAL(8):: x, y + 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=x, y=y, z=null) + DO e=1, self%numNodes + READ(10, *) n, x, y + ALLOCATE(meshNode2DCart:: self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/x, y, 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) + + ALLOCATE(meshEdge2DCart:: self%edges(e)%obj) + + 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(meshVol2DCartTria:: 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(meshVol2DCartQuad:: 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 readMesh2DCartGmsh + + !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(meshVol2DCartQuad) + !Element A is a quadrilateral + SELECT TYPE(elemB) + TYPE IS(meshVol2DCartQuad) + !Element B is a quadrilateral + CALL connectedQuadQuad(elemA, elemB) + + TYPE IS(meshVol2DCartTria) + !Element B is a triangle + CALL connectedQuadTria(elemA, elemB) + + END SELECT + + TYPE IS(meshVol2DCartTria) + !Element A is a Triangle + SELECT TYPE(elemB) + TYPE IS(meshVol2DCartQuad) + !Element B is a quadrilateral + CALL connectedQuadTria(elemB, elemA) + + TYPE IS(meshVol2DCartTria) + !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(meshEdge2DCart) + SELECT TYPE(elemA) + TYPE IS(meshVol2DCartQuad) + !Element A is a quadrilateral + CALL connectedQuadEdge(elemA, elemB) + + TYPE IS(meshVol2DCartTria) + !Element A is a triangle + CALL connectedTriaEdge(elemA, elemB) + + END SELECT + + END SELECT + + END SUBROUTINE connectedVolEdge + + SUBROUTINE connectedQuadQuad(elemA, elemB) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCartQuad), 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(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCartTria), 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(meshVol2DCartTria), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCartTria), 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(meshVol2DCartQuad), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCart), 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(meshVol2DCartTria), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCart), 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(meshVol2DCartQuad) + 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(meshVol2DCartTria) + 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 moduleMesh2DCartRead diff --git a/src/modules/mesh/2DCyl/makefile b/src/modules/mesh/2DCyl/makefile index 73cd926..fd8e453 100644 --- a/src/modules/mesh/2DCyl/makefile +++ b/src/modules/mesh/2DCyl/makefile @@ -1,11 +1,11 @@ -all : moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.o +all : moduleMesh2DCyl.o moduleMesh2DCylBoundary.o moduleMesh2DCylRead.o -moduleMeshCyl.o: moduleMeshCyl.f90 +moduleMesh2DCyl.o: moduleMesh2DCyl.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMeshCylBoundary.o: moduleMeshCyl.o moduleMeshCylBoundary.f90 +moduleMesh2DCylBoundary.o: moduleMesh2DCyl.o moduleMesh2DCylBoundary.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ -moduleMeshCylRead.o: moduleMeshCyl.o moduleMeshCylBoundary.o moduleMeshCylRead.f90 +moduleMesh2DCylRead.o: moduleMesh2DCyl.o moduleMesh2DCylBoundary.o moduleMesh2DCylRead.f90 $(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@ diff --git a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 similarity index 88% rename from src/modules/mesh/2DCyl/moduleMeshCyl.f90 rename to src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index ad377dc..da1b620 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -1,8 +1,8 @@ -!moduleMeshCyl: 2D axial symmetric extension of generic mesh from GMSH format. +!moduleMesh2DCyl: 2D axial symmetric extension of generic mesh from GMSH format. ! x == z ! y == r ! z == theta (unused) -MODULE moduleMeshCyl +MODULE moduleMesh2DCyl USE moduleMesh IMPLICIT NONE @@ -14,26 +14,26 @@ MODULE moduleMeshCyl 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 + TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCyl !Element coordinates REAL(8):: r = 0.D0, z = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNodeCyl - PROCEDURE, PASS:: getCoordinates => getCoordCyl + PROCEDURE, PASS:: init => initNode2DCyl + PROCEDURE, PASS:: getCoordinates => getCoord2DCyl - END TYPE meshNodeCyl + END TYPE meshNode2DCyl - TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdgeCyl + TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge2DCyl !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:: init => initEdge2DCyl + PROCEDURE, PASS:: getNodes => getNodes2DCyl PROCEDURE, PASS:: randPos => randPosEdge - END TYPE meshEdgeCyl + END TYPE meshEdge2DCyl !Boundary functions defined in the submodule Boundary INTERFACE @@ -84,15 +84,15 @@ MODULE moduleMeshCyl END INTERFACE - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVolCyl + TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCyl CONTAINS - PROCEDURE, PASS:: detJac => detJCyl - PROCEDURE, PASS:: invJac => invJCyl + PROCEDURE, PASS:: detJac => detJ2DCyl + PROCEDURE, PASS:: invJac => invJ2DCyl PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - END TYPE meshVolCyl + END TYPE meshVol2DCyl ABSTRACT INTERFACE PURE FUNCTION fPsi_interface(xi) RESULT(fPsi) @@ -108,8 +108,8 @@ MODULE moduleMeshCyl END FUNCTION dPsi_interface PURE SUBROUTINE partialDer_interface(self, dPsi, dz, dr) - IMPORT meshVolCyl - CLASS(meshVolCyl), INTENT(in):: self + IMPORT meshVol2DCyl + CLASS(meshVol2DCyl), INTENT(in):: self REAL(8), INTENT(in):: dPsi(1:,1:) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr @@ -118,7 +118,7 @@ MODULE moduleMeshCyl END INTERFACE !Quadrilateral volume element - TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylQuad + TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylQuad !Element coordinates REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 !Connectivity to nodes @@ -128,7 +128,7 @@ MODULE moduleMeshCyl REAL(8):: arNodes(1:4) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initVolQuadCyl + PROCEDURE, PASS:: init => initVolQuad2DCyl PROCEDURE, PASS:: randPos => randPosVolQuad PROCEDURE, PASS:: area => areaQuad PROCEDURE, NOPASS:: fPsi => fPsiQuad @@ -146,10 +146,10 @@ MODULE moduleMeshCyl PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: nextElement => nextElementQuad - END TYPE meshVolCylQuad + END TYPE meshVol2DCylQuad !Triangular volume element - TYPE, PUBLIC, EXTENDS(meshVolCyl):: meshVolCylTria + TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylTria !Element coordinates REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0 !Connectivity to nodes @@ -157,11 +157,9 @@ MODULE moduleMeshCyl !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:: init => initVolTria2DCyl PROCEDURE, PASS:: randPos => randPosVolTria PROCEDURE, PASS:: area => areaTria PROCEDURE, NOPASS:: fPsi => fPsiTria @@ -179,51 +177,51 @@ MODULE moduleMeshCyl PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: nextElement => nextElementTria - END TYPE meshVolCylTria + END TYPE meshVol2DCylTria CONTAINS !NODE FUNCTIONS !Inits node element - SUBROUTINE initNodeCyl(self, n, r) + SUBROUTINE initNode2DCyl(self, n, r) USE moduleSpecies USE moduleRefParam IMPLICIT NONE - CLASS(meshNodeCyl), INTENT(out):: self + CLASS(meshNode2DCyl), 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 + self%z = r(1)/L_ref + self%r = r(2)/L_ref !Node volume, to be determined in mesh self%v = 0.D0 !Allocates output: ALLOCATE(self%output(1:nSpecies)) - END SUBROUTINE initNodeCyl + END SUBROUTINE initNode2DCyl !Get coordinates from node - PURE FUNCTION getCoordCyl(self) RESULT(r) + PURE FUNCTION getCoord2DCyl(self) RESULT(r) IMPLICIT NONE - CLASS(meshNodeCyl), INTENT(in):: self + CLASS(meshNode2DCyl), INTENT(in):: self REAL(8):: r(1:3) r = (/self%z, self%r, 0.D0/) - END FUNCTION getCoordCyl + END FUNCTION getCoord2DCyl !EDGE FUNCTIONS !Inits edge element - SUBROUTINE initEdgeCyl(self, n, p, bt, physicalSurface) + SUBROUTINE initEdge2DCyl(self, n, p, bt, physicalSurface) USE moduleSpecies USE moduleBoundary USE moduleErrors IMPLICIT NONE - CLASS(meshEdgeCyl), INTENT(out):: self + CLASS(meshEdge2DCyl), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) INTEGER, INTENT(in):: bt @@ -265,7 +263,7 @@ MODULE moduleMeshCyl self%fBoundary(s)%apply => symmetryAxis CLASS DEFAULT - CALL criticalError("Boundary type not defined in this geometry", 'initEdgeCyl') + CALL criticalError("Boundary type not defined in this geometry", 'initEdge2DCyl') END SELECT @@ -274,14 +272,14 @@ MODULE moduleMeshCyl !Physical surface self%physicalSurface = physicalSurface - END SUBROUTINE initEdgeCyl + END SUBROUTINE initEdge2DCyl !Random position in quadrilateral volume FUNCTION randPosVolQuad(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: xii(1:3) REAL(8), ALLOCATABLE:: fPsi(:) @@ -299,23 +297,23 @@ MODULE moduleMeshCyl END FUNCTION randposVolQuad !Get nodes from edge - PURE FUNCTION getNodesCyl(self) RESULT(n) + PURE FUNCTION getNodes2DCyl(self) RESULT(n) IMPLICIT NONE - CLASS(meshEdgeCyl), INTENT(in):: self + CLASS(meshEdge2DCyl), INTENT(in):: self INTEGER, ALLOCATABLE:: n(:) ALLOCATE(n(1:2)) n = (/self%n1%n, self%n2%n /) - END FUNCTION getNodesCyl + END FUNCTION getNodes2DCyl !Calculates a random position in edge FUNCTION randPosEdge(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshEdgeCyl), INTENT(in):: self + CLASS(meshEdge2DCyl), INTENT(in):: self REAL(8):: rnd REAL(8):: r(1:3) REAL(8):: p1(1:2), p2(1:2) @@ -332,11 +330,11 @@ MODULE moduleMeshCyl !VOLUME FUNCTIONS !QUAD FUNCTIONS !Inits quadrilateral element - SUBROUTINE initVolQuadCyl(self, n, p) + SUBROUTINE initVolQuad2DCyl(self, n, p) USE moduleRefParam IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(out):: self + CLASS(meshVol2DCylQuad), INTENT(out):: self INTEGER, INTENT(in):: n INTEGER, INTENT(in):: p(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 @@ -365,14 +363,14 @@ MODULE moduleMeshCyl CALL OMP_INIT_LOCK(self%lock) - END SUBROUTINE initVolQuadCyl + END SUBROUTINE initVolQuad2DCyl !Computes element area PURE SUBROUTINE areaQuad(self) USE moduleConstParam, ONLY: PI8 IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(inout):: self + CLASS(meshVol2DCylQuad), INTENT(inout):: self REAL(8):: r, xi(1:3) REAL(8):: detJ REAL(8):: fPsi(1:4) @@ -453,7 +451,7 @@ MODULE moduleMeshCyl PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: dPsi(1:,1:) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr @@ -469,7 +467,7 @@ MODULE moduleMeshCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), 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) @@ -502,7 +500,7 @@ MODULE moduleMeshCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) REAL(8):: r, xi(1:3) @@ -558,7 +556,7 @@ MODULE moduleMeshCyl USE moduleSpecies IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), INTENT(in):: self CLASS(particle), INTENT(in):: part TYPE(outputNode), POINTER:: vertex REAL(8):: w_p(1:4) @@ -593,7 +591,7 @@ MODULE moduleMeshCyl PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), 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 @@ -620,7 +618,7 @@ MODULE moduleMeshCyl PURE FUNCTION getNodesQuad(self) RESULT(n) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), INTENT(in):: self INTEGER, ALLOCATABLE:: n(:) ALLOCATE(n(1:4)) @@ -632,7 +630,7 @@ MODULE moduleMeshCyl PURE FUNCTION phy2logQuad(self,r) RESULT(xN) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), 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) @@ -662,7 +660,7 @@ MODULE moduleMeshCyl SUBROUTINE nextElementQuad(self, xi, nextElement) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(in):: self + CLASS(meshVol2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) CLASS(*), POINTER, INTENT(out):: nextElement REAL(8):: xiArray(1:4) @@ -687,15 +685,14 @@ MODULE moduleMeshCyl !TRIA ELEMENT !Init tria element - SUBROUTINE initVolTriaCyl(self, n, p) + SUBROUTINE initVolTria2DCyl(self, n, p) USE moduleRefParam IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(out):: self + CLASS(meshVol2DCylTria), 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 @@ -716,31 +713,18 @@ MODULE moduleMeshCyl 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) - 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 + END SUBROUTINE initVolTria2DCyl !Random position in quadrilateral volume FUNCTION randPosVolTria(self) RESULT(r) USE moduleRandom IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: xii(1:3) REAL(8), ALLOCATABLE:: fPsi(:) @@ -762,7 +746,7 @@ MODULE moduleMeshCyl USE moduleConstParam, ONLY: PI IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(inout):: self + CLASS(meshVol2DCylTria), INTENT(inout):: self REAL(8):: r, xi(1:3) REAL(8):: detJ REAL(8):: fPsi(1:3) @@ -837,7 +821,7 @@ MODULE moduleMeshCyl PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr) IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: dPsi(1:,1:) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr @@ -853,7 +837,7 @@ MODULE moduleMeshCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), 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) @@ -883,7 +867,7 @@ MODULE moduleMeshCyl USE moduleConstParam, ONLY: PI2 IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:3) @@ -938,7 +922,7 @@ MODULE moduleMeshCyl USE moduleSpecies IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self CLASS(particle), INTENT(in):: part TYPE(outputNode), POINTER:: vertex REAL(8):: w_p(1:3) @@ -968,7 +952,7 @@ MODULE moduleMeshCyl PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) REAL(8):: dPsi(1:2,1:3) REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates @@ -995,7 +979,7 @@ MODULE moduleMeshCyl PURE FUNCTION getNodesTria(self) RESULT(n) IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self INTEGER, ALLOCATABLE:: n(:) ALLOCATE(n(1:3)) @@ -1007,7 +991,7 @@ MODULE moduleMeshCyl PURE FUNCTION phy2logTria(self,r) RESULT(xi) IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: xi(1:3) REAL(8):: invJ(1:2,1:2), detJ @@ -1027,7 +1011,7 @@ MODULE moduleMeshCyl SUBROUTINE nextElementTria(self, xi, nextElement) IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(in):: self + CLASS(meshVol2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) CLASS(*), POINTER, INTENT(out):: nextElement REAL(8):: xiArray(1:3) @@ -1049,10 +1033,10 @@ MODULE moduleMeshCyl !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS !Computes element Jacobian determinant - PURE FUNCTION detJCyl(self, xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ2DCyl(self, xi, dPsi_in) RESULT(dJ) IMPLICIT NONE - CLASS(meshVolCyl), INTENT(in):: self + CLASS(meshVol2DCyl), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) REAL(8), ALLOCATABLE:: dPsi(:,:) @@ -1070,13 +1054,13 @@ MODULE moduleMeshCyl CALL self%partialDer(dPsi, dz, dr) dJ = dz(1)*dr(2)-dz(2)*dr(1) - END FUNCTION detJCyl + END FUNCTION detJ2DCyl !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJCyl(self,xi,dPsi_in) RESULT(invJ) + PURE FUNCTION invJ2DCyl(self,xi,dPsi_in) RESULT(invJ) IMPLICIT NONE - CLASS(meshVolCyl), INTENT(in):: self + CLASS(meshVol2DCyl), INTENT(in):: self REAL(8), INTENT(in):: xi(1:3) REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:) REAL(8), ALLOCATABLE:: dPsi(:,:) @@ -1095,6 +1079,6 @@ MODULE moduleMeshCyl invJ(1,:) = (/ dr(2), -dz(2) /) invJ(2,:) = (/ -dr(1), dz(1) /) - END FUNCTION invJCyl + END FUNCTION invJ2DCyl -END MODULE moduleMeshCyl +END MODULE moduleMesh2DCyl diff --git a/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylBoundary.f90 similarity index 94% rename from src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 rename to src/modules/mesh/2DCyl/moduleMesh2DCylBoundary.f90 index 5410007..7a1c581 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCylBoundary.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCylBoundary.f90 @@ -1,6 +1,6 @@ -!moduleMeshCylBoundary: Boundary functions for cylindrical coordinates -SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary - USE moduleMeshCyl +!moduleMesh2DCylBoundary: Boundary functions for cylindrical coordinates +SUBMODULE (moduleMesh2DCyl) moduleMesh2DCylBoundary + USE moduleMesh2DCyl CONTAINS SUBROUTINE reflection(edge, part) @@ -13,7 +13,7 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary !TODO: Try to do this without select SELECT TYPE(edge) - TYPE IS(meshEdgeCyl) + 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) @@ -53,7 +53,7 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary REAL(8):: d !Distance from particle to edge SELECT TYPE(edge) - TYPE IS(meshEdgeCyl) + 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) @@ -124,7 +124,7 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary !Reflects particle in the edge SELECT TYPE(edge) - TYPE IS(meshEdgeCyl) + 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) @@ -161,4 +161,4 @@ SUBMODULE (moduleMeshCyl) moduleMeshCylBoundary END SUBROUTINE symmetryAxis -END SUBMODULE moduleMeshCylBoundary +END SUBMODULE moduleMesh2DCylBoundary diff --git a/src/modules/mesh/2DCyl/moduleMeshCylRead.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 similarity index 88% rename from src/modules/mesh/2DCyl/moduleMeshCylRead.f90 rename to src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 index a3ee6db..085c3c3 100644 --- a/src/modules/mesh/2DCyl/moduleMeshCylRead.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 @@ -1,11 +1,11 @@ -MODULE moduleMeshCylRead +MODULE moduleMesh2DCylRead USE moduleMesh - USE moduleMeshCyl + USE moduleMesh2DCyl - TYPE, EXTENDS(meshGeneric):: meshCylGeneric + TYPE, EXTENDS(meshGeneric):: mesh2DCylGeneric CONTAINS - PROCEDURE, PASS:: init => initCylMesh - PROCEDURE, PASS:: readMesh => readMeshCylGmsh + PROCEDURE, PASS:: init => init2DCylMesh + PROCEDURE, PASS:: readMesh => readMesh2DCylGmsh END TYPE @@ -16,12 +16,12 @@ MODULE moduleMeshCylRead CONTAINS !Init mesh - SUBROUTINE initCylMesh(self, meshFormat) + SUBROUTINE init2DCylMesh(self, meshFormat) USE moduleMesh USE moduleErrors IMPLICIT NONE - CLASS(meshCylGeneric), INTENT(out):: self + CLASS(mesh2DCylGeneric), INTENT(out):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat SELECT CASE(meshFormat) @@ -31,18 +31,18 @@ MODULE moduleMeshCylRead self%printEM => printEMGmsh CASE DEFAULT - CALL criticalError("Mesh type " // meshFormat // " not supported.", "initCylMesh") + CALL criticalError("Mesh type " // meshFormat // " not supported.", "init2DCylMesh") END SELECT - END SUBROUTINE initCylMesh + END SUBROUTINE init2DCylMesh !Read mesh from gmsh file - SUBROUTINE readMeshCylGmsh(self, filename) + SUBROUTINE readMesh2DCylGmsh(self, filename) USE moduleBoundary IMPLICIT NONE - CLASS(meshCylGeneric), INTENT(inout):: self + CLASS(mesh2DCylGeneric), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename REAL(8):: r, z INTEGER:: p(1:4) @@ -68,8 +68,8 @@ MODULE moduleMeshCylRead !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 /)) + ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj) + CALL self%nodes(n)%obj%init(n, (/z, r, 0.D0 /)) END DO !Skips comments @@ -103,7 +103,7 @@ MODULE moduleMeshCylRead !Associate boundary condition procedure. bt = getBoundaryId(boundaryType) - ALLOCATE(meshEdgeCyl:: self%edges(e)%obj) + ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj) CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType) @@ -118,13 +118,13 @@ MODULE moduleMeshCylRead CASE (2) !Triangular element READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3) - ALLOCATE(meshVolCylTria:: self%vols(e)%obj) + ALLOCATE(meshVol2DCylTria:: 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) + ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj) CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4)) END SELECT @@ -154,7 +154,7 @@ MODULE moduleMeshCylRead END DO - END SUBROUTINE readMeshCylGmsh + END SUBROUTINE readMesh2DCylGmsh !Selects type of elements to build connection SUBROUTINE connectedVolVol(elemA, elemB) @@ -164,27 +164,27 @@ MODULE moduleMeshCylRead CLASS(meshVol), INTENT(inout):: elemB SELECT TYPE(elemA) - TYPE IS(meshVolCylQuad) + TYPE IS(meshVol2DCylQuad) !Element A is a quadrilateral SELECT TYPE(elemB) - TYPE IS(meshVolCylQuad) + TYPE IS(meshVol2DCylQuad) !Element B is a quadrilateral CALL connectedQuadQuad(elemA, elemB) - TYPE IS(meshVolCylTria) + TYPE IS(meshVol2DCylTria) !Element B is a triangle CALL connectedQuadTria(elemA, elemB) END SELECT - TYPE IS(meshVolCylTria) + TYPE IS(meshVol2DCylTria) !Element A is a Triangle SELECT TYPE(elemB) - TYPE IS(meshVolCylQuad) + TYPE IS(meshVol2DCylQuad) !Element B is a quadrilateral CALL connectedQuadTria(elemB, elemA) - TYPE IS(meshVolCylTria) + TYPE IS(meshVol2DCylTria) !Element B is a triangle CALL connectedTriaTria(elemA, elemB) @@ -201,13 +201,13 @@ MODULE moduleMeshCylRead CLASS(meshEdge), INTENT(inout):: elemB SELECT TYPE(elemB) - CLASS IS(meshEdgeCyl) + CLASS IS(meshEdge2DCyl) SELECT TYPE(elemA) - TYPE IS(meshVolCylQuad) + TYPE IS(meshVol2DCylQuad) !Element A is a quadrilateral CALL connectedQuadEdge(elemA, elemB) - TYPE IS(meshVolCylTria) + TYPE IS(meshVol2DCylTria) !Element A is a triangle CALL connectedTriaEdge(elemA, elemB) @@ -220,8 +220,8 @@ MODULE moduleMeshCylRead SUBROUTINE connectedQuadQuad(elemA, elemB) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemB + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1) .AND. & @@ -264,8 +264,8 @@ MODULE moduleMeshCylRead SUBROUTINE connectedQuadTria(elemA, elemB) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -356,8 +356,8 @@ MODULE moduleMeshCylRead SUBROUTINE connectedTriaTria(elemA, elemB) IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemB + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -429,8 +429,8 @@ MODULE moduleMeshCylRead SUBROUTINE connectedQuadEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA - CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB + CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -501,8 +501,8 @@ MODULE moduleMeshCylRead SUBROUTINE connectedTriaEdge(elemA, elemB) IMPLICIT NONE - CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA - CLASS(meshEdgeCyl), INTENT(inout), TARGET:: elemB + CLASS(meshVol2DCylTria), INTENT(inout), TARGET:: elemA + CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB !Check direction 1 IF (.NOT. ASSOCIATED(elemA%e1)) THEN @@ -564,7 +564,7 @@ MODULE moduleMeshCylRead INTEGER, ALLOCATABLE:: n(:) SELECT TYPE(elem) - TYPE IS(meshVolCylQuad) + TYPE IS(meshVol2DCylQuad) nNodes = 4 ALLOCATE(localK(1:nNodes,1:nNodes)) localK = elem%elemK() @@ -572,7 +572,7 @@ MODULE moduleMeshCylRead n = (/ elem%n1%n, elem%n2%n, & elem%n3%n, elem%n4%n /) - TYPE IS(meshVolCylTria) + TYPE IS(meshVol2DCylTria) nNodes = 3 ALLOCATE(localK(1:nNodes,1:nNodes)) localK = elem%elemK() @@ -596,4 +596,4 @@ MODULE moduleMeshCylRead END SUBROUTINE constructGlobalK -END MODULE moduleMeshCylRead +END MODULE moduleMesh2DCylRead diff --git a/src/modules/mesh/makefile b/src/modules/mesh/makefile index 741a43e..0eeaa23 100644 --- a/src/modules/mesh/makefile +++ b/src/modules/mesh/makefile @@ -1,8 +1,11 @@ -all: moduleMesh.o 2DCyl.o 1DRad.o 1DCart.o +all: moduleMesh.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 2DCyl.o: $(MAKE) -C 2DCyl all +2DCart.o: + $(MAKE) -C 2DCart all + 1DCart.o: $(MAKE) -C 1DCart all diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index 2a6e334..a96bbc5 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -365,7 +365,6 @@ MODULE moduleInput !Gets the basename of the folder CALL config%get(object // '.folder', baseName, found) - PRINT *, baseName IF (found) THEN folder = baseName @@ -654,7 +653,8 @@ MODULE moduleInput !Read the geometry (mesh) for the case SUBROUTINE readGeometry(config) USE moduleMesh - USE moduleMeshCylRead, ONLY: meshCylGeneric + USE moduleMesh2DCylRead, ONLY: mesh2DCylGeneric + USE moduleMesh2DCartRead, ONLY: mesh2DCartGeneric USE moduleMesh1DCartRead, ONLY: mesh1DCartGeneric USE moduleMesh1DRadRead, ONLY: mesh1DRadGeneric USE moduleErrors @@ -672,7 +672,11 @@ MODULE moduleInput SELECT CASE(geometryType) CASE ("2DCyl") !Creates a 2D cylindrical mesh - ALLOCATE(meshCylGeneric:: mesh) + ALLOCATE(mesh2DCylGeneric:: mesh) + + CASE ("2DCart") + !Creates a 2D cylindrical mesh + ALLOCATE(mesh2DCartGeneric:: mesh) CASE ("1DCart") !Creates a 1D cartesian mesh diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index dfd5d4a..4b0124e 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -67,15 +67,31 @@ MODULE moduleSolver REAL(8):: tau, tauSp SELECT CASE(pusherType) + !2D Cylindrical CASE('2DCylNeutral') - self%pushParticle => pushCylNeutral + self%pushParticle => push2DCylNeutral CASE('2DCylCharged') - self%pushParticle => pushCylCharged + self%pushParticle => push2DCylCharged + + !2D Cartesian + CASE('2DCartNeutral') + self%pushParticle => push2DCartNeutral + + CASE('2DCartCharged') + self%pushParticle => push2DCartCharged + + !1D Cartesian + CASE('1DCartNeutral') + self%pushParticle => push1DCartNeutral CASE('1DCartCharged') self%pushParticle => push1DCartCharged + !1D Radial + CASE('1DRadNeutral') + self%pushParticle => push1DRadNeutral + CASE('1DRadCharged') self%pushParticle => push1DRadCharged @@ -147,7 +163,7 @@ MODULE moduleSolver END SUBROUTINE doPushes !Push one particle. Boris pusher for 2D Cyl Neutral particle - PURE SUBROUTINE pushCylNeutral(part, tauIn) + PURE SUBROUTINE push2DCylNeutral(part, tauIn) USE moduleSpecies IMPLICIT NONE @@ -181,10 +197,10 @@ MODULE moduleSolver !Copy temporal particle to particle part=part_temp - END SUBROUTINE pushCylNeutral + END SUBROUTINE push2DCylNeutral !Push one particle. Boris pusher for 2D Cyl Charged particle - PURE SUBROUTINE pushCylCharged(part, tauIn) + PURE SUBROUTINE push2DCylCharged(part, tauIn) USE moduleSpecies USE moduleEM IMPLICIT NONE @@ -222,7 +238,83 @@ MODULE moduleSolver !Copy temporal particle to particle part=part_temp - END SUBROUTINE pushCylCharged + END SUBROUTINE push2DCylCharged + + !Push neutral particles in 2D cartesian coordinates + PURE SUBROUTINE push2DCartNeutral(part, tauIn) + USE moduleSPecies + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn + TYPE(particle):: part_temp + + part_temp = part + + !x + part_temp%v(1) = part%v(1) + part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn + + !y + part_temp%v(2) = part%v(2) + part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn + + part_temp%n_in = .FALSE. + + part = part_temp + + END SUBROUTINE push2DCartNeutral + + !Push charged particles in 2D cartesian coordinates + PURE SUBROUTINE push2DCartCharged(part, tauIn) + USE moduleSPecies + USE moduleEM + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn + TYPE(particle):: part_temp + REAL(8):: qmEFt(1:3) + + part_temp = part + !Get the electric field at particle position + qmEFt = part_temp%qm*gatherElecField(part_temp)*tauIn + + !x + part_temp%v(1) = part%v(1) + qmEFt(1) + part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn + + !y + part_temp%v(2) = part%v(2) + qmEFt(2) + part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn + + part_temp%n_in = .FALSE. + + part = part_temp + + END SUBROUTINE push2DCartCharged + + !Push neutral particles in 1D cartesian coordinates + PURE SUBROUTINE push1DCartNeutral(part, tauIn) + USE moduleSPecies + USE moduleEM + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn + TYPE(particle):: part_temp + + part_temp = part + + !x + part_temp%v(1) = part%v(1) + part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn + + part_temp%n_in = .FALSE. + + part = part_temp + + END SUBROUTINE push1DCartNeutral !Push charged particles in 1D cartesian coordinates PURE SUBROUTINE push1DCartCharged(part, tauIn) @@ -249,6 +341,41 @@ MODULE moduleSolver END SUBROUTINE push1DCartCharged + !Push one particle. Boris pusher for 1D Radial Neutral particle + PURE SUBROUTINE push1DRadNeutral(part, tauIn) + USE moduleSpecies + USE moduleEM + IMPLICIT NONE + + TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn + REAL(8):: v_p_oh_star(1:2) + TYPE(particle):: part_temp + REAL(8):: x_new, y_new, r, sin_alpha, cos_alpha + + part_temp = part + !r,theta + v_p_oh_star(1) = part%v(1) + x_new = part%r(1) + v_p_oh_star(1)*tauIn + v_p_oh_star(2) = part%v(2) + y_new = v_p_oh_star(2)*tauIn + r = DSQRT(x_new**2+y_new**2) + part_temp%r(1) = 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(1) = cos_alpha*v_p_oh_star(1)+sin_alpha*v_p_oh_star(2) + part_temp%v(2) = -sin_alpha*v_p_oh_star(1)+cos_alpha*v_p_oh_star(2) + part_temp%n_in = .FALSE. !Assume particle is outside until cell is found + !Copy temporal particle to particle + part=part_temp + + END SUBROUTINE push1DRadNeutral + !Push one particle. Boris pusher for 1D Radial Charged particle PURE SUBROUTINE push1DRadCharged(part, tauIn) USE moduleSpecies @@ -263,7 +390,6 @@ MODULE moduleSolver REAL(8):: qmEFt(1:3)!charge*tauIn*EF/mass part_temp = part - !Time step for the species !Get electric field at particle position qmEFt = part_temp%qm*gatherElecField(part_temp)*tauMin !r,theta