From 757916e55092a522b5b765ead3df63c4ee16d60c Mon Sep 17 00:00:00 2001 From: Juha Ruokolainen Date: Tue, 27 Aug 2024 12:50:16 +0300 Subject: [PATCH] Several small fixes. --- fem/src/MeshUtils.F90 | 24 ++++++------ fem/src/RadiationFactors.F90 | 2 +- fem/src/SolverUtils.F90 | 38 +++++++++---------- fem/src/modules/MagnetoDynamics2D.F90 | 12 ++++-- fem/tests/mgdyn_anisotropic_cond/currents.f90 | 6 +-- .../currents.f90 | 6 +-- 6 files changed, 48 insertions(+), 40 deletions(-) diff --git a/fem/src/MeshUtils.F90 b/fem/src/MeshUtils.F90 index 904ec5224c..f6ede65fc4 100644 --- a/fem/src/MeshUtils.F90 +++ b/fem/src/MeshUtils.F90 @@ -10956,7 +10956,6 @@ SUBROUTINE AddProjectorWeakGeneric() END DO END IF - DO ind=1,BMesh1 % NumberOfBulkElements ! Optionally save the submesh for specified element, for visualization and debugging @@ -11139,7 +11138,6 @@ SUBROUTINE AddProjectorWeakGeneric() ElemCands = 0 ElemHits = 0 - DO indM=1,BMesh2 % NumberOfBulkElements ElementM => BMesh2 % Elements(indM) @@ -12823,16 +12821,20 @@ FUNCTION SegmentOriginDistance2(x1,y1,x2,y2) RESULT ( r2 ) REAL(KIND=dp) :: x1,y1,x2,y2,r2 REAL(KIND=dp) :: q,xc,yc - q = ( x1*(x1-x2) + y1*(y1-y2) ) / & - SQRT((x1**2+y1**2) * ((x1-x2)**2+(y1-y2)**2)) - IF( q <= 0.0_dp ) THEN - r2 = x1**2 + y1**2 - ELSE IF( q >= 1.0_dp ) THEN - r2 = x2**2 + y2**2 + IF(x1==0.AND.y1==0) THEN + r2 = 0; ELSE - xc = x1 + q * (x2-x1) - yc = y1 + q * (y2-y1) - r2 = xc**2 + yc**2 + q = ( x1*(x1-x2) + y1*(y1-y2) ) / & + SQRT((x1**2+y1**2) * ((x1-x2)**2+(y1-y2)**2)) + IF( q <= 0.0_dp ) THEN + r2 = x1**2 + y1**2 + ELSE IF( q >= 1.0_dp ) THEN + r2 = x2**2 + y2**2 + ELSE + xc = x1 + q * (x2-x1) + yc = y1 + q * (y2-y1) + r2 = xc**2 + yc**2 + END IF END IF END FUNCTION SegmentOriginDistance2 diff --git a/fem/src/RadiationFactors.F90 b/fem/src/RadiationFactors.F90 index 44636c3a89..56cc75eeb6 100644 --- a/fem/src/RadiationFactors.F90 +++ b/fem/src/RadiationFactors.F90 @@ -1280,7 +1280,7 @@ SUBROUTINE CalculateGebhartFactors() END IF ! Scale matrix to unit diagonals - Diag = SQRT(1._dp/ABS(Diag)) + Diag = SQRT(1._dp/MAX(ABS(Diag),1.0d-12)) DO i=1,RadiationSurfaces DO j=G % Rows(i),G % Rows(i+1)-1 G % Values(j) = G % Values(j)*Diag(i)*Diag(G % Cols(j)) diff --git a/fem/src/SolverUtils.F90 b/fem/src/SolverUtils.F90 index 6828d09339..6ece18a1c1 100644 --- a/fem/src/SolverUtils.F90 +++ b/fem/src/SolverUtils.F90 @@ -762,7 +762,12 @@ SUBROUTINE UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, & np = mGetElementDOFs(pIndexes,Element) np = MIN(np, n) - NormalIndexes(1:np) = NT % BoundaryReorder(pIndexes(1:np)) + DO i=1,np + j = pIndexes(i) + IF(j>0 .AND. j<= SIZE(NT % BoundaryReorder)) THEN + NormalIndexes(i) = NT % BoundaryReorder(j) + END IF + END DO CALL RotateMatrix( LocalStiffMatrix, LocalForce, n, dim, NDOFs, & NormalIndexes, NT % BoundaryNormals, NT % BoundaryTangent1, NT % BoundaryTangent2 ) @@ -816,7 +821,7 @@ SUBROUTINE UpdateGlobalEquationsVec( Gmtr, Lmtr, Gvec, Lvec, n, & ! Local variables INTEGER :: dim, i,j,k,np - INTEGER :: Ind(n*NDOFs),pIndexes(64) + INTEGER :: NormalIndexes(n),pIndexes(64) REAL(KIND=dp) :: Vals(n*NDOFs) !DIR$ ATTRIBUTES ALIGN:64::Ind, Vals @@ -848,22 +853,20 @@ SUBROUTINE UpdateGlobalEquationsVec( Gmtr, Lmtr, Gvec, Lvec, n, & END IF IF ( Rotate ) THEN - Ind = 0 + NormalIndexes = 0 np = mGetElementDOFs(pIndexes,Element) np = MIN(np,n) DO i=1,np j = pIndexes(i) - IF(j > 0 .AND. j <= SIZE(NT % BoundaryReorder) ) THEN - Ind(i) = NT % BoundaryReorder(j) - ELSE - Ind(i) = j + IF(j>0 .AND. j<= SIZE(NT % BoundaryReorder)) THEN + NormalIndexes(i) = NT % BoundaryReorder(j) END IF END DO ! TODO: See that RotateMatrix is vectorized - CALL RotateMatrix( Lmtr, Lvec, n, dim, NDOFs, Ind, NT % BoundaryNormals, & + CALL RotateMatrix( Lmtr, Lvec, n, dim, NDOFs, NormalIndexes, NT % BoundaryNormals, & NT % BoundaryTangent1, NT % BoundaryTangent2 ) !IF ( Rotate .AND. NormalTangentialNOFNodes > 0 .AND. ndofs>=dim) THEN @@ -989,7 +992,13 @@ SUBROUTINE UpdateGlobalForce(ForceVector, LocalForce, n, & np = mGetElementDOFs(pIndexes,Element) np = MIN(np,n) - NormalIndexes(1:np) = NT % BoundaryReorder(pIndexes(1:np)) + + DO i=1,np + j = pIndexes(i) + IF(j>0 .AND. j<= SIZE(NT % BoundaryReorder)) THEN + NormalIndexes(i) = NT % BoundaryReorder(j) + END IF + END DO CALL RotateMatrix( LocalStiffMatrix, LocalForce, n, dim, NDOFs, & NormalIndexes, NT % BoundaryNormals, NT % BoundaryTangent1, NT % BoundaryTangent2 ) @@ -5017,8 +5026,7 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & Conditional = ActiveCond(BC) Element => Mesh % Elements(t) - IF ( Element % BoundaryInfo % Constraint /= & - Model % BCs(BC) % Tag ) CYCLE + IF ( Element % BoundaryInfo % Constraint /= Model % BCs(BC) % Tag ) CYCLE Model % CurrentElement => Element IF ( ActivePart(BC) ) THEN @@ -5045,8 +5053,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & END IF END IF - - BLOCK INTEGER,ALLOCATABLE :: ChildBCs(:) INTEGER,POINTER :: BCInds(:) @@ -5215,7 +5221,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & END IF END DO - ! Check the boundaries and body forces for possible single nodes BCs that are used to fixed ! the domain for undetermined equations. The loop is slower than optimal in the case that there is ! a large amount of different boundaries that have a node to set. @@ -5321,7 +5326,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & END DO END IF - !------------------------------------------------------------------------------ ! Take care of the matrix entries of passive elements !------------------------------------------------------------------------------ @@ -5376,7 +5380,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & DEALLOCATE(PassPerm,NodeIndexes) END IF - ! Check the boundaries and body forces for possible single nodes BCs that must have a constant ! value on that boundary / body force. !-------------------------------------------------------------------------------------------- @@ -5578,7 +5581,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & //I2S(SIZE( A % Cols )),Level=8) END IF - ! Check the boundaries for a curve bc. !-------------------------------------------------------------------------------------------- BLOCK @@ -5781,7 +5783,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & END BLOCK - ! Check the boundaries and body forces for possible single nodes BCs that must have a constant ! value on that boundary / body force. !-------------------------------------------------------------------------------------------- @@ -5887,7 +5888,6 @@ SUBROUTINE SetDirichletBoundaries( Model, A, b, Name, DOF, NDOFs, Perm, & //I2S(DirCount),Level=12) END IF - !------------------------------------------------------------------------------ CONTAINS diff --git a/fem/src/modules/MagnetoDynamics2D.F90 b/fem/src/modules/MagnetoDynamics2D.F90 index 7d4ecaa884..a0c6ba012b 100755 --- a/fem/src/modules/MagnetoDynamics2D.F90 +++ b/fem/src/modules/MagnetoDynamics2D.F90 @@ -3308,7 +3308,7 @@ SUBROUTINE BulkAssembly() BAtIp(8) = AIMAG(imag_value) imag_value = CMPLX(BatIp(1), BatIp(3), KIND=dp) imag_value2 = CMPLX(BatIp(2), BatIp(4), KIND=dp) - BMagnAtIP = SQRT(ABS(imag_value**2._dp) + ABS(imag_value2**2._dp)) + BMagnAtIP = SQRT(ABS(imag_value*imag_value) + ABS(imag_value2*imag_value2)) END IF IF (LorentzForceCompute) THEN @@ -3320,8 +3320,14 @@ SUBROUTINE BulkAssembly() By = CMPLX(BatIp(2), BatIp(4), KIND=dp) Jz = CMPLX(BatIp(7), BatIp(8), KIND=dp) - LorentzForceDensX = ModelDepth * Weight * By / Jz * ABS(Jz)**2._dp - LorentzForceDensY = -ModelDepth * Weight * Bx / Jz * ABS(Jz)**2._dp + IF(Jz/=0) THEN + LorentzForceDensX = ModelDepth * Weight * By / Jz * ABS(Jz)*ABS(Jz) + LorentzForceDensY = -ModelDepth * Weight * Bx / Jz * ABS(Jz)*ABS(Jz) + ELSE + LorentzForceDensX = 0 + LorentzForceDensY = 0 + END IF + BodyLorentzForcesRe(1, BodyId) = BodyLorentzForcesRe(1, BodyId) + & REAL(LorentzForceDensX) BodyLorentzForcesRe(2, BodyId) = BodyLorentzForcesRe(2, BodyId) + & diff --git a/fem/tests/mgdyn_anisotropic_cond/currents.f90 b/fem/tests/mgdyn_anisotropic_cond/currents.f90 index bbcae66226..1021830536 100644 --- a/fem/tests/mgdyn_anisotropic_cond/currents.f90 +++ b/fem/tests/mgdyn_anisotropic_cond/currents.f90 @@ -57,7 +57,7 @@ FUNCTION currdens1( model, n, args) RESULT(curr) r = sqrt(x**2 + z**2) - theta = atan(x/z) + theta = atan2(x,z) curr = currentInToroidR(r, r0, y, turns, I) * sin(theta) @@ -87,7 +87,7 @@ FUNCTION currdens2( model, n, args) RESULT(curr) r = sqrt(x**2 + z**2) - theta = atan(x/z) + theta = atan2(x,z) curr = currentInToroidY(r, r0, y, turns, I) @@ -117,7 +117,7 @@ FUNCTION currdens3( model, n, args) RESULT(curr) r = sqrt(x**2 + z**2) - theta = atan(x/z) + theta = atan2(x,z) curr = currentInToroidR(r, r0, y, turns, I) * cos(theta) diff --git a/fem/tests/mgdyn_lamstack_lowfreq_transient/currents.f90 b/fem/tests/mgdyn_lamstack_lowfreq_transient/currents.f90 index 31e3c8fe09..1166f49a10 100644 --- a/fem/tests/mgdyn_lamstack_lowfreq_transient/currents.f90 +++ b/fem/tests/mgdyn_lamstack_lowfreq_transient/currents.f90 @@ -57,7 +57,7 @@ FUNCTION currdens1( model, n, args) RESULT(curr) r = sqrt(x**2 + z**2) - theta = atan(x/z) + theta = atan2(x,z) curr = currentInToroidR(r, r0, y, turns, I) * sin(theta) @@ -87,7 +87,7 @@ FUNCTION currdens2( model, n, args) RESULT(curr) r = sqrt(x**2 + z**2) - theta = atan(x/z) + theta = atan2(x,z) curr = currentInToroidY(r, r0, y, turns, I) @@ -117,7 +117,7 @@ FUNCTION currdens3( model, n, args) RESULT(curr) r = sqrt(x**2 + z**2) - theta = atan(x/z) + theta = atan2(x,z) curr = currentInToroidR(r, r0, y, turns, I) * cos(theta)