From 6d29a12e91c319f46f53b4951bde7cbd40dff872 Mon Sep 17 00:00:00 2001 From: Peter Raback Date: Tue, 29 Oct 2024 18:42:01 +0200 Subject: [PATCH] The CircuitsAndDynamics solver was not complaining at all even if 'W' was not found. This may end up lost time in debugging. Now at least a Warning is given. --- fem/src/CircuitUtils.F90 | 23 +++++++++++ fem/src/Lists.F90 | 7 +++- fem/src/modules/CircuitsAndDynamics.F90 | 51 +++++++++++++++++++------ 3 files changed, 67 insertions(+), 14 deletions(-) diff --git a/fem/src/CircuitUtils.F90 b/fem/src/CircuitUtils.F90 index 9c39b424a1..8c708484f6 100644 --- a/fem/src/CircuitUtils.F90 +++ b/fem/src/CircuitUtils.F90 @@ -157,6 +157,29 @@ SUBROUTINE GetWPotential(Wbase) END SUBROUTINE GetWPotential !------------------------------------------------------------------------------ + +!------------------------------------------------------------------------------ + SUBROUTINE GetWPotentialVar(pVar) +!------------------------------------------------------------------------------ + IMPLICIT NONE + + TYPE(Variable_t), POINTER :: pVar + + pVar => VariableGet( CurrentModel % Mesh % Variables,'W Potential') + IF(.NOT. ASSOCIATED(pVar) ) THEN + pVar => VariableGet( CurrentModel % Mesh % Variables,'W') + END IF + IF(ASSOCIATED(pVar)) THEN + CALL Info('GetWPotentialVar','Using gradient of field to define direction: '& + //TRIM(pVar % Name),Level=7) + ELSE + CALL Warn('GetWPotentialVar','Could not obtain variable for potential "W"') + END IF +!------------------------------------------------------------------------------ + END SUBROUTINE GetWPotentialVar +!------------------------------------------------------------------------------ + + !------------------------------------------------------------------------------ SUBROUTINE AddComponentsToBodyLists() !------------------------------------------------------------------------------ diff --git a/fem/src/Lists.F90 b/fem/src/Lists.F90 index 430ab4e941..cb59026771 100644 --- a/fem/src/Lists.F90 +++ b/fem/src/Lists.F90 @@ -8700,8 +8700,11 @@ FUNCTION ListGetElementVectorSolution( Handle, Basis, Element, Found, GaussPoint Val3D = 0.0_dp - IF( .NOT. ASSOCIATED( Handle % Variable ) ) RETURN - + IF( .NOT. ASSOCIATED( Handle % Variable ) ) THEN + IF(PRESENT(Found)) Found = .FALSE. + RETURN + END IF + IF( PRESENT( dofs ) ) THEN Ldofs = dofs ELSE diff --git a/fem/src/modules/CircuitsAndDynamics.F90 b/fem/src/modules/CircuitsAndDynamics.F90 index 9648628158..a08d0f445a 100644 --- a/fem/src/modules/CircuitsAndDynamics.F90 +++ b/fem/src/modules/CircuitsAndDynamics.F90 @@ -574,7 +574,8 @@ SUBROUTINE Add_stranded(Element,Tcoef,Comp,nn,nd,dt,CompParams) REAL(KIND=dp) :: WBasis(nd,3), RotWBasis(nd,3) INTEGER :: dim, ncdofs,q TYPE(VariableHandle_t), SAVE :: Wvec_h - + TYPE(Variable_t), POINTER, SAVE :: Wpot + LOGICAL :: PiolaVersion = .FALSE. SAVE CSymmetry, dim, First, InitHandle @@ -622,6 +623,8 @@ SUBROUTINE Add_stranded(Element,Tcoef,Comp,nn,nd,dt,CompParams) CALL ListInitElementVariable(Wvec_h, CoilWVecVarname) PiolaVersion = GetLogical( ASolver % Values, 'Use Piola Transform', Found ) + + CALL GetWPotentialVar(Wpot) END IF PS => Asolver % Variable % Perm @@ -644,14 +647,14 @@ SUBROUTINE Add_stranded(Element,Tcoef,Comp,nn,nd,dt,CompParams) ncdofs=nd IF (dim == 3) THEN - + ! We can choose the base per component. CoilUseWvec = GetLogical(CompParams, 'Coil Use W Vector', Found) IF (.NOT. Found) CoilUseWvec = CoilUseWvec0 IF (.NOT. CoilUseWvec) THEN !CALL GetLocalSolution(Wbase, 'w') ! when W Potential solver is used, 'w' is not enough. - CALL GetWPotential(WBase) + CALL GetLocalSolution( Wbase,UElement=Element,UVariable=Wpot, Found=Found) END IF ncdofs=nd-nn @@ -777,7 +780,8 @@ SUBROUTINE Add_massive(Element,Tcoef,Comp,nn,nd,dt,crt) REAL(KIND=dp) :: wBase(nn), gradv(3), WBasis(nd,3), RotWBasis(nd,3) INTEGER :: ncdofs,q - + TYPE(Variable_t), POINTER, SAVE :: Wpot + SAVE CSymmetry, dim, First IF (First) THEN @@ -785,6 +789,8 @@ SUBROUTINE Add_massive(Element,Tcoef,Comp,nn,nd,dt,crt) CSymmetry = ( CurrentCoordinateSystem() == AxisSymmetric .OR. & CurrentCoordinateSystem() == CylindricSymmetric ) dim = CoordinateSystemDimension() + + CALL GetWPotentialVar(Wpot) END IF ASolver => CurrentModel % Asolver @@ -835,7 +841,8 @@ SUBROUTINE Add_massive(Element,Tcoef,Comp,nn,nd,dt,crt) ncdofs=nd IF (dim == 3) THEN - CALL GetLocalSolution(Wbase, 'w') + !CALL GetLocalSolution(Wbase, 'w') + CALL GetLocalSolution( Wbase,UElement=Element,UVariable=Wpot, Found=Found) ncdofs=nd-nn END IF @@ -987,7 +994,8 @@ SUBROUTINE Add_foil_winding(Element,Tcoef,Comp,nn,nd,dt) REAL(KIND=dp) :: wBase(nn), gradv(3), WBasis(nd,3), RotWBasis(nd,3), & RotMLoc(3,3), RotM(3,3,nn) INTEGER :: i,ncdofs,q - + TYPE(Variable_t), POINTER, SAVE :: Wpot + SAVE CSymmetry, dim, First IF (First) THEN @@ -995,6 +1003,8 @@ SUBROUTINE Add_foil_winding(Element,Tcoef,Comp,nn,nd,dt) CSymmetry = ( CurrentCoordinateSystem() == AxisSymmetric .OR. & CurrentCoordinateSystem() == CylindricSymmetric ) dim = CoordinateSystemDimension() + + CALL GetWPotentialVar(Wpot) END IF ASolver => CurrentModel % Asolver @@ -1020,7 +1030,8 @@ SUBROUTINE Add_foil_winding(Element,Tcoef,Comp,nn,nd,dt) ncdofs=nd IF (dim == 3) THEN - CALL GetLocalSolution(Wbase, 'w') + CALL GetLocalSolution( Wbase,UElement=Element,UVariable=Wpot, Found=Found) + !CALL GetLocalSolution(Wbase, 'w') CALL GetElementRotM(Element, RotM, nn) ncdofs=nd-nn END IF @@ -1769,7 +1780,8 @@ SUBROUTINE Add_stranded(Element,Tcoef,Comp,nn,nd,CompParams) LOGICAL :: PiolaVersion = .FALSE. TYPE(VariableHandle_t), SAVE :: Wvec_h - + TYPE(Variable_t), POINTER, SAVE :: Wpot + SAVE CSymmetry, dim, First IF (First) THEN @@ -1804,6 +1816,7 @@ SUBROUTINE Add_stranded(Element,Tcoef,Comp,nn,nd,CompParams) IF(.NOT. Found) CoilWVecVarname = 'W Vector E' CALL ListInitElementVariable(Wvec_h, CoilWVecVarname) + CALL GetWPotentialVar(Wpot) END IF ASolver => CurrentModel % Asolver @@ -1827,7 +1840,8 @@ SUBROUTINE Add_stranded(Element,Tcoef,Comp,nn,nd,CompParams) IF(.NOT. Found) CoilUseWvec = CoilUseWvec0 IF (.NOT. CoilUseWvec) THEN - CALL GetWPotential(WBase) + CALL GetLocalSolution( Wbase,UElement=Element,UVariable=Wpot, Found=Found) + !CALL GetWPotential(WBase) END IF END IF @@ -1874,7 +1888,10 @@ SUBROUTINE Add_stranded(Element,Tcoef,Comp,nn,nd,CompParams) CALL GetEdgeBasis(Element,WBasis,RotWBasis,Basis,dBasisdx) END IF IF (CoilUseWvec) THEN - w = ListGetElementVectorSolution( Wvec_h, Basis, Element, dofs = dim ) + w = ListGetElementVectorSolution( Wvec_h, Basis, Element, Found = Found, dofs = dim ) + IF(.NOT. Found ) THEN + CALL Fatal('Add_stranded','Could not find coil current density!') + END IF ELSE w = -MATMUL(WBase(1:nn), dBasisdx(1:nn,:)) END IF @@ -1964,7 +1981,9 @@ SUBROUTINE Add_massive(Element,Tcoef,Comp,nn,nd) INTEGER :: ncdofs,q REAL(KIND=dp) :: ModelDepth COMPLEX(KIND=dp) :: Permittivity(nn), localP + TYPE(Variable_t), POINTER, SAVE :: Wpot + SAVE CSymmetry, dim, First IF (First) THEN @@ -1972,6 +1991,8 @@ SUBROUTINE Add_massive(Element,Tcoef,Comp,nn,nd) CSymmetry = ( CurrentCoordinateSystem() == AxisSymmetric .OR. & CurrentCoordinateSystem() == CylindricSymmetric ) dim = CoordinateSystemDimension() + + CALL GetWPotentialVar(Wpot) END IF ASolver => CurrentModel % Asolver @@ -1990,7 +2011,8 @@ SUBROUTINE Add_massive(Element,Tcoef,Comp,nn,nd) ncdofs=nd IF (dim == 3) THEN - CALL GetWPotential(WBase) + !CALL GetWPotential(WBase) + CALL GetLocalSolution( Wbase,UElement=Element,UVariable=Wpot, Found=Found) ncdofs=nd-nn END IF @@ -2209,7 +2231,9 @@ SUBROUTINE Add_foil_winding(Element,Tcoef,Comp,nn,nd,CompParams) RotMLoc(3,3), RotM(3,3,nn) REAL(KIND=dp) :: Jvec(3) INTEGER :: i,ncdofs,q + TYPE(Variable_t), POINTER, SAVE :: Wpot + SAVE CSymmetry, dim, First, InitHandle, InitJHandle IF( First ) THEN @@ -2243,6 +2267,8 @@ SUBROUTINE Add_foil_winding(Element,Tcoef,Comp,nn,nd,CompParams) END IF IF(.NOT. Found) CoilWVecVarname = 'W Vector E' CALL ListInitElementVariable(Wvec_h, CoilWVecVarname) + + CALL GetWPotentialVar(Wpot) END IF ASolver => CurrentModel % Asolver @@ -2266,7 +2292,8 @@ SUBROUTINE Add_foil_winding(Element,Tcoef,Comp,nn,nd,CompParams) IF (.NOT. CoilUseWvec) THEN !CALL GetLocalSolution(Wbase, 'w') - CALL GetWPotential(WBase) + !CALL GetWPotential(WBase) + CALL GetLocalSolution( Wbase,UElement=Element,UVariable=Wpot, Found=Found) END IF FoilUseJvec = GetLogical(CompParams, 'Foil Winding Use J Vector', Found)