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.
This commit is contained in:
Jorge Gonzalez 2023-01-07 10:47:18 +01:00
commit 1c5b887a6d
4 changed files with 25 additions and 40 deletions

View file

@ -34,7 +34,7 @@
{"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": {

View file

@ -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

View file

@ -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

View file

@ -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