Skip to content

Commit

Permalink
[skip ci] Fix for projection in case Projector % InvPerm has zero val…
Browse files Browse the repository at this point in the history
…ues.
  • Loading branch information
raback committed Oct 11, 2024
1 parent e080d89 commit 89b896c
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 21 deletions.
31 changes: 20 additions & 11 deletions fem/src/modules/SaveData/SaveProjection.F90
Original file line number Diff line number Diff line change
Expand Up @@ -82,8 +82,9 @@ SUBROUTINE SaveProjection( Model,Solver,dt,Transient )
INTEGER, POINTER :: UnitPerm(:)
REAL(KIND=dp) :: Nrm
LOGICAL :: Found, Normalize, ToSlave, ToMaster
INTEGER, POINTER :: ActiveProjs(:)

CALL Info('SaveProjection','Creating selected projected values as fields')
CALL Info('SaveProjection','Creating selected projected values as fields',Level=8)

Params => GetSolverParams()

Expand All @@ -105,7 +106,7 @@ SUBROUTINE SaveProjection( Model,Solver,dt,Transient )
DO i=1,NoVar
VarName = ListGetString( Params,'Variable '//I2S(i), Found )
Var => VariableGet( Model % Variables, TRIM(VarName) )
CALL info('SaveProjection','Doing variable: '//TRIM(VarName))
CALL info('SaveProjection','Doing variable: '//TRIM(VarName),Level=8)

TargetName = ListGetString( Params,'Target Variable '//I2S(i), Found )
IF(.NOT. Found) TargetName = 'Projection '//TRIM(VarName)
Expand All @@ -132,6 +133,8 @@ SUBROUTINE SaveProjection( Model,Solver,dt,Transient )
TargetVar => VariableGet( Model % Variables, TRIM(TargetName) )
END IF

ActiveProjs => ListGetIntegerArray( Params,'Active Projectors '//I2S(i),Found )

! Do additive projection!
CALL ProjectToVariable()
Nrm = Nrm + SUM(TargetVar % Values**2)
Expand All @@ -156,7 +159,7 @@ SUBROUTINE ProjectToVariable()
REAL(KIND=dp) :: r1
REAL(KIND=dp), POINTER :: Values(:)
REAL(KIND=dp), ALLOCATABLE :: Weight(:)
LOGICAL, ALLOCATABLE :: IsInvInvPerm(:)
LOGICAL, POINTER :: IsInvInvPerm(:)

dofs = Var % Dofs
TargetVar % Values = 0.0_dp
Expand All @@ -165,19 +168,25 @@ SUBROUTINE ProjectToVariable()
ALLOCATE(Weight(SIZE(TargetVar % Values)))
Weight = 0.0_dp
END IF

! Go through all the projectors.
! There could be perhaps reason to skip some, but this will do for now.
DO bc=1,Model % NumberOfBCs
IF(ASSOCIATED(ActiveProjs)) THEN
IF(.NOT. ANY(ActiveProjs == bc)) CYCLE
END IF

A => CurrentModel % BCs(bc) % PMatrix
IF(.NOT. ASSOCIATED(A) ) THEN
A => Solver % MortarBCs(bc) % Projector
END IF

IF(.NOT. ASSOCIATED(A)) CYCLE
n = A % NumberOfRows
IF(n==0) CYCLE

CALL Info('SaveProjection','Doing projection for BC '//I2S(bc)//' of size '//I2S(n),Level=20)

Rows => A % Rows
Cols => A % Cols
Values => A % Values
Expand All @@ -191,12 +200,14 @@ SUBROUTINE ProjectToVariable()
IsInvInvPerm = .FALSE.
DO i=1,n
pi = A % InvPerm(i)
IsInvInvPerm(pi) = .TRUE.
IF(pi > 0 ) IsInvInvPerm(pi) = .TRUE.
END DO


DO k=1,dofs
DO i=1,n
pi = A % InvPerm(i)
IF(pi == 0) CYCLE
acti = IsInvInvPerm(pi)
IF(ASSOCIATED(Var % Perm)) pi = Var % Perm(pi)
IF(pi==0) CYCLE
Expand Down Expand Up @@ -232,18 +243,16 @@ SUBROUTINE ProjectToVariable()
IF(Normalize .AND. .NOT. actj) Weight(pj) = Weight(pj) + Values(j)
END IF
END IF
END DO

END DO
END DO
END DO

IF(ALLOCATED(IsInvInvPerm)) DEALLOCATE(IsInvInvPerm)
DEALLOCATE(IsInvInvPerm)
END DO

IF(Normalize) THEN
WHERE(ABS(Weight) > EPSILON(r1) )
TargetVar % Values = TargetVar % Values / Weight
END WHERE
END WHERE
END IF

END SUBROUTINE ProjectToVariable
Expand Down
12 changes: 2 additions & 10 deletions fem/tests/MortarSaveProjection/case.sif
Original file line number Diff line number Diff line change
Expand Up @@ -21,18 +21,10 @@ Simulation

Output Intervals = 0

Post File = case.vtu
Ascii Output = True
! Post File = case.vtu
! Ascii Output = True
End

Constants
Gravity(4) = 0 -1 0 9.82
Stefan Boltzmann = 5.67e-08
Permittivity of Vacuum = 8.8542e-12
Boltzmann Constant = 1.3807e-23
Unit Charge = 1.602e-19
End

Body 1
Target Bodies(1) = 1
Name = "Body1"
Expand Down

0 comments on commit 89b896c

Please sign in to comment.