Tentative steps for adjoint rhs

ElmerCSC · Jan 8, 2025 · c7120f1 · c7120f1
1 parent dc58b4f
commit c7120f1
Show file tree

Hide file tree

Showing 3 changed files with 68 additions and 22 deletions.
diff --git a/fem/src/DefUtils.F90 b/fem/src/DefUtils.F90
@@ -3511,6 +3511,10 @@ RECURSIVE SUBROUTINE DefaultInitialize( USolver, UseConstantBulk )
 
      CALL InitializeToZero( Solver % Matrix, Solver % Matrix % RHS )
 
+     IF(ASSOCIATED(Solver % Matrix % RhsAdjoint) ) THEN
+       Solver % Matrix % RhsAdjoint = 0.0_dp
+     END IF
+
      IF( ALLOCATED(Solver % Matrix % ConstrainedDOF) ) THEN
        Solver % Matrix % ConstrainedDOF = .FALSE.
      END IF
@@ -3606,6 +3610,13 @@ RECURSIVE SUBROUTINE DefaultStart( USolver )
 
        Solver % LocalSystemMode = 1
      END IF
+
+     IF(ListGetLogical( Params,'Solve Adjoint Equation',Found ) ) THEN
+       IF(.NOT. ASSOCIATED( Solver % Matrix % RhsAdjoint ) ) THEN
+         ALLOCATE( Solver % Matrix % RhsAdjoint(SIZE(Solver % Matrix % Rhs)))
+       END IF     
+       CALL ListAddLogical( Params,'Constraint Modes Analysis Frozen',.TRUE.)
+     END IF      
 
 !------------------------------------------------------------------------------
    END SUBROUTINE DefaultStart
@@ -3655,25 +3666,34 @@ RECURSIVE SUBROUTINE DefaultFinish( USolver )
            REAL(KIND=dp) :: Norm
            REAL(KIND=dp), ALLOCATABLE :: xtmp(:), btmp(:)
 
-           CALL ListAddLogical(Params,'Constraint Modes Analysis Frozen',.FALSE.)
-           CALL ListAddLogical(Params,'Skip Compute Nonlinear Change',.TRUE.)
-
            n = SIZE(Solver % Matrix % rhs)
-           ALLOCATE(xtmp(n),btmp(n))
-
-           btmp = Solver % Matrix % rhs
-           xtmp = Solver % Variable % Values 
+           CALL ListAddLogical(Params,'Skip Compute Nonlinear Change',.TRUE.)
 
-           Solver % Matrix % rhs = 0.0_dp
-           CALL SolveSystem( Solver % Matrix, ParMatrix, Solver % Matrix % rhs, &
-               Solver % Variable % Values, Norm, Solver % Variable % DOFs,Solver )
-
-           CALL ListAddLogical(Params,'Constraint Modes Analysis Frozen',.TRUE.)
+           IF( ListGetLogical(Params,'Solve Adjoint Equation', Found ) ) THEN
+
+             CALL SolveSystem( Solver % Matrix, ParMatrix, Solver % Matrix % rhsAdjoint, &
+                 Solver % Variable % ConstraintModes(1,:), Norm, Solver % Variable % DOFs,Solver )
+
+           ELSE           
+             ALLOCATE(xtmp(n),btmp(n))
+             btmp = Solver % Matrix % rhs
+             xtmp = Solver % Variable % Values 
+
+             Solver % Matrix % rhs = 0.0_dp
+
+             CALL ListAddLogical(Params,'Constraint Modes Analysis Frozen',.FALSE.)
+
+             CALL SolveSystem( Solver % Matrix, ParMatrix, Solver % Matrix % rhs, &
+                 Solver % Variable % Values, Norm, Solver % Variable % DOFs,Solver )
+
+             CALL ListAddLogical(Params,'Constraint Modes Analysis Frozen',.TRUE.)
+
+             Solver % Matrix % rhs = btmp 
+             Solver % Variable % Values = xtmp
+             DEALLOCATE(xtmp,btmp)
+           END IF
+
            CALL ListAddLogical(Params,'Skip Compute Nonlinear Change',.FALSE.)
-
-           Solver % Matrix % rhs = btmp 
-           Solver % Variable % Values = xtmp
-           DEALLOCATE(xtmp,btmp)
          END BLOCK
        END IF
 
@@ -4498,6 +4518,7 @@ SUBROUTINE DefaultUpdateForceR( F, UElement, USolver, BulkUpdate )
   END SUBROUTINE DefaultUpdateForceR
 !------------------------------------------------------------------------------
 
+
 
 !> Adds the elementwise contribution the right-hand-side of the complex valued matrix equation 
 !------------------------------------------------------------------------------
@@ -4692,6 +4713,8 @@ END SUBROUTINE DefaultUpdateTimeForceC
 
 
 
+
+
 !------------------------------------------------------------------------------
   SUBROUTINE DefaultUpdatePrecR( M, UElement, USolver ) 
 !------------------------------------------------------------------------------

