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

@@ -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()
 !------------------------------------------------------------------------------

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

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,14 +780,17 @@ 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
       First = .FALSE.
       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,14 +994,17 @@ 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
       First = .FALSE.
       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,14 +1981,18 @@ 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
       First = .FALSE.
       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)