1!/*****************************************************************************/ 2! * 3! * Elmer/Ice, a glaciological add-on to Elmer 4! * http://elmerice.elmerfem.org 5! * 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! ****************************************************************************** 24! * 25! * Authors: 26! * Email: 27! * Web: http://elmerice.elmerfem.org 28! * 29! * Original Date: 30! * 31! ***************************************************************************** 32! Compute the integrated gradient of the Cost function for the 33! Adjoint Inverse Problem 34! 35! Serial/Parallel(not halo?) and 2D/3D 36! 37! Need: 38! - Navier-stokes and Adjoint Problems Solutions 39! - Name of Beta variable and Grad Variable 40! - Power formulation if 'Slip Coef'=10^Beta and optimization on Beta 41! - Beta2 formulation if 'Slip Coef'=Beta^2 and optimization on Beta 42! - Lambda : the regularization coefficient (Jr=0.5 Lambda (dBeta/dx)^2) 43! 44! If DJDBeta should be set to zero for floating ice shelves the following is 45! to be used (default is not to do this): 46! - FreeSlipShelves (logical, default false) 47! - mask name (string, default GroundedMask) 48! Note that if FreeSlipShelves is true not only is DJDBeta set to zero for 49! floating ice, but also the regularisation term NodalRegb. 50! 51! ***************************************************************************** 52SUBROUTINE DJDBeta_Adjoint( Model,Solver,dt,TransientSimulation ) 53!------------------------------------------------------------------------------ 54!****************************************************************************** 55 USE DefUtils 56 IMPLICIT NONE 57!------------------------------------------------------------------------------ 58 TYPE(Solver_t) :: Solver 59 TYPE(Model_t) :: Model 60 61 REAL(KIND=dp) :: dt 62 LOGICAL :: TransientSimulation 63! 64 TYPE(ValueList_t), POINTER :: BC 65 66 CHARACTER(LEN=MAX_NAME_LEN) :: SolverName,NeumannSolName,AdjointSolName 67 CHARACTER(LEN=MAX_NAME_LEN) :: VarSolName,GradSolName 68 TYPE(Element_t),POINTER :: Element 69 TYPE(Variable_t), POINTER :: PointerToVariable, BetaVariable, VeloSolN,VeloSolD 70 TYPE(ValueList_t), POINTER :: SolverParams 71 TYPE(Nodes_t) :: ElementNodes 72 TYPE(GaussIntegrationPoints_t) :: IntegStuff 73 74 REAL(KIND=dp), POINTER :: VariableValues(:),VelocityN(:),VelocityD(:),BetaValues(:) 75 INTEGER, POINTER :: Permutation(:), VeloNPerm(:),VeloDPerm(:),BetaPerm(:),NodeIndexes(:) 76 77 real(kind=dp),allocatable :: VisitedNode(:),db(:),Basis(:),dBasisdx(:,:) 78 real(kind=dp),allocatable :: nodalbetab(:),NodalRegb(:) 79 real(kind=dp) :: betab 80 real(kind=dp) :: u,v,w,SqrtElementMetric,s 81 real(kind=dp) :: Lambda 82 REAL(KIND=dp) :: Normal(3),Tangent(3),Tangent2(3),Vect(3) 83 84 integer :: i,j,k,e,t,n,NMAX,NActiveNodes,DIM 85 integer :: p,q 86 87 logical :: PowerFormulation,Beta2Formulation 88 Logical :: Firsttime=.true.,Found,stat,UnFoundFatal=.TRUE. 89 Logical :: NormalTangential1,NormalTangential2 90 91 ! Variables for setting DJDBeta to zero for ice shelves. 92 TYPE(Variable_t), POINTER :: PointerToMask => NULL() 93 REAL(KIND=dp), POINTER :: MaskValues(:) => NULL() 94 INTEGER, POINTER :: MaskPerm(:) => NULL() 95 LOGICAL :: FreeSlipShelves 96 CHARACTER(LEN=MAX_NAME_LEN) :: MaskName 97 98 save FreeSlipShelves,MaskName 99 save SolverName,NeumannSolName,AdjointSolName,VarSolName,GradSolName 100 save VisitedNode,db,Basis,dBasisdx,nodalbetab,NodalRegb 101 save Firsttime,DIM,Lambda 102 save ElementNodes 103 save PowerFormulation,Beta2Formulation 104 105 If (Firsttime) then 106 107 DIM = CoordinateSystemDimension() 108 WRITE(SolverName, '(A)') 'DJDBeta_Adjoint' 109 110 CALL Info(SolverName,'***********************',level=0) 111 CALL Info(SolverName,' This solver has been replaced by:',level=0) 112 CALL Info(SolverName,' AdjointStokes_GradientBetaSolver ',level=0) 113 CALL Info(SolverName,' See documentation under: ',level=0) 114 CALL Info(SolverName,' elmerice/Solvers/Documentation ',level=0) 115 CALL Info(SolverName,'***********************',level=0) 116 CALL FATAL(SolverName,' Use new solver !!') 117 118 NMAX=Solver % Mesh % NumberOfNodes 119 allocate(VisitedNode(NMAX),db(NMAX), & 120 nodalbetab(Model % MaxElementNodes),& 121 NodalRegb(Model % MaxElementNodes),& 122 Basis(Model % MaxElementNodes), & 123 dBasisdx(Model % MaxElementNodes,3)) 124 125!!!!!!!!!!! get Solver Variables 126 SolverParams => GetSolverParams() 127 128 NeumannSolName = GetString( SolverParams,'Flow Solution Name', Found) 129 IF(.NOT.Found) THEN 130 CALL WARN(SolverName,'Keyword >Neumann Solution Name< not found in section >Solver<') 131 CALL WARN(SolverName,'Taking default value >Flow Solution<') 132 WRITE(NeumannSolName,'(A)') 'Flow Solution' 133 END IF 134 AdjointSolName = GetString( SolverParams,'Adjoint Solution Name', Found) 135 IF(.NOT.Found) THEN 136 CALL WARN(SolverName,'Keyword >Adjoint Solution Name< not found in section >Solver<') 137 CALL WARN(SolverName,'Taking default value >Adjoint<') 138 WRITE(AdjointSolName,'(A)') 'Adjoint' 139 END IF 140 VarSolName = GetString( SolverParams,'Optimized Variable Name', Found) 141 IF(.NOT.Found) THEN 142 CALL WARN(SolverName,'Keyword >Optimized Variable Name< not found in section >Solver<') 143 CALL WARN(SolverName,'Taking default value >Beta<') 144 WRITE(VarSolName,'(A)') 'Beta' 145 END IF 146 GradSolName = GetString( SolverParams,'Gradient Variable Name', Found) 147 IF(.NOT.Found) THEN 148 CALL WARN(SolverName,'Keyword >Gradient Variable Name< not found in section >Solver<') 149 CALL WARN(SolverName,'Taking default value >DJDB<') 150 WRITE(GradSolName,'(A)') 'DJDB' 151 END IF 152 PowerFormulation=GetLogical( SolverParams, 'PowerFormulation', Found) 153 IF(.NOT.Found) THEN 154 CALL WARN(SolverName,'Keyword >PowerFormulation< not found in section >Equation<') 155 CALL WARN(SolverName,'Taking default value >FALSE<') 156 PowerFormulation=.FALSE. 157 END IF 158 Beta2Formulation=GetLogical( SolverParams, 'Beta2Formulation', Found) 159 IF(.NOT.Found) THEN 160 CALL WARN(SolverName,'Keyword >Beta2Formulation< not found in section >Equation<') 161 CALL WARN(SolverName,'Taking default value >FALSE<') 162 Beta2Formulation=.FALSE. 163 END IF 164 IF (PowerFormulation.and.Beta2Formulation) then 165 WRITE(Message,'(A)') 'Can t be PowerFormulation and Beta2Formulation in the same time' 166 CALL FATAL(SolverName,Message) 167 End if 168 Lambda = GetConstReal( SolverParams,'Lambda', Found) 169 IF(.NOT.Found) THEN 170 CALL WARN(SolverName,'Keyword >Lambda< not found in section >Solver<') 171 CALL WARN(SolverName,'Taking default value Lambda=0.0') 172 Lambda = 0.0 173 End if 174 175 FreeSlipShelves=GetLogical( SolverParams, 'FreeSlipShelves', Found) 176 IF(.NOT.Found) THEN 177 CALL WARN(SolverName,'Keyword >FreeSlipShelves< not found in solver params') 178 CALL WARN(SolverName,'Taking default value >FALSE<') 179 FreeSlipShelves=.FALSE. 180 END IF 181 IF (FreeSlipShelves) THEN 182 MaskName = GetString( SolverParams,'mask name', Found) 183 IF(.NOT.Found) THEN 184 CALL WARN(SolverName,'Keyword >mask name< not found in solver section') 185 CALL WARN(SolverName,'Taking default value >GroundedMask<') 186 WRITE(MaskName,'(A)') 'GroundedMask' 187 END IF 188 END IF 189 190!!! End of First visit 191 Firsttime=.false. 192 Endif 193 194 PointerToVariable => VariableGet( Model % Mesh % Variables, GradSolName,UnFoundFatal=UnFoundFatal) 195 VariableValues => PointerToVariable % Values 196 Permutation => PointerToVariable % Perm 197 VariableValues=0._dp 198 199 BetaVariable => VariableGet( Model % Mesh % Variables, VarSolName,UnFoundFatal=UnFoundFatal) 200 BetaValues => BetaVariable % Values 201 BetaPerm => BetaVariable % Perm 202 203 VeloSolN => VariableGet( Model % Mesh % Variables, NeumannSolName,UnFoundFatal=UnFoundFatal) 204 VelocityN => VeloSolN % Values 205 VeloNPerm => VeloSolN % Perm 206 207 VeloSolD => VariableGet( Model % Mesh % Variables, AdjointSolName,UnFoundFatal=UnFoundFatal) 208 VelocityD => VeloSolD % Values 209 VeloDPerm => VeloSolD % Perm 210 211 IF (FreeSlipShelves) THEN 212 PointerToMask => VariableGet( Model % Variables, MaskName, UnFoundFatal=.TRUE.) 213 MaskValues => PointerToMask % Values 214 MaskPerm => PointerToMask % Perm 215 END IF 216 217 VisitedNode=0.0_dp 218 db=0.0_dp 219 220 DO e=1,Solver % NumberOfActiveElements 221 Element => GetActiveElement(e) 222 CALL GetElementNodes( ElementNodes ) 223 n = GetElementNOFNodes() 224 NodeIndexes => Element % NodeIndexes 225 226 ! Compute Nodal Value of DJDBeta 227 BC => GetBC() 228 if (.NOT.ASSOCIATED(BC)) CYCLE 229 230 NormalTangential1 = GetLogical( BC, & 231 'Normal-Tangential Velocity', Found ) 232 IF (.NOT.Found) then 233 NormalTangential1 = GetLogical( BC, & 234 'Normal-Tangential '//trim(NeumannSolName), Found) 235 END IF 236 NormalTangential2 = GetLogical( BC, & 237 'Normal-Tangential '//trim(AdjointSolName), Found) 238 IF (NormalTangential1.NEQV.NormalTangential2) then 239 WRITE(Message,'(A,I1,A,I1)') & 240 'NormalTangential Velocity is : ',NormalTangential1, & 241 'But NormalTangential Adjoint is : ',NormalTangential2 242 CALL FATAL(SolverName,Message) 243 ENDIF 244 IF (.NOT.NormalTangential1) then 245 WRITE(Message,'(A)') & 246 'ALWAYS USE Normal-Tangential COORDINATES with SlipCoef 2=SlipCoef 3' 247 CALL FATAL(SolverName,Message) 248 ENDIF 249 250 VisitedNode(NodeIndexes(1:n))=VisitedNode(NodeIndexes(1:n))+1.0_dp 251 252 ! Compute Integrated Nodal Value of DJDBeta 253 nodalbetab=0.0_dp 254 NodalRegb=0.0_dp 255 256 257 IntegStuff = GaussPoints( Element ) 258 DO t=1,IntegStuff % n 259 U = IntegStuff % u(t) 260 V = IntegStuff % v(t) 261 W = IntegStuff % w(t) 262 stat = ElementInfo(Element,ElementNodes,U,V,W,SqrtElementMetric, & 263 Basis,dBasisdx ) 264 265 s = SqrtElementMetric * IntegStuff % s(t) 266 267 ! compute gradient from Stokes and adjoint computation 268 ! follow the computation of the stiffMatrix as done in the NS solver 269 Normal = NormalVector( Element, ElementNodes, u,v,.TRUE. ) 270 SELECT CASE( Element % TYPE % DIMENSION ) 271 CASE(1) 272 Tangent(1) = Normal(2) 273 Tangent(2) = -Normal(1) 274 Tangent(3) = 0.0_dp 275 Tangent2 = 0.0_dp 276 CASE(2) 277 CALL TangentDirections( Normal, Tangent, Tangent2 ) 278 END SELECT 279 280 281 betab=0.0_dp 282 Do p=1,n 283 Do q=1,n 284 285 Do i=2,dim 286 SELECT CASE(i) 287 CASE(2) 288 Vect = Tangent 289 CASE(3) 290 Vect = Tangent2 291 END SELECT 292 293 Do j=1,DIM 294 Do k=1,DIM 295 betab = betab + s * Basis(q) * Basis(p) * Vect(j) * Vect(k) * & 296 (- VelocityN((DIM+1)*(VeloNPerm(NodeIndexes(q))-1)+k) * & 297 VelocityD((DIM+1)*(VeloDPerm(NodeIndexes(p))-1)+j)) 298 End Do !on k 299 End Do !on j 300 End Do !on i 301 End Do !on q 302 End Do !on p 303 304 nodalbetab(1:n)=nodalbetab(1:n)+betab*Basis(1:n) 305 306 If (Lambda /= 0.0) then 307 NodalRegb(1:n)=NodalRegb(1:n)+& 308 s*Lambda*(SUM(dBasisdx(1:n,1)*BetaValues(BetaPerm(NodeIndexes(1:n))))*dBasisdx(1:n,1)) 309 IF (DIM.eq.3) then 310 NodalRegb(1:n)=NodalRegb(1:n)+& 311 s*Lambda*(SUM(dBasisdx(1:n,2)*BetaValues(BetaPerm(NodeIndexes(1:n))))*dBasisdx(1:n,2)) 312 End if 313 End if 314 End DO ! on IPs 315 316 IF (PowerFormulation) then 317 nodalbetab(1:n)=nodalbetab(1:n)*(10**(BetaValues(BetaPerm(NodeIndexes(1:n)))))*log(10.0) 318 ENDIF 319 IF (Beta2Formulation) then 320 nodalbetab(1:n)=nodalbetab(1:n)*2.0_dp*BetaValues(BetaPerm(NodeIndexes(1:n))) 321 END IF 322 323 ! Set regularisation to zero for floating points 324 IF (FreeSlipShelves) THEN 325 DO t=1,n 326 IF ( MaskValues(MaskPerm(NodeIndexes(t))).LT.0.0_dp ) THEN 327 NodalRegb(t) = 0.0_dp 328 END IF 329 END DO 330 END IF 331 332 db(NodeIndexes(1:n)) = db(NodeIndexes(1:n)) + nodalbetab(1:n) + NodalRegb(1:n) 333 End do ! on elements 334 335 Do t=1,Solver % Mesh % NumberOfNodes 336 if (VisitedNode(t).lt.1.0_dp) cycle 337 VariableValues(Permutation(t)) = db(t) 338 IF (FreeSlipShelves) THEN 339 IF ( MaskValues(MaskPerm(t)).LT.0.0_dp ) THEN 340 VariableValues(Permutation(t)) = 0.0_dp 341 END IF 342 END IF 343 End do 344 345 IF (FreeSlipShelves) THEN 346 NULLIFY(PointerToMask) 347 NULLIFY(MaskValues) 348 NULLIFY(MaskPerm) 349 END IF 350 351 Return 352 353CONTAINS 354 355 function calcNorm(v) result(v2) 356 implicit none 357 real(kind=dp) :: v(3),v2 358 359 v2=v(1)*v(1)+v(2)*v(2)+v(3)*v(3) 360 end function calcNorm 361 362 function scalar(v1,v2) result(vr) 363 implicit none 364 real(kind=dp) :: v2(3),v1(3),vr 365 366 vr=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) 367 end function scalar 368 369!------------------------------------------------------------------------------ 370END SUBROUTINE DJDBeta_Adjoint 371!------------------------------------------------------------------------------ 372 373 374