diff --git a/fem/src/Types.F90 b/fem/src/Types.F90
@@ -219,6 +219,8 @@ MODULE Types
     REAL(KIND=dp), POINTER CONTIG :: RHS(:)=>NULL(),BulkRHS(:)=>NULL(),RHS_im(:)=>NULL(),Force(:,:)=>NULL()
     REAL(KIND=dp), POINTER CONTIG :: BulkResidual(:)=>NULL()
 
+    REAL(KIND=dp), POINTER CONTIG :: RhsAdjoint(:)=>NULL()
+
     REAL(KIND=dp),  POINTER CONTIG :: Values(:)=>NULL(), ILUValues(:)=>NULL(), &
                DiagScaling(:) => NULL(), TValues(:) => NULL(), Values_im(:) => NULL()
 

diff --git a/fem/src/modules/MagnetoDynamics2D.F90 b/fem/src/modules/MagnetoDynamics2D.F90
@@ -407,17 +407,17 @@ SUBROUTINE CalculateLumpedTransient(Torque)
 
    REAL(KIND=dp) :: torq,TorqArea,IMoment,IA, &
        rinner,router,rmean,rdiff,ctorq,detJ,Weight,&
-       Bp,Br,Bx,By,x,y,r,rho
+       Bp,Br,Bx,By,x,y,r,rho,wtorq,Bp0,Br0,Bx0,By0
    REAL(KIND=dp), ALLOCATABLE :: a(:),u(:),POT(:),dPOT(:), &
-       pPot(:),Density(:)
+       pPot(:),Density(:),TorqSens(:)
    REAL(KIND=dp), POINTER :: Basis(:), dBasisdx(:,:)
    LOGICAL, ALLOCATABLE :: TorqueElem(:)
    INTEGER :: i,bfid,n,nd,nbf,PrevComm,NoSlices,NoTimes
    LOGICAL :: Found, Stat
    TYPE(ValueList_t),POINTER::Params
    TYPE(GaussIntegrationPoints_t) :: IP
    TYPE(Nodes_t) :: Nodes
-   LOGICAL :: CalcTorque, CalcPot, CalcInert
+   LOGICAL :: CalcTorque, CalcPot, CalcInert, AdjointTorque 
    LOGICAL :: ThisTorque, ThisPot, ThisInert, Parallel, HaveRange
    LOGICAL :: Visited = .FALSE.
 
@@ -457,13 +457,15 @@ SUBROUTINE CalculateLumpedTransient(Torque)
 
    IF(.NOT. (CalcTorque .OR. CalcPot .OR. CalcInert) ) RETURN
 
+   AdjointTorque = ListGetLogical( Solver % Values,'Solve Adjoint Equation', Found )
+
    Parallel = ( ParEnv % PEs > 1 )
 
    nbf = Model % NumberOfBodyForces
    IF(.NOT. Visited ) THEN
      n = Model % Mesh % MaxElementDofs
      ALLOCATE( a(nbf), u(nbf), POT(n), dPOT(n), pPot(n), Density(n), &
-         Basis(n), dBasisdx(n,3) )
+         Basis(n), dBasisdx(n,3), TorqSens(n) )
    END IF
 
    IF( CalcTorque ) THEN
@@ -611,14 +613,33 @@ SUBROUTINE CalculateLumpedTransient(Torque)
          U(bfid) = U(bfid) + Weight * SUM(dPot(1:nd)*Basis(1:nd))
        END IF
 
-       IF( ThisTorque ) THEN                      
+       IF( ThisTorque ) THEN
+         wtorq = weight * r / (PI*4.0d-7*rdiff) 
+
          Bx =  SUM(POT(1:nd)*dBasisdx(1:nd,2))
          By = -SUM(POT(1:nd)*dBasisdx(1:nd,1))
          Br =  x/r*Bx + y/r*By
          Bp = -y/r*Bx + x/r*By
-         Torq = Torq + Weight * r*Br*Bp / (PI*4.0d-7*rdiff)
+
+         Torq = Torq + wtorq * Br*Bp 
          TorqArea = TorqArea + Weight
+
+         ! This is tentative addition for solving the adjoint equation. 
+         IF( AdjointTorque ) THEN
+           DO j=1,nd
+             Bx0 =  dBasisdx(j,2)
+             By0 = -dBasisdx(j,1)
+             Br0 =  x/r*Bx0 + y/r*By0
+             Bp0 = -y/r*Bx0 + x/r*By0             
+             TorqSens(j) = wtorq * (Br0*Bp+Br*Bp0)              
+           END DO
+
+           CALL UpdateGlobalForce( Solver % Matrix % RhsAdjoint, &
+               TorqSens(1:nd), nd, Solver % Variable % DOFs, &
+               Solver % Variable % Perm(Element % NodeIndexes(1:nd)) )
+         END IF
        END IF
+
 
        IF( ThisInert ) THEN
          IF( r < rmean ) THEN