From 1c5b887a6d7b0221efd7fcceb072330406f0d7f4 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sat, 7 Jan 2023 10:47:18 +0100 Subject: [PATCH] Small bugfix when testing examples While testing the examples distributed with the code, a few errors were found and fixed, mostly related with the K matrix in 1D geometry and reading values from initial conditions for species. --- runs/1D_Cathode/inputRad.json | 4 +- src/modules/init/moduleInput.f90 | 47 +++++++------------- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 6 +-- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 8 ++-- 4 files changed, 25 insertions(+), 40 deletions(-) diff --git a/runs/1D_Cathode/inputRad.json b/runs/1D_Cathode/inputRad.json index 34a09ce..824c0a8 100644 --- a/runs/1D_Cathode/inputRad.json +++ b/runs/1D_Cathode/inputRad.json @@ -31,10 +31,10 @@ ]} ], "boundaryEM": [ - {"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1} + {"name": "Cathode", "type": "dirichlet", "potential": 0.0, "physicalSurface": 1} ], "inject": [ - {"name": "Plasma Cat e", "species": "Electron", "flow": 2.64e-2, "units": "A", "v": 180000.0, "T": [ 2300.0, 2300.0, 2300.0], + {"name": "Plasma Cat e", "species": "Electron", "flow": 2.64e-5, "units": "A", "v": 180000.0, "T": [ 2300.0, 2300.0, 2300.0], "velDist": ["Maxwellian", "Maxwellian", "Maxwellian"], "n": [ 1, 0, 0], "physicalSurface": 1} ], "solver": { diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index ff2cc4b..60f33ec 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -331,7 +331,8 @@ MODULE moduleInput REAL(8), ALLOCATABLE, DIMENSION(:,:):: velocity INTEGER, ALLOCATABLE, DIMENSION(:):: nodes INTEGER:: nNodes - REAL(8), ALLOCATABLE, DIMENSION(:):: source, fPsi + REAL(8), ALLOCATABLE, DIMENSION(:):: sourceScalar + REAL(8), ALLOCATABLE, DIMENSION(:,:):: sourceArray !Density at the volume centroid REAL(8):: densityCen !Mean velocity and temperature at particle position @@ -361,14 +362,10 @@ MODULE moduleInput !Density at centroid of cell nNodes = mesh%cells(e)%obj%nNodes nodes = mesh%cells(e)%obj%getNodes(nNodes) - ALLOCATE(fPsi(1:nNodes)) - fPsi = mesh%cells(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/), nNodes) - ALLOCATE(source(1:nNodes)) - DO j = 1, nNodes - source(j) = density(nodes(j)) - - END DO - densityCen = DOT_PRODUCT(fPsi, source) + ALLOCATE(sourceScalar(1:nNodes)) + ALLOCATE(sourceArray(1:nNodes, 1:3)) + sourceScalar = density(nodes) + densityCen = mesh%cells(e)%obj%gatherF((/ 0.D0, 0.D0, 0.D0 /), nNodes, sourceScalar) !Calculate number of particles nNewPart = INT(densityCen * (mesh%cells(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) @@ -380,37 +377,23 @@ MODULE moduleInput partNew%r = mesh%cells(e)%obj%randPos() partNew%Xi = mesh%cells(e)%obj%phy2log(partNew%r) !Get mean velocity at particle position - fPsi = mesh%cells(e)%obj%fPsi(partNew%Xi, nNodes) - DO j = 1, nNodes - source(j) = velocity(nodes(j), 1) - - END DO - velocityXi(1) = DOT_PRODUCT(fPsi, source) - DO j = 1, nNodes - source(j) = velocity(nodes(j), 2) - - END DO - velocityXi(2) = DOT_PRODUCT(fPsi, source) - DO j = 1, nNodes - source(j) = velocity(nodes(j), 3) - - END DO - velocityXi(3) = DOT_PRODUCT(fPsi, source) + sourceArray(:,1) = velocity((nodes), 1) + sourceArray(:,2) = velocity((nodes), 2) + sourceArray(:,3) = velocity((nodes), 3) + velocityXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceArray) velocityXi = velocityXi / v_ref !Get temperature at particle position - DO j = 1, nNodes - source(j) = temperature(nodes(j)) - - END DO - temperatureXi = DOT_PRODUCT(fPsi, source) + sourceScalar = temperature(nodes) + temperatureXi = mesh%cells(e)%obj%gatherF(partNew%Xi, nNodes, sourceScalar) temperatureXi = temperatureXi / T_ref vTh = DSQRT(temperatureXi / species(sp)%obj%m) partNew%v(1) = velocityXi(1) + vTh*randomMaxwellian() partNew%v(2) = velocityXi(2) + vTh*randomMaxwellian() partNew%v(3) = velocityXi(3) + vTh*randomMaxwellian() + partNew%vol = e IF (doubleMesh) THEN partNew%volColl = findCellBrute(meshColl, partNew%r) @@ -419,7 +402,9 @@ MODULE moduleInput partNew%volColl = partNew%vol END IF + partNew%n_in = .TRUE. + partNew%weight = species(sp)%obj%weight !Assign particle to temporal list of particles @@ -434,7 +419,7 @@ MODULE moduleInput END DO - DEALLOCATE(source) + DEALLOCATE(sourceScalar, sourceArray) END DO diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 82f43c2..2f8b7b8 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -324,7 +324,7 @@ MODULE moduleMesh1DCart REAL(8):: Xi(1:3) REAL(8):: dPsi(1:3, 1:2) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3,1:3), detJ + REAL(8):: invJ(1:3, 1:3), detJ INTEGER:: l localK = 0.D0 @@ -337,8 +337,8 @@ MODULE moduleMesh1DCart pDer = self%partialDer(2, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) - localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & wSeg(l)/detJ END DO diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index c8ff414..81a8864 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -326,7 +326,7 @@ MODULE moduleMesh1DRad REAL(8):: dPsi(1:3, 1:2) REAL(8):: pDer(1:3, 1:3) REAL(8):: r - REAL(8):: invJ(1:3,1:3), detJ + REAL(8):: invJ(1:3, 1:3), detJ INTEGER:: l localK = 0.D0 @@ -339,9 +339,9 @@ MODULE moduleMesh1DRad pDer = self%partialDer(2, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) - r = self%gatherF(Xi, 4, self%r) - localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + r = self%gatherF(Xi, 2, self%r) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & r*wSeg(l)/detJ END DO