From c7120f136db4a2f5c5b94f4e2a2972efe5bcce76 Mon Sep 17 00:00:00 2001 From: Peter R Date: Wed, 8 Jan 2025 16:11:28 +0200 Subject: [PATCH] Tentative steps for adjoint rhs --- fem/src/DefUtils.F90 | 55 +++++++++++++++++++-------- fem/src/Types.F90 | 2 + fem/src/modules/MagnetoDynamics2D.F90 | 33 +++++++++++++--- 3 files changed, 68 insertions(+), 22 deletions(-) diff --git a/fem/src/DefUtils.F90 b/fem/src/DefUtils.F90 index 49ccd61907..6e9425e39b 100644 --- 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 index 1be06e81e6..ef4e2232bc 100644 --- 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 index e1f7ce50eb..d8a39f26d2 100755 --- a/fem/src/modules/MagnetoDynamics2D.F90 +++ b/fem/src/modules/MagnetoDynamics2D.F90 @@ -407,9 +407,9 @@ 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 @@ -417,7 +417,7 @@ SUBROUTINE CalculateLumpedTransient(Torque) 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