![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 26478 for issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
- Timestamp:
- 10/13/21 07:56:16 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp ¶
r26475 r26478 52 52 case FSApproximationEnum: 53 53 return CreateKMatrixFS(element); 54 case MLHOApproximationEnum: 55 // a more accurate option, but integrate in the vertical direction numerically. 56 // return CreateKMatrixMLHOVerticalIntergrated(element); 57 return CreateKMatrixMLHO(element); 54 58 case NoneApproximationEnum: 55 59 return NULL; … … 215 219 xDelete<IssmDouble>(xyz_list); 216 220 return Ke; 221 }/*}}}*/ 222 ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/ 223 224 /* Check if ice in element */ 225 if(!element->IsIceInElement()) return NULL; 226 227 /*Intermediaries */ 228 bool incomplete_adjoint; 229 IssmDouble Jdet,mu_prime,n,thickness,mu,effmu; 230 IssmDouble *xyz_list = NULL; 231 IssmDouble viscosity[9]; //9 mu for different integrand 232 int domaintype; 233 int dim=2; 234 235 IssmDouble eb1i,eb1j,esh1i,esh1j,eb2i,eb2j,esh2i,esh2j; 236 IssmDouble epsilon[5],epsilonbase[5],epsilonshear[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ 237 IssmDouble e1b[2], e2b[2], e1sh[2], e2sh[2]; 238 IssmDouble vxshear, vyshear; 239 240 Element* basalelement; 241 242 /*Get basal element*/ 243 element->FindParam(&domaintype,DomainTypeEnum); 244 switch(domaintype){ 245 case Domain2DhorizontalEnum: 246 basalelement = element; 247 break; 248 case Domain3DEnum: case Domain2DverticalEnum: 249 _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 250 break; 251 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 252 } 253 /*Fetch number of nodes and dof for this finite element*/ 254 int numnodes = basalelement->GetNumberOfNodes(); 255 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes); // like SSA 256 IssmDouble* basis = xNew<IssmDouble>(numnodes); // like SSA 257 258 /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ 259 element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); 260 StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); 261 ElementMatrix* Ke=analysis->CreateKMatrix(basalelement); 262 delete analysis; 263 if(incomplete_adjoint) return Ke; 264 265 /*Retrieve all inputs and parameters*/ 266 element->GetVerticesCoordinates(&xyz_list); 267 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); //vertically integrated vx 268 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); //vertically integrated vy 269 Input* vxbase_input = element->GetInput(VxBaseEnum); _assert_(vxbase_input); 270 Input* vybase_input = element->GetInput(VyBaseEnum); _assert_(vybase_input); 271 Input* vxshear_input = element->GetInput(VxShearEnum); _assert_(vxshear_input); 272 Input* vyshear_input = element->GetInput(VyShearEnum); _assert_(vyshear_input); 273 Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); 274 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); 275 276 /* Start looping on the number of gaussian points: */ 277 Gauss* gauss = element->NewGauss(5); 278 Gauss* gauss_base = basalelement->NewGauss(); 279 while(gauss->next()){ 280 gauss->SynchronizeGaussBase(gauss_base); 281 282 element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); 283 basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); 284 basalelement->NodalFunctions(basis, gauss_base); 285 286 thickness_input->GetInputValue(&thickness, gauss); 287 n_input->GetInputValue(&n,gauss); 288 289 vxshear_input->GetInputValue(&vxshear,gauss); 290 vyshear_input->GetInputValue(&vyshear,gauss); 291 292 element->material->ViscosityMLHOAdjoint(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input); 293 294 effmu = 2.0*(1-n)/2.0/n; 295 296 element->StrainRateHO(&epsilonbase[0],xyz_list,gauss,vxbase_input, vybase_input); 297 element->StrainRateHO(&epsilonshear[0],xyz_list,gauss,vxshear_input, vyshear_input); 298 299 e1b[0] = 2*epsilonbase[0]+epsilonbase[1]; e1b[1] = epsilonbase[2]; 300 e2b[1] = epsilonbase[0]+2*epsilonbase[1]; e2b[0] = epsilonbase[2]; 301 302 e1sh[0] = 2*epsilonshear[0]+epsilonshear[1]; e1sh[1] = epsilonshear[2]; 303 e2sh[1] = epsilonshear[0]+2*epsilonshear[1]; e2sh[0] = epsilonshear[2]; 304 305 for(int i=0;i<numnodes;i++){ 306 for(int j=0;j<numnodes;j++){ 307 eb1i = e1b[0]*dbasis[0*numnodes+i]+e1b[1]*dbasis[1*numnodes+i]; 308 eb1j = e1b[0]*dbasis[0*numnodes+j]+e1b[1]*dbasis[1*numnodes+j]; 309 eb2i = e2b[0]*dbasis[0*numnodes+i]+e2b[1]*dbasis[1*numnodes+i]; 310 eb2j = e2b[0]*dbasis[0*numnodes+j]+e2b[1]*dbasis[1*numnodes+j]; 311 esh1i = e1sh[0]*dbasis[0*numnodes+i]+e1sh[1]*dbasis[1*numnodes+i]; 312 esh1j = e1sh[0]*dbasis[0*numnodes+j]+e1sh[1]*dbasis[1*numnodes+j]; 313 esh2i = e2sh[0]*dbasis[0*numnodes+i]+e2sh[1]*dbasis[1*numnodes+i]; 314 esh2j = e2sh[0]*dbasis[0*numnodes+j]+e2sh[1]*dbasis[1*numnodes+j]; 315 316 Ke->values[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb1i+viscosity[1]*(esh1j*eb1i+eb1j*esh1i)+viscosity[2]*esh1j*esh1i); 317 Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1j*basis[i]); 318 Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb2i+viscosity[1]*(esh1j*eb2i+eb1j*esh2i)+viscosity[2]*esh1j*esh2i); 319 Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1j*basis[i]); 320 321 Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1i*basis[j]); 322 Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb1i+viscosity[4]*(esh1j*eb1i+eb1j*esh1i)+viscosity[5]*esh1j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(eb1i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(esh1i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); 323 Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2i*basis[j]); 324 Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb2i+viscosity[4]*(esh1j*eb2i+eb1j*esh2i)+viscosity[5]*esh1j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*eb2i*basis[j]+vyshear*eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2i*basis[j]+vyshear*esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vyshear*basis[j]*basis[i]); 325 326 327 Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb1i+viscosity[1]*(esh2j*eb1i+eb2j*esh1i)+viscosity[2]*esh2j*esh1i); 328 Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2j*basis[i]); 329 Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb2i+viscosity[1]*(esh2j*eb2i+eb2j*esh2i)+viscosity[2]*esh2j*esh2i); 330 Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2j*basis[i]); 331 332 333 Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1i*basis[j]); 334 Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb1i+viscosity[4]*(esh2j*eb1i+eb2j*esh1i)+viscosity[5]*esh2j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vyshear*eb1i*basis[j]+vxshear*eb2j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2j*basis[i]+vyshear*esh1i*basis[j])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); 335 Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2i*basis[j]); 336 Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb2i+viscosity[4]*(esh2j*eb2i+eb2j*esh2i)+viscosity[5]*esh2j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(eb2i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(esh2i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vyshear*vyshear*basis[j]*basis[i]); 337 } 338 } 339 } 340 341 /*Transform Coordinate System*/ 342 /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/ 343 344 /*Clean up and return*/ 345 delete gauss; 346 delete gauss_base; 347 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; 348 xDelete<IssmDouble>(xyz_list); 349 xDelete<IssmDouble>(dbasis); 350 xDelete<IssmDouble>(basis); 351 return Ke; 352 }/*}}}*/ 353 ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHOVerticalIntergrated(Element* element){/*{{{*/ 354 355 /* Check if ice in element */ 356 if(!element->IsIceInElement()) return NULL; 357 358 /*Intermediaries */ 359 bool incomplete_adjoint; 360 IssmDouble Jdet,mu_prime,n,thickness,mu,effmu; 361 IssmDouble *xyz_list = NULL; 362 IssmDouble zeta, epsilon_eff; 363 int domaintype; 364 int dim=2; 365 366 IssmDouble e1phi1i, e1phi1j, e2phi1i, e2phi1j, e1phi2i, e1phi2j, e2phi2i, e2phi2j; 367 IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ 368 IssmDouble e1[3],e2[3]; 369 370 Element* basalelement; 371 372 /*Get basal element*/ 373 element->FindParam(&domaintype,DomainTypeEnum); 374 switch(domaintype){ 375 case Domain2DhorizontalEnum: 376 basalelement = element; 377 break; 378 case Domain3DEnum: case Domain2DverticalEnum: 379 _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 380 break; 381 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 382 } 383 /*Fetch number of nodes and dof for this finite element*/ 384 int numnodes = basalelement->GetNumberOfNodes(); 385 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes); // like SSA 386 IssmDouble* basis = xNew<IssmDouble>(numnodes); // like SSA 387 388 /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ 389 element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); 390 StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); 391 ElementMatrix* Ke=analysis->CreateKMatrix(basalelement); 392 delete analysis; 393 if(incomplete_adjoint) return Ke; 394 395 /*Retrieve all inputs and parameters*/ 396 element->GetVerticesCoordinates(&xyz_list); 397 Input* vxbase_input = element->GetInput(VxBaseEnum); _assert_(vxbase_input); 398 Input* vybase_input = element->GetInput(VyBaseEnum); _assert_(vybase_input); 399 Input* vxshear_input = element->GetInput(VxShearEnum); _assert_(vxshear_input); 400 Input* vyshear_input = element->GetInput(VyShearEnum); _assert_(vyshear_input); 401 Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); 402 Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); 403 404 /* Start looping on the number of gaussian points: */ 405 Gauss* gauss = element->NewGauss(10); 406 Gauss* gauss_base = basalelement->NewGauss(); 407 408 GaussSeg* gauss_seg=new GaussSeg(5); 409 410 while(gauss->next()){ 411 gauss->SynchronizeGaussBase(gauss_base); 412 413 element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); 414 basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); 415 basalelement->NodalFunctions(basis, gauss_base); 416 417 thickness_input->GetInputValue(&thickness, gauss); 418 n_input->GetInputValue(&n,gauss); 419 420 effmu = 2*(1-n)/2.0/n; 421 422 /* Get the integration in the vertical direction */ 423 gauss_seg->Reset(); 424 while(gauss_seg->next()){ 425 /*Compute zeta for gauss_seg point (0=surface, 1=base)*/ 426 zeta=0.5*(gauss_seg->coord1+1); 427 428 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy (for a given zeta)*/ 429 element->StrainRateMLHO(&epsilon[0],xyz_list,gauss, 430 vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); 431 epsilon_eff=sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] 432 + epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]); 433 434 /*Get viscosity at zeta*/ 435 element->material->GetViscosity(&mu,epsilon_eff,gauss); 436 /*the adjoint viscosity with zeta dependent term*/ 437 mu = mu /epsilon_eff/epsilon_eff; 438 439 e1[0] = 2.0*epsilon[0]+epsilon[1]; 440 e1[1] = epsilon[2]; 441 e1[2] = epsilon[3]; 442 443 e2[0] = epsilon[2]; 444 e2[1] = epsilon[0]+2.0*epsilon[1]; 445 e2[2] = epsilon[4]; 446 447 for(int i=0;i<numnodes;i++){ 448 for(int j=0;j<numnodes;j++){ 449 e1phi1i = e1[0]*dbasis[0*numnodes+i]+e1[1]*dbasis[1*numnodes+i]; 450 e2phi1i = e2[0]*dbasis[0*numnodes+i]+e2[1]*dbasis[1*numnodes+i]; 451 e1phi1j = e1[0]*dbasis[0*numnodes+j]+e1[1]*dbasis[1*numnodes+j]; 452 e2phi1j = e2[0]*dbasis[0*numnodes+j]+e2[1]*dbasis[1*numnodes+j]; 453 454 e1phi2i = (1-pow(zeta, n+1))*(e1[0]*dbasis[0*numnodes+i]+e1[1]*dbasis[1*numnodes+i])+(n+1)/thickness*pow(zeta, n)*e1[2]*basis[i]; 455 e2phi2i = (1-pow(zeta, n+1))*(e2[0]*dbasis[0*numnodes+i]+e2[1]*dbasis[1*numnodes+i])+(n+1)/thickness*pow(zeta, n)*e2[2]*basis[i]; 456 e1phi2j = (1-pow(zeta, n+1))*(e1[0]*dbasis[0*numnodes+j]+e1[1]*dbasis[1*numnodes+j])+(n+1)/thickness*pow(zeta, n)*e1[2]*basis[j]; 457 e2phi2j = (1-pow(zeta, n+1))*(e2[0]*dbasis[0*numnodes+j]+e2[1]*dbasis[1*numnodes+j])+(n+1)/thickness*pow(zeta, n)*e2[2]*basis[j]; 458 459 Ke->values[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e1phi1j*thickness*0.5*gauss_seg->weight; 460 Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e1phi1j*thickness*0.5*gauss_seg->weight; 461 Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e1phi1j*thickness*0.5*gauss_seg->weight; 462 Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e1phi1j*thickness*0.5*gauss_seg->weight; 463 464 Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e1phi2j*thickness*0.5*gauss_seg->weight; 465 Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e1phi2j*thickness*0.5*gauss_seg->weight; 466 Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e1phi2j*thickness*0.5*gauss_seg->weight; 467 Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e1phi2j*thickness*0.5*gauss_seg->weight; 468 469 Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e2phi1j*thickness*0.5*gauss_seg->weight; 470 Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e2phi1j*thickness*0.5*gauss_seg->weight; 471 Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e2phi1j*thickness*0.5*gauss_seg->weight; 472 Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e2phi1j*thickness*0.5*gauss_seg->weight; 473 474 Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e2phi2j*thickness*0.5*gauss_seg->weight; 475 Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e2phi2j*thickness*0.5*gauss_seg->weight; 476 Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e2phi2j*thickness*0.5*gauss_seg->weight; 477 Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e2phi2j*thickness*0.5*gauss_seg->weight; 478 } 479 } 480 } 481 } 482 /*Transform Coordinate System*/ 483 /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/ 484 485 /*Clean up and return*/ 486 delete gauss; 487 delete gauss_base; 488 delete gauss_seg; 489 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; 490 xDelete<IssmDouble>(xyz_list); 491 xDelete<IssmDouble>(dbasis); 492 xDelete<IssmDouble>(basis); 493 return Ke; 217 494 }/*}}}*/ 218 495 ElementMatrix* AdjointHorizAnalysis::CreateKMatrixL1L2(Element* element){/*{{{*/ … … 335 612 case FSApproximationEnum: 336 613 return CreatePVectorFS(element); 614 case MLHOApproximationEnum: 615 return CreatePVectorMLHO(element); 337 616 case NoneApproximationEnum: 338 617 return NULL; … … 820 1099 821 1100 }/*}}}*/ 1101 ElementVector* AdjointHorizAnalysis::CreatePVectorMLHO(Element* element){/*{{{*/ 1102 1103 /*Intermediaries*/ 1104 int domaintype; 1105 Element* basalelement; 1106 1107 /* Check if ice in element */ 1108 if(!element->IsIceInElement()) return NULL; 1109 1110 /*Get basal element*/ 1111 element->FindParam(&domaintype,DomainTypeEnum); 1112 switch(domaintype){ 1113 case Domain2DhorizontalEnum: 1114 basalelement = element; 1115 break; 1116 case Domain2DverticalEnum: 1117 if(!element->IsOnBase()) return NULL; 1118 basalelement = element->SpawnBasalElement(); 1119 break; 1120 case Domain3DEnum: 1121 if(!element->IsOnBase()) return NULL; 1122 basalelement = element->SpawnBasalElement(true); 1123 break; 1124 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 1125 } 1126 1127 /*Intermediaries */ 1128 int num_responses,i; 1129 IssmDouble Jdet,obs_velocity_mag,velocity_mag; 1130 IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; 1131 IssmDouble scalex,scaley,scale,S,thickness,n; 1132 int *responses = NULL; 1133 IssmDouble *xyz_list = NULL; 1134 1135 /*Fetch number of nodes and dof for this finite element*/ 1136 int numnodes = basalelement->GetNumberOfNodes(); 1137 1138 /*Initialize Element vector and vectors*/ 1139 ElementVector* pe = basalelement->NewElementVector(MLHOApproximationEnum); 1140 IssmDouble* basis = xNew<IssmDouble>(numnodes); 1141 1142 /*Retrieve all inputs and parameters*/ 1143 basalelement->GetVerticesCoordinates(&xyz_list); 1144 basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum); 1145 basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum); 1146 DatasetInput* weights_input = basalelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 1147 Input* vx_input = basalelement->GetInput(VxSurfaceEnum); _assert_(vx_input); 1148 Input* vxobs_input = basalelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 1149 Input* thickness_input=basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); 1150 Input* n_input =element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); 1151 Input* vy_input=NULL; 1152 Input* vyobs_input=NULL; 1153 if(domaintype!=Domain2DverticalEnum){ 1154 vy_input = basalelement->GetInput(VySurfaceEnum); _assert_(vy_input); 1155 vyobs_input = basalelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 1156 } 1157 IssmDouble epsvel = 2.220446049250313e-16; 1158 IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ 1159 1160 /*Get Surface if required by one response*/ 1161 Input* S_input = NULL; 1162 for(int resp=0;resp<num_responses;resp++){ 1163 if(responses[resp]==SurfaceAverageVelMisfitEnum){ 1164 S_input = element->GetInput(SurfaceAreaEnum); _assert_(S_input); break; 1165 } 1166 } 1167 1168 /* Start looping on the number of gaussian points: */ 1169 Gauss* gauss=basalelement->NewGauss(4); 1170 while(gauss->next()){ 1171 1172 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); 1173 basalelement->NodalFunctions(basis, gauss); 1174 1175 thickness_input->GetInputValue(&thickness,gauss); 1176 n_input->GetInputValue(&n,gauss); 1177 vx_input->GetInputValue(&vx,gauss); 1178 vxobs_input->GetInputValue(&vxobs,gauss); 1179 if(domaintype!=Domain2DverticalEnum){ 1180 vy_input->GetInputValue(&vy,gauss); 1181 vyobs_input->GetInputValue(&vyobs,gauss); 1182 } 1183 /*Loop over all requested responses*/ 1184 for(int resp=0;resp<num_responses;resp++){ 1185 weights_input->GetInputValue(&weight,gauss,responses[resp]); 1186 1187 switch(responses[resp]){ 1188 case SurfaceAbsVelMisfitEnum: 1189 /* 1190 * 1 [ 2 2 ] 1191 * J = --- | (u - u ) + (v - v ) | 1192 * 2 [ obs obs ] 1193 * 1194 * dJ 1195 * DU = - -- = (u - u ) 1196 * du obs 1197 */ 1198 for(i=0;i<numnodes;i++){ 1199 if(domaintype!=Domain2DverticalEnum){ 1200 dux=vxobs-vx; 1201 duy=vyobs-vy; 1202 pe->values[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 1203 pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; 1204 pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; 1205 pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; 1206 } 1207 else { 1208 _error_("2D vertical is not implemented for MLHO"); 1209 } 1210 } 1211 break; 1212 case SurfaceRelVelMisfitEnum: 1213 /* 1214 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] 1215 * J = --- | ------------- (u - u ) + ------------- (v - v ) | 1216 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] 1217 * obs obs 1218 * 1219 * dJ \bar{v}^2 1220 * DU = - -- = ------------- (u - u ) 1221 * du (u + eps)^2 obs 1222 * obs 1223 */ 1224 for(i=0;i<numnodes;i++){ 1225 if(domaintype!=Domain2DverticalEnum){ 1226 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 1227 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 1228 dux=scalex*(vxobs-vx); 1229 duy=scaley*(vyobs-vy); 1230 pe->values[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 1231 pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; 1232 pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; 1233 pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; 1234 } 1235 else{ 1236 _error_("2D vertical is not implemented for MLHO"); 1237 } 1238 } 1239 break; 1240 case SurfaceLogVelMisfitEnum: 1241 /* 1242 * [ vel + eps ] 2 1243 * J = 4 \bar{v}^2 | log ( ----------- ) | 1244 * [ vel + eps ] 1245 * obs 1246 * 1247 * dJ 2 * log(...) 1248 * DU = - -- = - 4 \bar{v}^2 ------------- u 1249 * du vel^2 + eps 1250 * 1251 */ 1252 for(i=0;i<numnodes;i++){ 1253 if(domaintype!=Domain2DverticalEnum){ 1254 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 1255 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 1256 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 1257 dux=scale*vx; 1258 duy=scale*vy; 1259 pe->values[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 1260 pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; 1261 pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; 1262 pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; 1263 } 1264 else{ 1265 _error_("2D vertical is not implemented for MLHO"); 1266 } 1267 } 1268 break; 1269 case SurfaceAverageVelMisfitEnum: 1270 /* 1271 * 1 2 2 1272 * J = --- sqrt( (u - u ) + (v - v ) ) 1273 * S obs obs 1274 * 1275 * dJ 1 1 1276 * DU = - -- = - --- ----------- * 2 (u - u ) 1277 * du S 2 sqrt(...) obs 1278 */ 1279 S_input->GetInputValue(&S,gauss); 1280 for(i=0;i<numnodes;i++){ 1281 if(domaintype!=Domain2DverticalEnum){ 1282 scale=1./(S*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 1283 dux=scale*(vxobs-vx); 1284 duy=scale*(vyobs-vy); 1285 pe->values[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 1286 pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; 1287 pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; 1288 pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; 1289 } 1290 else{ 1291 _error_("2D vertical is not implemented for MLHO"); 1292 } 1293 } 1294 break; 1295 case SurfaceLogVxVyMisfitEnum: 1296 /* 1297 * 1 [ |u| + eps 2 |v| + eps 2 ] 1298 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | 1299 * 2 [ |u |+ eps |v |+ eps ] 1300 * obs obs 1301 * dJ 1 u 1 1302 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ 1303 * du |u| + eps |u| u + eps 1304 */ 1305 for(i=0;i<numnodes;i++){ 1306 if(domaintype!=Domain2DverticalEnum){ 1307 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1308 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1309 pe->values[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 1310 pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; 1311 pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; 1312 pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; 1313 } 1314 else{ 1315 _error_("2D vertical is not implemented for MLHO"); 1316 } 1317 } 1318 break; 1319 case DragCoefficientAbsGradientEnum: 1320 /*Nothing in P vector*/ 1321 break; 1322 case ThicknessAbsGradientEnum: 1323 /*Nothing in P vector*/ 1324 break; 1325 case ThicknessAlongGradientEnum: 1326 /*Nothing in P vector*/ 1327 break; 1328 case ThicknessAcrossGradientEnum: 1329 /*Nothing in P vector*/ 1330 break; 1331 case RheologyBbarAbsGradientEnum: 1332 /*Nothing in P vector*/ 1333 break; 1334 case RheologyBAbsGradientEnum: 1335 /*Nothing in P vector*/ 1336 break; 1337 case RheologyBInitialguessMisfitEnum: 1338 /*Nothing in P vector*/ 1339 break; 1340 default: 1341 _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); 1342 } 1343 } 1344 } 1345 1346 /*Transform coordinate system*/ 1347 //if(domaintype!=Domain2DverticalEnum) basalelement->TransformLoadVectorCoord(pe,XYEnum); 1348 /*Clean up and return*/ 1349 xDelete<int>(responses); 1350 xDelete<IssmDouble>(xyz_list); 1351 xDelete<IssmDouble>(basis); 1352 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; 1353 delete gauss; 1354 return pe; 1355 }/*}}}*/ 822 1356 ElementVector* AdjointHorizAnalysis::CreatePVectorSSA(Element* element){/*{{{*/ 823 1357 … … 1129 1663 case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_interp,control_index); break; 1130 1664 case HOApproximationEnum: GradientJDragHO( element,gradient,control_interp,control_index); break; 1665 case MLHOApproximationEnum: GradientJDragMLHO( element,gradient,control_interp,control_index); break; 1131 1666 case FSApproximationEnum: GradientJDragFS( element,gradient,control_interp,control_index); break; 1132 1667 case NoneApproximationEnum: /*Gradient is 0*/ break; … … 1149 1684 case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_interp,control_index); break; 1150 1685 case HOApproximationEnum: GradientJBbarHO( element,gradient,control_interp,control_index); break; 1686 case MLHOApproximationEnum: GradientJBbarMLHO( element,gradient,control_interp,control_index); break; 1151 1687 case FSApproximationEnum: GradientJBbarFS( element,gradient,control_interp,control_index); break; 1152 1688 case NoneApproximationEnum: /*Gradient is 0*/ break; … … 1158 1694 case SSAApproximationEnum: GradientJBSSA(element,gradient,control_interp,control_index); break; 1159 1695 case HOApproximationEnum: GradientJBHO( element,gradient,control_interp,control_index); break; 1696 // case MLHOApproximationEnum: GradientJBMLHO( element,gradient,control_interp,control_index); break; 1160 1697 case FSApproximationEnum: GradientJBFS( element,gradient,control_interp,control_index); break; 1161 1698 case NoneApproximationEnum: /*Gradient is 0*/ break; … … 1266 1803 /*WARNING: We use SSA as an estimate for now*/ 1267 1804 this->GradientJBbarSSA(element,gradient,control_interp,control_index); 1805 }/*}}}*/ 1806 void AdjointHorizAnalysis::GradientJBbarMLHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/ 1807 1808 if(control_interp!=P1Enum) _error_("not implemented yet..."); 1809 /*Intermediaries*/ 1810 int domaintype,dim; 1811 Element* basalelement; 1812 1813 /*Get basal element*/ 1814 element->FindParam(&domaintype,DomainTypeEnum); 1815 switch(domaintype){ 1816 case Domain2DhorizontalEnum: 1817 basalelement = element; 1818 dim = 2; 1819 break; 1820 case Domain2DverticalEnum: 1821 if(!element->IsOnBase()) return; 1822 basalelement = element->SpawnBasalElement(); 1823 dim = 1; 1824 break; 1825 case Domain3DEnum: 1826 if(!element->IsOnBase()) return; 1827 basalelement = element->SpawnBasalElement(true); 1828 dim = 2; 1829 break; 1830 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 1831 } 1832 1833 /*Intermediaries*/ 1834 IssmDouble Jdet,weight; 1835 IssmDouble thickness,dmudB,n; 1836 IssmDouble dvx[3],dvy[3],dadjbx[3],dadjby[3],dadjshx[3],dadjshy[3]; 1837 IssmDouble *xyz_list= NULL; 1838 1839 /*Fetch number of vertices for this finite element*/ 1840 int numvertices = basalelement->GetNumberOfVertices(); 1841 1842 /*Initialize some vectors*/ 1843 IssmDouble* basis = xNew<IssmDouble>(numvertices); 1844 IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices); 1845 int* vertexpidlist = xNew<int>(numvertices); 1846 IssmDouble zeta; 1847 IssmDouble e1[3],e2[3], phishx[3], phishy[3]; 1848 IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ 1849 IssmDouble adjshx, adjshy; 1850 1851 /*Retrieve all inputs we will be needing: */ 1852 basalelement->GetVerticesCoordinates(&xyz_list); 1853 basalelement->GradientIndexing(&vertexpidlist[0],control_index); 1854 Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); 1855 Input* n_input = basalelement->GetInput(MaterialsRheologyNEnum); _assert_(n_input); 1856 Input* vxbase_input = basalelement->GetInput(VxBaseEnum); _assert_(vxbase_input); 1857 Input* vybase_input = basalelement->GetInput(VyBaseEnum); _assert_(vybase_input); 1858 Input* vxshear_input = basalelement->GetInput(VxShearEnum); _assert_(vxshear_input); 1859 Input* vyshear_input = basalelement->GetInput(VyShearEnum); _assert_(vyshear_input); 1860 1861 Input* adjointbx_input = basalelement->GetInput(AdjointxBaseEnum); _assert_(adjointbx_input); 1862 Input* adjointby_input = basalelement->GetInput(AdjointyBaseEnum); _assert_(adjointby_input); 1863 Input* adjointshx_input = basalelement->GetInput(AdjointxShearEnum); _assert_(adjointshx_input); 1864 Input* adjointshy_input = basalelement->GetInput(AdjointyShearEnum); _assert_(adjointshy_input); 1865 1866 /* Start looping on the number of gaussian points: */ 1867 Gauss* gauss=basalelement->NewGauss(5); 1868 GaussSeg* gauss_seg=new GaussSeg(2); 1869 1870 while(gauss->next()){ 1871 1872 thickness_input->GetInputValue(&thickness,gauss); 1873 n_input->GetInputValue(&n,gauss); 1874 adjointbx_input->GetInputDerivativeValue(&dadjbx[0],xyz_list,gauss); 1875 adjointby_input->GetInputDerivativeValue(&dadjby[0],xyz_list,gauss); 1876 adjointshx_input->GetInputDerivativeValue(&dadjshx[0],xyz_list,gauss); 1877 adjointshy_input->GetInputDerivativeValue(&dadjshy[0],xyz_list,gauss); 1878 adjointshx_input->GetInputValue(&adjshx, gauss); 1879 adjointshy_input->GetInputValue(&adjshy, gauss); 1880 1881 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); 1882 basalelement->NodalFunctionsP1(basis,gauss); 1883 1884 /* Get the integration in the vertical direction */ 1885 gauss_seg->Reset(); 1886 while(gauss_seg->next()){ 1887 zeta=0.5*(gauss_seg->coord1+1); 1888 1889 basalelement->StrainRateMLHO(&epsilon[0],xyz_list,gauss, 1890 vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); 1891 1892 basalelement->dViscositydBMLHO(&dmudB,dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); 1893 1894 e1[0] = 2.0*epsilon[0]+epsilon[1]; 1895 e1[1] = epsilon[2]; 1896 e1[2] = epsilon[3]; 1897 1898 e2[0] = epsilon[2]; 1899 e2[1] = epsilon[0]+2.0*epsilon[1]; 1900 e2[2] = epsilon[4]; 1901 1902 phishx[0] = dadjbx[0] + (1-pow(zeta, n+1))*dadjshx[0]; 1903 phishx[1] = dadjbx[1] + (1-pow(zeta, n+1))*dadjshx[1]; 1904 phishx[2] = (n+1)/thickness*pow(zeta, n)*adjshx; 1905 phishy[0] = dadjby[0] + (1-pow(zeta, n+1))*dadjshy[0]; 1906 phishy[1] = dadjby[1] + (1-pow(zeta, n+1))*dadjshy[1]; 1907 phishy[2] = (n+1)/thickness*pow(zeta, n)*adjshy; 1908 1909 /*Build gradient vector (actually -dJ/dB): */ 1910 for(int i=0;i<numvertices;i++){ 1911 for (int j=0;j<3;j++){ 1912 ge[i]+=(-dmudB)*2*(e1[j]*phishx[j]+e2[j]*phishy[j])*Jdet*gauss->weight*basis[i]*thickness*0.5*gauss_seg->weight; 1913 _assert_(!xIsNan<IssmDouble>(ge[i])); 1914 } 1915 } 1916 } 1917 } 1918 if(control_interp==P1Enum){ 1919 gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); 1920 } 1921 else{ 1922 _error_("not supported"); 1923 } 1924 1925 /*Clean up and return*/ 1926 xDelete<IssmDouble>(xyz_list); 1927 xDelete<IssmDouble>(basis); 1928 xDelete<IssmDouble>(ge); 1929 xDelete<int>(vertexpidlist); 1930 delete gauss; 1931 delete gauss_seg; 1932 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; 1268 1933 }/*}}}*/ 1269 1934 void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/ … … 1580 2245 delete gauss; 1581 2246 }/*}}}*/ 2247 void AdjointHorizAnalysis::GradientJBMLHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/ 2248 _error_("not implemented yet..."); 2249 }/*}}}*/ 1582 2250 void AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/ 1583 2251 … … 1977 2645 delete friction; 1978 2646 }/*}}}*/ 2647 void AdjointHorizAnalysis::GradientJDragMLHO(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/ 2648 2649 /*return if floating (gradient is 0)*/ 2650 if(element->IsAllFloating()) return; 2651 2652 /*Intermediaries*/ 2653 int domaintype,dim; 2654 Element* basalelement; 2655 2656 /*Get basal element*/ 2657 element->FindParam(&domaintype,DomainTypeEnum); 2658 switch(domaintype){ 2659 case Domain2DhorizontalEnum: 2660 basalelement = element; 2661 dim = 2; 2662 break; 2663 case Domain2DverticalEnum: 2664 if(!element->IsOnBase()) return; 2665 basalelement = element->SpawnBasalElement(); 2666 dim = 1; 2667 break; 2668 case Domain3DEnum: 2669 if(!element->IsOnBase()) return; 2670 basalelement = element->SpawnBasalElement(); 2671 dim = 2; 2672 break; 2673 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2674 } 2675 2676 /*Intermediaries*/ 2677 IssmDouble Jdet,weight; 2678 IssmDouble drag,dalpha2dk; 2679 IssmDouble vx,vy,lambda,mu; 2680 IssmDouble *xyz_list= NULL; 2681 2682 /*Fetch number of vertices for this finite element*/ 2683 int numvertices = basalelement->GetNumberOfVertices(); 2684 2685 /*Initialize some vectors*/ 2686 IssmDouble* basis = xNew<IssmDouble>(numvertices); 2687 IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices); 2688 int* vertexpidlist = xNew<int>(numvertices); 2689 2690 /*Build friction element, needed later: */ 2691 Friction* friction=new Friction(basalelement,dim); 2692 2693 /*Retrieve all inputs we will be needing: */ 2694 basalelement->GetVerticesCoordinates(&xyz_list); 2695 basalelement->GradientIndexing(&vertexpidlist[0],control_index); 2696 Input* vx_input = basalelement->GetInput(VxBaseEnum); _assert_(vx_input); 2697 Input* vy_input = basalelement->GetInput(VyBaseEnum); _assert_(vy_input); 2698 Input* adjointx_input = basalelement->GetInput(AdjointxBaseEnum); _assert_(adjointx_input); 2699 Input* adjointy_input = basalelement->GetInput(AdjointyBaseEnum); _assert_(adjointy_input); 2700 2701 /* get the friction law: 1- Budd, 11-Schoof*/ 2702 int frictionlaw;element->FindParam(&frictionlaw, FrictionLawEnum); 2703 Input* dragcoeff_input = NULL; 2704 switch(frictionlaw) { 2705 case 1: 2706 dragcoeff_input = basalelement->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); 2707 break; 2708 case 2: 2709 case 11: 2710 dragcoeff_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoeff_input); 2711 break; 2712 default: 2713 _error_("Friction law "<< frictionlaw <<" not supported in the inversion."); 2714 } 2715 2716 /* Start looping on the number of gaussian points: */ 2717 Gauss* gauss=basalelement->NewGauss(4); 2718 while(gauss->next()){ 2719 2720 adjointx_input->GetInputValue(&lambda, gauss); 2721 adjointy_input->GetInputValue(&mu, gauss); 2722 vx_input->GetInputValue(&vx,gauss); 2723 vy_input->GetInputValue(&vy,gauss); 2724 dragcoeff_input->GetInputValue(&drag, gauss); 2725 2726 friction->GetAlphaComplement(&dalpha2dk,gauss); 2727 2728 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); 2729 basalelement->NodalFunctionsP1(basis,gauss); 2730 2731 /*Build gradient vector (actually -dJ/dD): */ 2732 if(control_interp==P1Enum){ 2733 for(int i=0;i<numvertices;i++){ 2734 ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i]; 2735 _assert_(!xIsNan<IssmDouble>(ge[i])); 2736 } 2737 } 2738 else if(control_interp==P0Enum){ 2739 ge[0]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight; 2740 _assert_(!xIsNan<IssmDouble>(ge[0])); 2741 } 2742 else{ 2743 _error_("not supported"); 2744 } 2745 } 2746 if(control_interp==P1Enum){ 2747 gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); 2748 } 2749 else if(control_interp==P0Enum){ 2750 gradient->SetValue(vertexpidlist[0],ge[0],ADD_VAL); 2751 } 2752 else{ 2753 _error_("not supported"); 2754 } 2755 2756 /*Clean up and return*/ 2757 xDelete<IssmDouble>(xyz_list); 2758 xDelete<IssmDouble>(basis); 2759 xDelete<IssmDouble>(ge); 2760 xDelete<int>(vertexpidlist); 2761 delete gauss; 2762 delete friction; 2763 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; 2764 }/*}}}*/ 1979 2765 void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_interp,int control_index){/*{{{*/ 1980 2766 … … 2480 3266 InputUpdateFromSolutionFS(solution,element); 2481 3267 } 3268 else if (approximation==MLHOApproximationEnum) { 3269 InputUpdateFromSolutionMLHO(solution, element); 3270 } 2482 3271 else{ 2483 3272 InputUpdateFromSolutionHoriz(solution,element); … … 2625 3414 xDelete<int>(doflist); 2626 3415 }/*}}}*/ 3416 void AdjointHorizAnalysis::InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element){/*{{{*/ 3417 int i; 3418 int* doflist=NULL; 3419 3420 int domaintype; 3421 element->FindParam(&domaintype,DomainTypeEnum); 3422 if (domaintype!=Domain2DhorizontalEnum) _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 3423 3424 /*Fetch number of nodes and dof for this finite element*/ 3425 int numnodes = element->GetNumberOfNodes(); 3426 int numdof = numnodes * 4; 3427 3428 /*Fetch dof list and allocate solution vectors*/ 3429 element->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); 3430 IssmDouble* values = xNew<IssmDouble>(numdof); 3431 IssmDouble* lambdax = xNew<IssmDouble>(numnodes); 3432 IssmDouble* lambday = xNew<IssmDouble>(numnodes); 3433 IssmDouble* lambdabx = xNew<IssmDouble>(numnodes); 3434 IssmDouble* lambdaby = xNew<IssmDouble>(numnodes); 3435 IssmDouble* lambdashx = xNew<IssmDouble>(numnodes); 3436 IssmDouble* lambdashy = xNew<IssmDouble>(numnodes); 3437 IssmDouble* n = xNew<IssmDouble>(numnodes); 3438 3439 /*Use the dof list to index into the solution vector: */ 3440 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]]; 3441 3442 /*Transform solution in Cartesian Space*/ 3443 if(domaintype!=Domain2DverticalEnum) element->TransformSolutionCoord(&values[0],XYEnum); 3444 3445 element->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum,0.); 3446 3447 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ 3448 for(i=0;i<numnodes;i++){ 3449 lambdabx[i] =values[i*4+0]; 3450 lambdashx[i]=values[i*4+1]; 3451 /*Check solution*/ 3452 if(xIsNan<IssmDouble>(lambdabx[i])) _error_("NaN found in solution vector"); 3453 if(xIsInf<IssmDouble>(lambdabx[i])) _error_("Inf found in solution vector"); 3454 if(xIsNan<IssmDouble>(lambdashx[i])) _error_("NaN found in solution vector"); 3455 if(xIsInf<IssmDouble>(lambdashx[i])) _error_("Inf found in solution vector"); 3456 /* adjoint for the surface velocity */ 3457 lambdax[i] = lambdabx[i] + lambdashx[i]*(n[i]+1)/(n[i]+2); 3458 3459 lambdaby[i] =values[i*4+2]; 3460 lambdashy[i]=values[i*4+3]; 3461 /*Check solution*/ 3462 if(xIsNan<IssmDouble>(lambdaby[i])) _error_("NaN found in solution vector"); 3463 if(xIsInf<IssmDouble>(lambdaby[i])) _error_("Inf found in solution vector"); 3464 if(xIsNan<IssmDouble>(lambdashy[i])) _error_("NaN found in solution vector"); 3465 if(xIsInf<IssmDouble>(lambdashy[i])) _error_("Inf found in solution vector"); 3466 /* adjoint for the surface velocity */ 3467 lambday[i] = lambdaby[i] + lambdashy[i]*(n[i]+1)/(n[i]+2); 3468 } 3469 3470 /*Add vx and vy as inputs to the tria element: */ 3471 element->AddInput(AdjointxBaseEnum,lambdabx,element->GetElementType()); 3472 element->AddInput(AdjointyBaseEnum,lambdaby,element->GetElementType()); 3473 element->AddInput(AdjointxShearEnum,lambdashx,element->GetElementType()); 3474 element->AddInput(AdjointyShearEnum,lambdashy,element->GetElementType()); 3475 3476 element->AddInput(AdjointxEnum,lambdax,element->GetElementType()); 3477 element->AddInput(AdjointyEnum,lambday,element->GetElementType()); 3478 3479 /*Free ressources:*/ 3480 xDelete<IssmDouble>(values); 3481 xDelete<IssmDouble>(lambdax); 3482 xDelete<IssmDouble>(lambday); 3483 xDelete<IssmDouble>(lambdabx); 3484 xDelete<IssmDouble>(lambdaby); 3485 xDelete<IssmDouble>(lambdashx); 3486 xDelete<IssmDouble>(lambdashy); 3487 xDelete<int>(doflist); 3488 }/*}}}*/ 2627 3489 void AdjointHorizAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ 2628 3490 /*Default, do nothing*/
Note:
See TracChangeset
for help on using the changeset viewer.
