I am doing frequency-transient simulation with WhitneyAVHarmonicSolver (complex magnetic vector potential), then computing JxB with MagnetoDynamicsCalcFields and using resulting force in FlowSolver (Navier-Stokes), but it seems that Elmer does not support this functionality.
Firstly, it seems that Elmer is not computing current density for complex fields JatIP(2,l): MagnetoDynamics.F90 line 6477:
Code: Select all
IF (Vdofs == 1) THEN
...
ELSE
DO l=1,3
FORCE(p,k+l) = FORCE(p,k+l)+s*SUM( REAL(CMat_ip(l,1:3)) * E(1,1:3) )*Basis(p) - &
s * SUM( AIMAG(CMat_ip(l,1:3)) * E(2,1:3) ) * Basis(p) + s*REAL(BodyForceCurrDens_ip(l))*Basis(p)
END DO
k = k+3
DO l=1,3
FORCE(p,k+l) = FORCE(p,k+l)+s*SUM( AIMAG(CMat_ip(l,1:3)) * E(1,1:3) )*Basis(p) + &
s * SUM( REAL(CMat_ip(l,1:3)) * E(2,1:3) ) * Basis(p) + s*AIMAG(BodyForceCurrDens_ip(l))*Basis(p)
END DO
k = k+3
END IF
END IF
Code: Select all
IF (Vdofs == 1) THEN
...
ELSE
DO l=1,3
JatIp(1,l) = SUM( REAL(CMat_ip(l,1:3)) * E(1,1:3) ) - &
SUM( AIMAG(CMat_ip(l,1:3)) * E(2,1:3) ) + REAL(BodyForceCurrDens_ip(l))
FORCE(p,k+l) = FORCE(p,k+l)+s*JatIp(1,l)*Basis(p)
END DO
k = k+3
DO l=1,3
JatIp(2,l) = SUM( AIMAG(CMat_ip(l,1:3)) * E(1,1:3) ) + &
SUM( REAL(CMat_ip(l,1:3)) * E(2,1:3) ) + AIMAG(BodyForceCurrDens_ip(l))
FORCE(p,k+l) = FORCE(p,k+l)+s*JatIp(2,l)*Basis(p)
END DO
k = k+3
Code: Select all
IF ( ASSOCIATED(JXB).OR.ASSOCIATED(EL_JXB)) THEN
IF (.NOT. ASSOCIATED(CD) .AND. .NOT. ASSOCIATED(EL_CD)) THEN
CALL Warn('MagnetoDynamicsCalcFields', 'Cannot Calculate JxB since Current Density is not calculated!')
ELSE
IF (Vdofs == 1) THEN
JXBatIP(1,:) = crossproduct(JatIP(1,:),B(1,:))
DO l=1,dim
FORCE(p,k+l) = FORCE(p,k+l)+s*JXBatIP(1,l)*Basis(p)
END DO
k = k+3
END IF
END IF
END IF
Code: Select all
IF ( ASSOCIATED(JXB).OR.ASSOCIATED(EL_JXB)) THEN
IF (.NOT. ASSOCIATED(CD) .AND. .NOT. ASSOCIATED(EL_CD)) THEN
CALL Warn('MagnetoDynamicsCalcFields', 'Cannot Calculate JxB since Current Density is not calculated!')
ELSE
IF (Vdofs == 1) THEN
JXBatIP(1,:) = crossproduct(JatIP(1,:),B(1,:))
DO l=1,dim
FORCE(p,k+l) = FORCE(p,k+l)+s*JXBatIP(1,l)*Basis(p)
END DO
k = k+3
ELSE
! |X1| |Y1| |Conj(X2*Y3-X3*Y2)|
! |X2| x |Y2| = |Conj(X3*Y1-X1*Y3)|
! |X3| |Y3| |Conj(X1*Y2-X2*Y1)|
JXBatIP(1,1) = JatIP(1,2)*B(1,3) - JatIP(1,3)*B(1,2) - JatIP(2,2)*B(2,3) + JatIP(2,3)*B(2,2)
JXBatIP(1,2) = JatIP(1,3)*B(1,1) - JatIP(1,1)*B(1,3) - JatIP(2,3)*B(2,1) + JatIP(2,1)*B(2,3)
JXBatIP(1,3) = JatIP(1,1)*B(1,2) - JatIP(1,2)*B(1,1) - JatIP(2,1)*B(2,2) + JatIP(2,2)*B(2,1)
JXBatIP(2,1) = - JatIP(1,2)*B(2,3) - JatIP(2,2)*B(1,3) + JatIP(1,3)*B(2,2) + JatIP(2,3)*B(1,2)
JXBatIP(2,2) = - JatIP(1,3)*B(2,1) - JatIP(2,3)*B(1,1) + JatIP(1,1)*B(2,3) + JatIP(2,1)*B(1,3)
JXBatIP(2,3) = - JatIP(1,1)*B(2,2) - JatIP(2,1)*B(1,2) + JatIP(1,2)*B(2,1) + JatIP(2,2)*B(1,1)
Code: Select all
TYPE(Variable_t), POINTER :: JxB
REAL(KIND=dp), POINTER :: JxBvalues(:)
INTEGER, POINTER :: JxBperm(:)
INTEGER :: JxBdofs
JXB => NULL()
JxB => VariableGet( CurrentModel % Variables, 'JxB' )
IF ( .NOT.ASSOCIATED( JxB ) ) THEN
CALL Fatal( 'Differentials', 'Cannot access JxB from MagnetoDynamics.' )
ELSE
JxBvalues => JxB % Values
JxBperm => JxB % Perm
JxBdofs = JxB % Dofs
stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, Basis,dBasisdx )
L(1) = SUM( Basis(1:n)* JxBvalues( JxBdofs * (JxBperm(NodeIndexes)-1)+1 ) )
L(2) = SUM( Basis(1:n)* JxBvalues( JxBdofs * (JxBperm(NodeIndexes)-1)+2 ) )
L(3) = SUM( Basis(1:n)* JxBvalues( JxBdofs * (JxBperm(NodeIndexes)-1)+3 ) )
RETURN
END IF