diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index 79939cd..93fd4bd 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -19,6 +19,7 @@ MODULE moduleMesh0D PROCEDURE, PASS:: randPos => randPos0D PROCEDURE, NOPASS:: fPsi => fPsi0D PROCEDURE, PASS:: gatherEF => gatherEF0D + PROCEDURE, PASS:: gatherMF => gatherMF0D PROCEDURE, PASS:: elemK => elemK0D PROCEDURE, PASS:: elemF => elemF0D PROCEDURE, PASS:: phy2log => phy2log0D @@ -123,6 +124,17 @@ MODULE moduleMesh0D END FUNCTION gatherEF0D + PURE FUNCTION gatherMF0D(self, xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol0D), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: MF(1:3) + + MF = 0.D0 + + END FUNCTION gatherMF0D + PURE FUNCTION elemK0D(self) RESULT(localK) IMPLICIT NONE diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 24b7287..fe755fa 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -77,6 +77,7 @@ MODULE moduleMesh1DCart PROCEDURE, PASS:: elemF => elemFSegm PROCEDURE, NOPASS:: inside => insideSegm PROCEDURE, PASS:: gatherEF => gatherEFSegm + PROCEDURE, PASS:: gatherMF => gatherMFSegm PROCEDURE, PASS:: getNodes => getNodesSegm PROCEDURE, PASS:: phy2log => phy2logSegm PROCEDURE, PASS:: nextElement => nextElementSegm @@ -388,6 +389,28 @@ MODULE moduleMesh1DCart END FUNCTION gatherEFSegm + PURE FUNCTION gatherMFSegm(self, xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol1DCartSegm), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: fPsi(1:2) + REAL(8):: MF_Nodes(1:2, 1:3) + REAL(8):: MF(1:3) + REAL(8):: invJ + + MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1) /) + MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2) /) + MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3) /) + + fPsi = self%fPsi(xi) + MF = MATMUL(fPsi, MF_Nodes) + + END FUNCTION gatherMFSegm + !Get nodes from 1D volume PURE FUNCTION getNodesSegm(self) RESULT(n) IMPLICIT NONE diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 7567ab9..6fc3c66 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -78,6 +78,7 @@ MODULE moduleMesh1DRad PROCEDURE, PASS:: elemF => elemFRad PROCEDURE, NOPASS:: inside => insideRad PROCEDURE, PASS:: gatherEF => gatherEFRad + PROCEDURE, PASS:: gatherMF => gatherMFRad PROCEDURE, PASS:: getNodes => getNodesRad PROCEDURE, PASS:: phy2log => phy2logRad PROCEDURE, PASS:: nextElement => nextElementRad @@ -400,6 +401,28 @@ MODULE moduleMesh1DRad END FUNCTION gatherEFRad + PURE FUNCTION gatherMFRad(self, xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol1DRadSegm), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: fPsi(1:2) + REAL(8):: MF_Nodes(1:2, 1:3) + REAL(8):: MF(1:3) + REAL(8):: invJ + + MF_Nodes(1:2,1) = (/ self%n1%emData%B(1), & + self%n2%emData%B(1) /) + MF_Nodes(1:2,2) = (/ self%n1%emData%B(2), & + self%n2%emData%B(2) /) + MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), & + self%n2%emData%B(3) /) + + fPsi = self%fPsi(xi) + MF = MATMUL(fPsi, MF_Nodes) + + END FUNCTION gatherMFRad + !Get nodes from 1D volume PURE FUNCTION getNodesRad(self) RESULT(n) IMPLICIT NONE diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index bc60a69..4b23ff1 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -85,6 +85,7 @@ MODULE moduleMesh2DCart PROCEDURE, PASS:: elemF => elemFQuad PROCEDURE, NOPASS:: inside => insideQuad PROCEDURE, PASS:: gatherEF => gatherEFQuad + PROCEDURE, PASS:: gatherMF => gatherMFQuad PROCEDURE, PASS:: getNodes => getNodesQuad PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: nextElement => nextElementQuad @@ -114,6 +115,7 @@ MODULE moduleMesh2DCart PROCEDURE, PASS:: elemF => elemFTria PROCEDURE, NOPASS:: inside => insideTria PROCEDURE, PASS:: gatherEF => gatherEFTria + PROCEDURE, PASS:: gatherMF => gatherMFTria PROCEDURE, PASS:: getNodes => getNodesTria PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: nextElement => nextElementTria @@ -505,6 +507,33 @@ MODULE moduleMesh2DCart END FUNCTION gatherEFQuad + PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol2DCartQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: fPsi(1:4) + REAL(8):: MF_Nodes(1:4,1:3) + REAL(8):: MF(1:3) + + MF_Nodes(1:4,1) = (/self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1), & + self%n4%emData%B(1) /) + MF_Nodes(1:4,2) = (/self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2), & + self%n4%emData%B(2) /) + MF_Nodes(1:4,3) = (/self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3), & + self%n4%emData%B(3) /) + + fPsi = self%fPsi(xi) + MF = MATMUL(fPsi(:), MF_Nodes) + + END FUNCTION gatherMFQuad + !Gets nodes from quadrilateral element PURE FUNCTION getNodesQuad(self) RESULT(n) IMPLICIT NONE @@ -816,6 +845,30 @@ MODULE moduleMesh2DCart END FUNCTION gatherEFTria + PURE FUNCTION gatherMFTria(self,xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol2DCartTria), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: fPsi(1:3) + REAL(8):: MF_Nodes(1:3,1:3) + REAL(8):: MF(1:3) + + MF_Nodes(1:3,1) = (/self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1) /) + MF_Nodes(1:3,2) = (/self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2) /) + MF_Nodes(1:3,3) = (/self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3) /) + + fPsi = self%fPsi(xi) + MF = MATMUL(fPsi, MF_Nodes) + + END FUNCTION gatherMFTria + !Gets node indexes from triangular element PURE FUNCTION getNodesTria(self) RESULT(n) IMPLICIT NONE diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 62933d8..1e4e580 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -86,6 +86,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: elemF => elemFQuad PROCEDURE, NOPASS:: inside => insideQuad PROCEDURE, PASS:: gatherEF => gatherEFQuad + PROCEDURE, PASS:: gatherMF => gatherMFQuad PROCEDURE, PASS:: getNodes => getNodesQuad PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: nextElement => nextElementQuad @@ -115,6 +116,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: elemF => elemFTria PROCEDURE, NOPASS:: inside => insideTria PROCEDURE, PASS:: gatherEF => gatherEFTria + PROCEDURE, PASS:: gatherMF => gatherMFTria PROCEDURE, PASS:: getNodes => getNodesTria PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: nextElement => nextElementTria @@ -526,6 +528,33 @@ MODULE moduleMesh2DCyl END FUNCTION gatherEFQuad + PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol2DCylQuad), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: fPsi(1:4) + REAL(8):: MF_Nodes(1:4,1:3) + REAL(8):: MF(1:3) + + MF_Nodes(1:4,1) = (/self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1), & + self%n4%emData%B(1) /) + MF_Nodes(1:4,2) = (/self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2), & + self%n4%emData%B(2) /) + MF_Nodes(1:4,3) = (/self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3), & + self%n4%emData%B(3) /) + + fPsi = self%fPsi(xi) + MF = MATMUL(fPsi(:), MF_Nodes) + + END FUNCTION gatherMFQuad + !Gets nodes from quadrilateral element PURE FUNCTION getNodesQuad(self) RESULT(n) IMPLICIT NONE @@ -845,6 +874,30 @@ MODULE moduleMesh2DCyl END FUNCTION gatherEFTria + PURE FUNCTION gatherMFTria(self,xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol2DCylTria), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: fPsi(1:3) + REAL(8):: MF_Nodes(1:3,1:3) + REAL(8):: MF(1:3) + + MF_Nodes(1:3,1) = (/self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1) /) + MF_Nodes(1:3,2) = (/self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2) /) + MF_Nodes(1:3,3) = (/self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3) /) + + fPsi = self%fPsi(xi) + MF = MATMUL(fPsi, MF_Nodes) + + END FUNCTION gatherMFTria + !Gets node indexes from triangular element PURE FUNCTION getNodesTria(self) RESULT(n) IMPLICIT NONE diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 477373f..d6a6391 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -78,6 +78,7 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: elemF => elemFTetra PROCEDURE, NOPASS:: inside => insideTetra PROCEDURE, PASS:: gatherEF => gatherEFTetra + PROCEDURE, PASS:: gatherMF => gatherMFTetra PROCEDURE, PASS:: getNodes => getNodesTetra PROCEDURE, PASS:: phy2log => phy2logTetra PROCEDURE, PASS:: nextElement => nextElementTetra @@ -498,6 +499,33 @@ MODULE moduleMesh3DCart END FUNCTION gatherEFTetra + PURE FUNCTION gatherMFTetra(self, xi) RESULT(MF) + IMPLICIT NONE + + CLASS(meshVol3DCartTetra), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: fPsi(1:4) + REAL(8):: MF_Nodes(1:4,1:3) + REAL(8):: MF(1:3) + + MF_Nodes(1:4,1) = (/self%n1%emData%B(1), & + self%n2%emData%B(1), & + self%n3%emData%B(1), & + self%n4%emData%B(1) /) + MF_Nodes(1:4,2) = (/self%n1%emData%B(2), & + self%n2%emData%B(2), & + self%n3%emData%B(2), & + self%n4%emData%B(2) /) + MF_Nodes(1:4,3) = (/self%n1%emData%B(3), & + self%n2%emData%B(3), & + self%n3%emData%B(3), & + self%n4%emData%B(3) /) + + fPsi = self%fPsi(xi) + MF = MATMUL(fPsi, MF_Nodes) + + END FUNCTION gatherMFTetra + PURE FUNCTION getNodesTetra(self) RESULT(n) IMPLICIT NONE diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 index 66fe345..e63a73c 100644 --- a/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 +++ b/src/modules/mesh/inout/gmsh2/moduleMeshOutputGmsh2.f90 @@ -200,6 +200,21 @@ MODULE moduleMeshOutputGmsh2 WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref END DO WRITE(20, "(A)") '$EndElementData' + + WRITE(20, "(A)") '$NodeData' + WRITE(20, "(A)") '1' + WRITE(20, "(A)") '"Magnetic Field (T)"' + WRITE(20, *) 1 + WRITE(20, *) time + WRITE(20, *) 3 + WRITE(20, *) t + WRITE(20, *) 3 + WRITE(20, *) self%numNodes + DO n=1, self%numNodes + WRITE(20, *) n, self%nodes(n)%obj%emData%B * B_ref + END DO + WRITE(20, "(A)") '$EndNodeData' + CLOSE(20) END IF diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 19283b5..9547304 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -166,6 +166,7 @@ MODULE moduleMesh PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi PROCEDURE, PASS:: scatter PROCEDURE(gatherEF_interface), DEFERRED, PASS:: gatherEF + PROCEDURE(gatherMF_interface), DEFERRED, PASS:: gatherMF PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF PROCEDURE, PASS:: findCell @@ -194,6 +195,14 @@ MODULE moduleMesh END FUNCTION gatherEF_interface + PURE FUNCTION gatherMF_interface(self, xi) RESULT(MF) + IMPORT:: meshVol + CLASS(meshVol), INTENT(in):: self + REAL(8), INTENT(in):: xi(1:3) + REAL(8):: MF(1:3) + + END FUNCTION gatherMF_interface + PURE FUNCTION getNodesVol_interface(self) RESULT(n) IMPORT:: meshVol CLASS(meshVol), INTENT(in):: self diff --git a/src/modules/moduleEM.f90 b/src/modules/moduleEM.f90 index 876093e..ce91858 100644 --- a/src/modules/moduleEM.f90 +++ b/src/modules/moduleEM.f90 @@ -62,6 +62,23 @@ MODULE moduleEM END FUNCTION gatherElecField + PURE FUNCTION gatherMagnField(part) RESULT(BField) + USE moduleSpecies + USE moduleMesh + IMPLICIT NONE + + TYPE(particle), INTENT(in):: part + REAl(8):: xi(1:3) !Logical coordinates of particle in element + REAL(8):: BField(1:3) + + BField = 0.D0 + + xi = part%xi + + BField = mesh%vols(part%vol)%obj%gatherMF(xi) + + END FUNCTION gatherMagnField + SUBROUTINE assembleSourceVector(vectorF) USE moduleMesh USE moduleRefParam diff --git a/src/modules/moduleInput.f90 b/src/modules/moduleInput.f90 index d250954..31e2588 100644 --- a/src/modules/moduleInput.f90 +++ b/src/modules/moduleInput.f90 @@ -155,6 +155,7 @@ MODULE moduleInput Vol_ref = L_ref**3 !reference volume EF_ref = qe*n_ref*L_ref/eps_0 !reference electric field Volt_ref = EF_ref*L_ref !reference voltage + B_ref = m_ref / (ti_ref * qe) !reference magnetic field END SUBROUTINE readReference @@ -177,8 +178,9 @@ MODULE moduleInput !simulation final and initial times in [t] REAL(8):: finalTime, initialTime CHARACTER(:), ALLOCATABLE:: pusherType, WSType, EMType + REAL(8):: B(1:3) INTEGER:: nTau, nSolver - INTEGER:: i + INTEGER:: i, n CHARACTER(2):: iString CHARACTER(1):: tString @@ -267,6 +269,31 @@ MODULE moduleInput !Read BC CALL readEMBoundary(config) + CASE("ConstantB") + !Read BC + CALL readEMBoundary(config) + !Read constant magnetic field + DO i = 1, 3 + WRITE(istring, '(i2)') i + CALL config%get(object // '.B(' // istring // ')', B(i), found) + IF (.NOT. found) THEN + CALL criticalError('Constant magnetic field not provided in direction ' // iString, 'readSolver') + + END IF + + END DO + + !Non-dimensional magnetic field + B = B / B_ref + + !Assign it to the nodes + DO n =1, mesh%numNodes + mesh%nodes(n)%obj%emData%B(1) = B(1) + mesh%nodes(n)%obj%emData%B(2) = B(2) + mesh%nodes(n)%obj%emData%B(3) = B(3) + + END DO + CASE DEFAULT CALL criticalError('EM Solver ' // EMType // ' not found', 'readSolver') diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90 index 4a3f4da..ecc5462 100644 --- a/src/modules/moduleOutput.f90 +++ b/src/modules/moduleOutput.f90 @@ -11,6 +11,7 @@ MODULE moduleOutput TYPE emNode CHARACTER(:), ALLOCATABLE:: type REAL(8):: phi + REAL(8):: B(1:3) END TYPE emNode diff --git a/src/modules/moduleRefParam.f90 b/src/modules/moduleRefParam.f90 index 9d2b8d9..2dacc21 100644 --- a/src/modules/moduleRefParam.f90 +++ b/src/modules/moduleRefParam.f90 @@ -3,7 +3,7 @@ MODULE moduleRefParam !Parameters that define the problem (inputs) REAL(8):: n_ref, m_ref, T_ref, r_ref, debye_ref, sigma_ref !Reference parameters for non-dimensional problem - REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref + REAL(8):: L_ref, v_ref, ti_ref, Vol_ref, EF_ref, Volt_ref, B_ref END MODULE moduleRefParam diff --git a/src/modules/moduleSolver.f90 b/src/modules/moduleSolver.f90 index 69132d2..6e23c6f 100644 --- a/src/modules/moduleSolver.f90 +++ b/src/modules/moduleSolver.f90 @@ -27,7 +27,7 @@ MODULE moduleSolver INTERFACE !Push a particle - PURE SUBROUTINE push_interafece(part, tauIn) + SUBROUTINE push_interafece(part, tauIn) USE moduleSpecies TYPE(particle), INTENT(inout):: part @@ -80,7 +80,10 @@ MODULE moduleSolver self%pushParticle => pushCartNeutral CASE('Electrostatic') - self%pushParticle => pushCartCharged + self%pushParticle => pushCartElectrostatic + + CASE('Electromagnetic') + self%pushParticle => pushCartElectromagnetic CASE DEFAULT CALL criticalError('Pusher ' // pusherType // ' not found for Cart','initPusher') @@ -88,22 +91,17 @@ MODULE moduleSolver END SELECT CASE('Cyl') - IF (self%dimen == 2) THEN - SELECT CASE(pusherType) - CASE('Neutral') - self%pushParticle => push2DCylNeutral + SELECT CASE(pusherType) + CASE('Neutral') + self%pushParticle => push2DCylNeutral - CASE('Electrostatic') - self%pushParticle => push2DCylCharged + CASE('Electrostatic') + self%pushParticle => push2DCylElectrostatic - CASE DEFAULT - CALL criticalError('Pusher ' // pusherType // ' not found for Cyl','initPusher') + CASE DEFAULT + CALL criticalError('Pusher ' // pusherType // ' not found for Cyl','initPusher') - END SELECT - - END IF - - END SELECT + END SELECT CASE('Rad') SELECT CASE(pusherType) @@ -111,7 +109,7 @@ MODULE moduleSolver self%pushParticle => push1DRadNeutral CASE('Electrostatic') - self%pushParticle => push1DRadCharged + self%pushParticle => push1DRadElectrostatic CASE DEFAULT CALL criticalError('Pusher ' // pusherType // ' not found for Rad','initPusher') @@ -134,7 +132,7 @@ MODULE moduleSolver CHARACTER(:), ALLOCATABLE:: EMType SELECT CASE(EMType) - CASE('Electrostatic') + CASE('Electrostatic','ConstantB') self%solveEM => solveElecField END SELECT @@ -191,56 +189,71 @@ MODULE moduleSolver 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 - - !z - part_temp%v(3) = part%v(3) - part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn - - part = part_temp + part%r = part%r + part%v*tauIn END SUBROUTINE pushCartNeutral - !Push charged particles in 3D cartesian coordinates - PURE SUBROUTINE pushCartCharged(part, tauIn) + PURE SUBROUTINE pushCartElectrostatic(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 + qmEFt = part%qm*gatherElecField(part)*tauIn - !x - part_temp%v(1) = part%v(1) + qmEFt(1) - part_temp%r(1) = part%r(1) + part_temp%v(1)*tauIn + !Update velocity + part%v = part%v + qmEFt - !y - part_temp%v(2) = part%v(2) + qmEFt(2) - part_temp%r(2) = part%r(2) + part_temp%v(2)*tauIn + !Update position + part%r = part%r + part%v*tauIn - !z - part_temp%v(3) = part%v(3) + qmEFt(3) - part_temp%r(3) = part%r(3) + part_temp%v(3)*tauIn + END SUBROUTINE pushCartElectrostatic - part = part_temp + SUBROUTINE pushCartElectromagnetic(part, tauIn) + USE moduleSPecies + USE moduleEM + USE moduleMath + IMPLICIT NONE - END SUBROUTINE pushCartCharged + TYPE(particle), INTENT(inout):: part + REAL(8), INTENT(in):: tauIn + REAL(8):: tauInHalf + REAL(8):: qmEFt(1:3) + REAL(8):: B(1:3), BNorm + REAL(8):: fn + REAL(8):: v_minus(1:3), v_prime(1:3), v_plus(1:3) + + tauInHalf = tauIn *0.5D0 + !Half of the force o f the electric field + qmEFt = part%qm*gatherElecField(part)*tauInHalf + + !Half step for electrostatic + v_minus = part%v + qmEFt + + !Full step rotation + B = gatherMagnField(part) + BNorm = NORM2(B) + IF (BNorm > 0.D0) THEN + fn = DTAN(part%qm * tauInHalf*BNorm) / Bnorm + v_prime = v_minus + fn * crossProduct(v_minus, B) + v_plus = v_minus + 2.D0 * fn / (1.D0 + fn**2 * B**2)*crossProduct(v_prime, B) + + PRINT *, v_minus, v_plus + + END IF + + !Half step for electrostatic + part%v = v_plus + qmEFt + + !Update position + part%r = part%r + part%v*tauIn + + END SUBROUTINE pushCartElectromagnetic !Push one particle. Boris pusher for 2D Cyl Neutral particle PURE SUBROUTINE push2DCylNeutral(part, tauIn) @@ -279,8 +292,8 @@ MODULE moduleSolver END SUBROUTINE push2DCylNeutral - !Push one particle. Boris pusher for 2D Cyl Charged particle - PURE SUBROUTINE push2DCylCharged(part, tauIn) + !Push one particle. Boris pusher for 2D Cyl Electrostatic particle + PURE SUBROUTINE push2DCylElectrostatic(part, tauIn) USE moduleSpecies USE moduleEM IMPLICIT NONE @@ -318,7 +331,7 @@ MODULE moduleSolver !Copy temporal particle to particle part=part_temp - END SUBROUTINE push2DCylCharged + END SUBROUTINE push2DCylElectrostatic !Push one particle. Boris pusher for 1D Radial Neutral particle PURE SUBROUTINE push1DRadNeutral(part, tauIn) @@ -355,8 +368,8 @@ MODULE moduleSolver END SUBROUTINE push1DRadNeutral - !Push one particle. Boris pusher for 1D Radial Charged particle - PURE SUBROUTINE push1DRadCharged(part, tauIn) + !Push one particle. Boris pusher for 1D Radial Electrostatic particle + PURE SUBROUTINE push1DRadElectrostatic(part, tauIn) USE moduleSpecies USE moduleEM IMPLICIT NONE @@ -391,7 +404,7 @@ MODULE moduleSolver !Copy temporal particle to particle part=part_temp - END SUBROUTINE push1DRadCharged + END SUBROUTINE push1DRadElectrostatic !Dummy pusher for 0D geometry PURE SUBROUTINE push0D(part, tauIn)