From dd1fca3fee5dad0bf31e7a93ef333cd5da335aa8 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Tue, 20 Dec 2022 15:51:43 +0100 Subject: [PATCH] Fix for node volume An issue in the node volume calculation in cylindrical coordinates was found. This was causing wrong conservation of current. Still to test with ALPHIE_Grid case. Still to check triangular element. Still to theck 1D radial geometry --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 15 +++++++++++++-- .../mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 | 10 ++++++---- 2 files changed, 19 insertions(+), 6 deletions(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 0a6b7f4..0ac01f3 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -249,7 +249,7 @@ MODULE moduleMesh2DCyl dr = self%r(2) - self%r(1) dz = self%z(2) - self%z(1) IF (dr /= 0.D0) THEN - r(2) = dr*DSQRT(rnd) + self%r(1) + r(2) = dr * DSQRT(rnd) + self%r(1) r(1) = dz * (r(2) - self%r(1))/dr + self%z(1) ELSE @@ -320,7 +320,18 @@ MODULE moduleMesh2DCyl fPsi = self%fPsi(xi) r = DOT_PRODUCT(fPsi,self%r) self%volume = r*detJ - self%arNodes = fPsi*r*detJ + xi = (/-5.D-1, -5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(1) = fPsi(1)*r*detJ + xi = (/ 5.D-1, -5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(2) = fPsi(2)*r*detJ + xi = (/ 5.D-1, 5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(3) = fPsi(3)*r*detJ + xi = (/-5.D-1, 5.D-1, 0.D0/) + r = DOT_PRODUCT(self%fPsi(xi),self%r) + self%arNodes(4) = fPsi(4)*r*detJ END SUBROUTINE areaQuad diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 2dcb7f6..7eade3e 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -151,8 +151,9 @@ MODULE moduleMeshOutputGmsh2 INTEGER:: n REAL(8):: time CHARACTER(:), ALLOCATABLE:: fileName - CHARACTER (LEN=iterationDigits):: tstring + CHARACTER (LEN=iterationDigits):: tString CHARACTER(:), ALLOCATABLE:: title + CHARACTER (LEN=2):: cString SELECT TYPE(self) TYPE IS(meshParticles) @@ -168,9 +169,9 @@ MODULE moduleMeshOutputGmsh2 IF (collOutput) THEN time = DBLE(t)*tauMin*ti_ref - WRITE(tstring, iterationFormat) t + WRITE(tString, iterationFormat) t - fileName='OUTPUT_' // tstring// '_Collisions.msh' + fileName='OUTPUT_' // tString// '_Collisions.msh' WRITE(*, "(6X,A15,A)") "Creating file: ", fileName OPEN (60, file = path // folder // '/' // fileName) @@ -178,7 +179,8 @@ MODULE moduleMeshOutputGmsh2 DO k = 1, nCollPairs DO c = 1, interactionMatrix(k)%amount - WRITE(title, "(5A,I2)") '"Pair ', interactionMatrix(k)%sp_i%name, '-', interactionMatrix(k)%sp_j%name, ' collision ', c + WRITE(cString, "(I2)") c + title = '"Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString CALL writeGmsh2HeaderElementData(60, title, t, time, 1, self%numVols) DO n=1, self%numVols WRITE(60, "(I6,I10)") n + numEdges, self%vols(n)%obj%tallyColl(k)%tally(c)