From 35bd61fda9088b235b6faaedcca61a3b34fbee05 Mon Sep 17 00:00:00 2001 From: Jorge Gonzalez Date: Wed, 21 Apr 2021 23:40:58 +0200 Subject: [PATCH] Common scatter subroutine All subroutines of scattering particle properties to the nodes of a volume have been converged into one in moduleMesh. --- src/modules/mesh/0D/moduleMesh0D.f90 | 25 ---- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 45 -------- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 45 -------- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 115 ------------------- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 115 ------------------- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 58 ---------- src/modules/mesh/moduleMesh.f90 | 36 +++++- src/modules/moduleInject.f90 | 2 +- src/modules/moduleInput.f90 | 6 +- 9 files changed, 39 insertions(+), 408 deletions(-) diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index f51efd2..79939cd 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -18,7 +18,6 @@ MODULE moduleMesh0D PROCEDURE, PASS:: getNodes => getNodes0D PROCEDURE, PASS:: randPos => randPos0D PROCEDURE, NOPASS:: fPsi => fPsi0D - PROCEDURE, PASS:: scatter => scatter0D PROCEDURE, PASS:: gatherEF => gatherEF0D PROCEDURE, PASS:: elemK => elemK0D PROCEDURE, PASS:: elemF => elemF0D @@ -113,30 +112,6 @@ MODULE moduleMesh0D END FUNCTION fPsi0D - SUBROUTINE scatter0D(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol0D), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: tensorS(1:3,1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - tensorS = outerProduct(part%v, part%v) - - node => self%n1 - sp = part%species%n - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatter0D - PURE FUNCTION gatherEF0D(self, xi) RESULT(EF) IMPLICIT NONE diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index e96150c..24b7287 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -75,9 +75,7 @@ MODULE moduleMesh1DCart PROCEDURE, PASS:: partialDer => partialDerSegm PROCEDURE, PASS:: elemK => elemKSegm PROCEDURE, PASS:: elemF => elemFSegm - PROCEDURE, NOPASS:: weight => weightSegm PROCEDURE, NOPASS:: inside => insideSegm - PROCEDURE, PASS:: scatter => scatterSegm PROCEDURE, PASS:: gatherEF => gatherEFSegm PROCEDURE, PASS:: getNodes => getNodesSegm PROCEDURE, PASS:: phy2log => phy2logSegm @@ -357,16 +355,6 @@ MODULE moduleMesh1DCart END FUNCTION elemFSegm - PURE FUNCTION weightSegm(xi) RESULT(w) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: w(1:2) - - w = fPsiSegm(xi) - - END FUNCTION weightSegm - PURE FUNCTION insideSegm(xi) RESULT(ins) IMPLICIT NONE @@ -378,39 +366,6 @@ MODULE moduleMesh1DCart END FUNCTION insideSegm - SUBROUTINE scatterSegm(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol1DCartSegm), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: w_p(1:2) - REAL(8):: tensorS(1:3,1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - node => self%n1 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n2 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatterSegm - !Gathers EF at position Xii PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF) IMPLICIT NONE diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 99b3331..7567ab9 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -76,9 +76,7 @@ MODULE moduleMesh1DRad PROCEDURE, PASS:: partialDer => partialDerRad PROCEDURE, PASS:: elemK => elemKRad PROCEDURE, PASS:: elemF => elemFRad - PROCEDURE, NOPASS:: weight => weightRad PROCEDURE, NOPASS:: inside => insideRad - PROCEDURE, PASS:: scatter => scatterRad PROCEDURE, PASS:: gatherEF => gatherEFRad PROCEDURE, PASS:: getNodes => getNodesRad PROCEDURE, PASS:: phy2log => phy2logRad @@ -369,16 +367,6 @@ MODULE moduleMesh1DRad END FUNCTION elemFRad - PURE FUNCTION weightRad(xi) RESULT(w) - IMPLICIT NONE - - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: w(1:2) - - w = fPsiRad(xi) - - END FUNCTION weightRad - PURE FUNCTION insideRad(xi) RESULT(ins) IMPLICIT NONE @@ -390,39 +378,6 @@ MODULE moduleMesh1DRad END FUNCTION insideRad - SUBROUTINE scatterRad(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol1DRadSegm), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: w_p(1:2) - REAL(8):: tensorS(1:3,1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - node => self%n1 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n2 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(2) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(2)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatterRad - !Gathers EF at position Xii PURE FUNCTION gatherEFRad(self, xi) RESULT(EF) IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 36a182a..ff17688 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -83,9 +83,7 @@ MODULE moduleMesh2DCart 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 @@ -114,9 +112,7 @@ MODULE moduleMesh2DCart 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 @@ -470,17 +466,6 @@ MODULE moduleMesh2DCart 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 @@ -493,54 +478,6 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Scatter properties of particle into element nodes - SUBROUTINE scatterQuad(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol2DCartQuad), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: w_p(1:4) - REAL(8):: tensorS(1:3,1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - node => self%n1 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n2 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(2) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(2)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n3 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(3) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(3)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n4 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(4) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(4)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(4)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatterQuad - !Gathers the electric field at position xi PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) IMPLICIT NONE @@ -841,17 +778,6 @@ MODULE moduleMesh2DCart 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 @@ -864,47 +790,6 @@ MODULE moduleMesh2DCart END FUNCTION insideTria - !Scatter properties of particles into element - SUBROUTINE scatterTria(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol2DCartTria), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: w_p(1:3) - REAL(8):: tensorS(1:3,1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - node => self%n1 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n2 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(2) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(2)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n3 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(3) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(3)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatterTria - !Gathers the electric field at position xi PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) IMPLICIT NONE diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index d4e06bd..dcd285d 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -84,9 +84,7 @@ MODULE moduleMesh2DCyl 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 @@ -115,9 +113,7 @@ MODULE moduleMesh2DCyl 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 @@ -491,17 +487,6 @@ MODULE moduleMesh2DCyl 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 @@ -514,54 +499,6 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Scatter properties of particle into element nodes - SUBROUTINE scatterQuad(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol2DCylQuad), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: w_p(1:4) - REAL(8):: tensorS(1:3,1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - node => self%n1 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n2 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(2) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(2)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n3 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(3) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(3)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n4 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(4) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(4)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(4)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatterQuad - !Gathers the electric field at position xi PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) IMPLICIT NONE @@ -871,17 +808,6 @@ MODULE moduleMesh2DCyl 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 @@ -894,47 +820,6 @@ MODULE moduleMesh2DCyl END FUNCTION insideTria - !Scatter properties of particles into element - SUBROUTINE scatterTria(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol2DCylTria), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: w_p(1:3) - REAL(8):: tensorS(1:3,1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - node => self%n1 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n2 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(2) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(2)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n3 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(3) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(3)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatterTria - !Gathers the electric field at position xi PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) IMPLICIT NONE diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 8c6e350..15aed4d 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -76,9 +76,7 @@ MODULE moduleMesh3DCart 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 @@ -459,15 +457,6 @@ MODULE moduleMesh3DCart 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 @@ -481,53 +470,6 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra - SUBROUTINE scatterTetra(self, part) - USE moduleMath - USE moduleSpecies - USE OMP_LIB - IMPLICIT NONE - - CLASS(meshVol3DCartTetra), INTENT(in):: self - CLASS(particle), INTENT(in):: part - REAL(8):: w_p(1:4) - REAL(8):: tensorS(1:3, 1:3) - CLASS(meshNode), POINTER:: node - INTEGER:: sp - - w_p = self%weight(part%xi) - tensorS = outerProduct(part%v, part%v) - - sp = part%species%n - node => self%n1 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(1) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(1)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(1)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n2 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(2) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(2)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(2)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n3 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(3) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(3)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(3)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - node => self%n4 - CALL OMP_SET_LOCK(node%lock) - node%output(sp)%den = node%output(sp)%den + part%weight*w_p(4) - node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*w_p(4)*part%v(:) - node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*w_p(4)*tensorS - CALL OMP_UNSET_LOCK(node%lock) - - END SUBROUTINE scatterTetra - PURE FUNCTION gatherEFTetra(self, xi) RESULT(EF) IMPLICIT NONE diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 1b1f70a..030f725 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -164,7 +164,7 @@ MODULE moduleMesh PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi - PROCEDURE(scatter_interface), DEFERRED, PASS:: scatter + PROCEDURE, PASS:: scatter PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF @@ -468,6 +468,40 @@ MODULE moduleMesh END SUBROUTINE resetOutput + !Scatters particle properties into vol nodes + SUBROUTINE scatter(self, part) + USE moduleMath + USE moduleSpecies + USE OMP_LIB + IMPLICIT NONE + + CLASS(meshVol), INTENT(inout):: self + CLASS(particle), INTENT(in):: part + REAL(8), ALLOCATABLE:: fPsi(:) + INTEGER, ALLOCATABLE:: volNodes(:) + REAL(8):: tensorS(1:3, 1:3) + INTEGER:: sp + INTEGER:: i, nNodes + CLASS(meshNode), POINTER:: node + + fPsi = self%fPsi(part%xi) + tensorS = outerProduct(part%v, part%v) + sp = part%species%n + volNodes = self%getNodes() + nNodes = SIZE(volNodes) + + DO i = 1, nNodes + node => mesh%nodes(volNodes(i))%obj + CALL OMP_SET_LOCK(node%lock) + node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i) + node%output(sp)%mom(:) = node%output(sp)%mom(:) + part%weight*fPsi(i)*part%v(:) + node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS + CALL OMP_UNSET_LOCK(node%lock) + + END DO + + END SUBROUTINE scatter + !Find next cell for particle RECURSIVE SUBROUTINE findCell(self, part, oldCell) USE moduleSpecies diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90 index 8c8f65d..9bfb35f 100644 --- a/src/modules/moduleInject.f90 +++ b/src/modules/moduleInject.f90 @@ -127,7 +127,7 @@ MODULE moduleInject et = et + 1 self%edges(et) = mesh%edges(e)%obj%n !Assign connectivity between injection edge and meshColl volume - IF (ASSOCIATED(meshForMCC, meshColl)) THEN + IF (doubleMesh) THEN nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos()) IF (nVolColl > 0) THEN mesh%edges(e)%obj%eColl => meshColl%vols(nVolColl)%obj diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index b420861..34e7507 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -345,11 +345,11 @@ MODULE moduleInput partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian() partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian() partNew%vol = e - IF (ASSOCIATED(meshForMCC, mesh)) THEN - partNew%volColl = partNew%vol + IF (doubleMesh) THEN + partNew%volColl = findCellBrute(meshColl, partNew%r) ELSE - partNew%volColl = findCellBrute(meshColl, partNew%r) + partNew%volColl = partNew%vol END IF partNew%n_in = .TRUE.