![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 18834
- Timestamp:
- 11/24/14 10:37:06 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
r18771 r18834 337 337 338 338 /*Intermediaries */ 339 int num_responses,i,d im;339 int num_responses,i,domaintype; 340 340 IssmDouble Jdet,obs_velocity_mag,velocity_mag; 341 341 IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; … … 344 344 IssmDouble *xyz_list_top = NULL; 345 345 346 /* Get problem dimension*/347 element->FindParam(&d im,DomainDimensionEnum);346 /* Get domaintype*/ 347 element->FindParam(&domaintype,DomainTypeEnum); 348 348 349 349 /*Fetch number of nodes and dof for this finite element*/ … … 353 353 /*Prepare coordinate system list*/ 354 354 int* cs_list = xNew<int>(vnumnodes+pnumnodes); 355 if(d im==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;355 if(domaintype==Domain2DverticalEnum) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum; 356 356 else for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum; 357 357 for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum; … … 367 367 Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 368 368 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); 369 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);370 369 Input* vxobs_input = element->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 371 Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 370 Input* vy_input=NULL; 371 Input* vyobs_input=NULL; 372 if (domaintype!=Domain2DverticalEnum) { 373 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); 374 Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 375 } 372 376 IssmDouble epsvel = 2.220446049250313e-16; 373 377 IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ … … 389 393 390 394 vx_input->GetInputValue(&vx,gauss); 391 vy_input->GetInputValue(&vy,gauss);392 395 vxobs_input->GetInputValue(&vxobs,gauss); 393 vyobs_input->GetInputValue(&vyobs,gauss); 396 if (domaintype!=Domain2DverticalEnum) { 397 vy_input->GetInputValue(&vy,gauss); 398 vyobs_input->GetInputValue(&vyobs,gauss); 399 } 394 400 395 401 /*Loop over all requested responses*/ … … 409 415 */ 410 416 for(i=0;i<vnumnodes;i++){ 411 dux=vxobs-vx;412 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];413 if(dim==3){417 if(domaintype!=Domain2DverticalEnum){ 418 dux=vxobs-vx; 419 pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 414 420 duy=vyobs-vy; 415 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 421 pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 422 } 423 else{ 424 dux=vxobs-vx; 425 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 416 426 } 417 427 } … … 430 440 */ 431 441 for(i=0;i<vnumnodes;i++){ 432 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 433 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 434 dux=scalex*(vxobs-vx); 435 duy=scaley*(vyobs-vy); 436 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 437 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 442 if(domaintype!=Domain2DverticalEnum){ 443 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 444 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 445 dux=scalex*(vxobs-vx); 446 duy=scaley*(vyobs-vy); 447 pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 448 pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 449 } 450 else{ 451 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 452 dux=scalex*(vxobs-vx); 453 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 454 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 455 } 438 456 } 439 457 break; … … 451 469 */ 452 470 for(i=0;i<vnumnodes;i++){ 453 if(d im==3){471 if(domaintype!=Domain2DverticalEnum){ 454 472 velocity_mag =sqrt(vx*vx+vy*vy)+epsvel; 455 473 obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel; … … 457 475 dux=scale*vx; 458 476 duy=scale*vy; 459 pe->values[i* dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];460 pe->values[i* dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i];477 pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 478 pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 461 479 } 462 480 else{ … … 465 483 scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag); 466 484 dux=scale*vx; 467 pe->values[i* dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i];485 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 468 486 } 469 487 } … … 480 498 */ 481 499 for(i=0;i<vnumnodes;i++){ 482 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 483 dux=scale*(vxobs-vx); 484 duy=scale*(vyobs-vy); 485 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 486 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 500 if (domaintype!=Domain2DverticalEnum){ 501 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 502 dux=scale*(vxobs-vx); 503 duy=scale*(vyobs-vy); 504 pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 505 pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 506 } 507 else{ 508 scale=1./(S*2*fabs(vx-vxobs)+epsvel); 509 dux=scale*(vxobs-vx); 510 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 511 } 487 512 } 488 513 break; … … 498 523 */ 499 524 for(i=0;i<vnumnodes;i++){ 500 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 501 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 502 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 503 pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 525 if(domaintype!=Domain2DverticalEnum){ 526 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 527 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 528 pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 529 pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 530 } 531 else{ 532 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 533 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 534 } 504 535 } 505 536 break; … … 542 573 543 574 /*Intermediaries */ 544 int num_responses,i ;575 int num_responses,i,domaintype; 545 576 IssmDouble Jdet,obs_velocity_mag,velocity_mag; 546 577 IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; … … 548 579 int *responses = NULL; 549 580 IssmDouble *xyz_list_top = NULL; 581 582 /* Get domaintype*/ 583 element->FindParam(&domaintype,DomainTypeEnum); 550 584 551 585 /*Fetch number of nodes and dof for this finite element*/ … … 562 596 Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 563 597 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); 564 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);565 598 Input* vxobs_input = element->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 566 Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 599 Input* vy_input=NULL; 600 Input* vyobs_input=NULL; 601 if(domaintype!=Domain2DverticalEnum){ 602 vy_input = element->GetInput(VyEnum); _assert_(vy_input); 603 vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 604 } 567 605 IssmDouble epsvel = 2.220446049250313e-16; 568 606 IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ … … 584 622 585 623 vx_input->GetInputValue(&vx,gauss); 586 vy_input->GetInputValue(&vy,gauss);587 624 vxobs_input->GetInputValue(&vxobs,gauss); 588 vyobs_input->GetInputValue(&vyobs,gauss); 589 625 if(domaintype!=Domain2DverticalEnum){ 626 vy_input->GetInputValue(&vy,gauss); 627 vyobs_input->GetInputValue(&vyobs,gauss); 628 } 590 629 /*Loop over all requested responses*/ 591 630 for(int resp=0;resp<num_responses;resp++){ … … 604 643 */ 605 644 for(i=0;i<numnodes;i++){ 606 dux=vxobs-vx; 607 duy=vyobs-vy; 608 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 609 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 645 if(domaintype!=Domain2DverticalEnum){ 646 dux=vxobs-vx; 647 duy=vyobs-vy; 648 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 649 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 650 } 651 else{ 652 dux=vxobs-vx; 653 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 654 } 610 655 } 611 656 break; … … 623 668 */ 624 669 for(i=0;i<numnodes;i++){ 625 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 626 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 627 dux=scalex*(vxobs-vx); 628 duy=scaley*(vyobs-vy); 629 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 630 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 670 if(domaintype!=Domain2DverticalEnum){ 671 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 672 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 673 dux=scalex*(vxobs-vx); 674 duy=scaley*(vyobs-vy); 675 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 676 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 677 } 678 else{ 679 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 680 dux=scalex*(vxobs-vx); 681 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 682 } 631 683 } 632 684 break; … … 644 696 */ 645 697 for(i=0;i<numnodes;i++){ 646 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 647 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 648 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 649 dux=scale*vx; 650 duy=scale*vy; 651 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 652 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 698 if(domaintype!=Domain2DverticalEnum){ 699 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 700 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 701 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 702 dux=scale*vx; 703 duy=scale*vy; 704 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 705 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 706 } 707 else{ 708 velocity_mag =fabs(vx)+epsvel; 709 obs_velocity_mag=fabs(vxobs)+epsvel; 710 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 711 dux=scale*vx; 712 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 713 } 653 714 } 654 715 break; … … 664 725 */ 665 726 for(i=0;i<numnodes;i++){ 666 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 667 dux=scale*(vxobs-vx); 668 duy=scale*(vyobs-vy); 669 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 670 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 727 if(domaintype!=Domain2DverticalEnum){ 728 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 729 dux=scale*(vxobs-vx); 730 duy=scale*(vyobs-vy); 731 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 732 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 733 } 734 else{ 735 scale=1./(S*2*fabs(vx-vxobs)+epsvel); 736 dux=scale*(vxobs-vx); 737 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 738 } 671 739 } 672 740 break; … … 682 750 */ 683 751 for(i=0;i<numnodes;i++){ 684 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 685 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 686 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 687 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 752 if(domaintype!=Domain2DverticalEnum){ 753 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 754 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 755 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 756 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 757 } 758 else{ 759 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 760 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 761 } 688 762 } 689 763 break; … … 701 775 break; 702 776 case RheologyBbarAbsGradientEnum: 777 /*Nothing in P vector*/ 778 break; 779 case RheologyBAbsGradientEnum: 703 780 /*Nothing in P vector*/ 704 781 break; … … 710 787 711 788 /*Transform coordinate system*/ 712 element->TransformLoadVectorCoord(pe,XYEnum);789 if(domaintype!=Domain2DverticalEnum) element->TransformLoadVectorCoord(pe,XYEnum); 713 790 714 791 /*Clean up and return*/ … … 731 808 case Domain2DhorizontalEnum: 732 809 basalelement = element; 810 break; 811 case Domain2DverticalEnum: 812 if(!element->IsOnBase()) return NULL; 813 basalelement = element->SpawnBasalElement(); 733 814 break; 734 815 case Domain3DEnum: … … 760 841 Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 761 842 Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); 762 Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input);763 843 Input* vxobs_input = basalelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); 764 Input* vyobs_input = basalelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 844 Input* vy_input=NULL; 845 Input* vyobs_input=NULL; 846 if(domaintype!=Domain2DverticalEnum){ 847 vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); 848 vyobs_input = basalelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); 849 } 765 850 IssmDouble epsvel = 2.220446049250313e-16; 766 851 IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ … … 782 867 783 868 vx_input->GetInputValue(&vx,gauss); 784 vy_input->GetInputValue(&vy,gauss);785 869 vxobs_input->GetInputValue(&vxobs,gauss); 786 vyobs_input->GetInputValue(&vyobs,gauss); 787 870 if(domaintype!=Domain2DverticalEnum){ 871 vy_input->GetInputValue(&vy,gauss); 872 vyobs_input->GetInputValue(&vyobs,gauss); 873 } 788 874 /*Loop over all requested responses*/ 789 875 for(int resp=0;resp<num_responses;resp++){ … … 802 888 */ 803 889 for(i=0;i<numnodes;i++){ 804 dux=vxobs-vx; 805 duy=vyobs-vy; 806 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 807 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 890 if(domaintype!=Domain2DverticalEnum){ 891 dux=vxobs-vx; 892 duy=vyobs-vy; 893 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 894 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 895 } 896 else { 897 dux=vxobs-vx; 898 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 899 } 808 900 } 809 901 break; … … 821 913 */ 822 914 for(i=0;i<numnodes;i++){ 823 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 824 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 825 dux=scalex*(vxobs-vx); 826 duy=scaley*(vyobs-vy); 827 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 828 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 915 if(domaintype!=Domain2DverticalEnum){ 916 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 917 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; 918 dux=scalex*(vxobs-vx); 919 duy=scaley*(vyobs-vy); 920 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 921 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 922 } 923 else{ 924 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; 925 dux=scalex*(vxobs-vx); 926 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 927 } 829 928 } 830 929 break; … … 842 941 */ 843 942 for(i=0;i<numnodes;i++){ 844 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 845 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 846 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 847 dux=scale*vx; 848 duy=scale*vy; 849 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 850 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 943 if(domaintype!=Domain2DverticalEnum){ 944 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; 945 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; 946 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 947 dux=scale*vx; 948 duy=scale*vy; 949 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 950 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 951 } 952 else{ 953 velocity_mag =fabs(vx)+epsvel; 954 obs_velocity_mag=fabs(vxobs)+epsvel; 955 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); 956 dux=scale*vx; 957 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 958 } 851 959 } 852 960 break; … … 862 970 */ 863 971 for(i=0;i<numnodes;i++){ 864 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 865 dux=scale*(vxobs-vx); 866 duy=scale*(vyobs-vy); 867 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 868 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 972 if(domaintype!=Domain2DverticalEnum){ 973 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); 974 dux=scale*(vxobs-vx); 975 duy=scale*(vyobs-vy); 976 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 977 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 978 } 979 else{ 980 scale=1./(S*2*fabs(vx-vxobs)+epsvel); 981 dux=scale*(vxobs-vx); 982 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 983 } 869 984 } 870 985 break; … … 880 995 */ 881 996 for(i=0;i<numnodes;i++){ 882 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 883 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 884 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 885 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 997 if(domaintype!=Domain2DverticalEnum){ 998 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 999 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); 1000 pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 1001 pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 1002 } 1003 else{ 1004 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); 1005 pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; 1006 } 886 1007 } 887 1008 break; … … 899 1020 break; 900 1021 case RheologyBbarAbsGradientEnum: 1022 /*Nothing in P vector*/ 1023 break; 1024 case RheologyBAbsGradientEnum: 901 1025 /*Nothing in P vector*/ 902 1026 break; … … 908 1032 909 1033 /*Transform coordinate system*/ 910 basalelement->TransformLoadVectorCoord(pe,XYEnum);1034 if(domaintype!=Domain2DverticalEnum) basalelement->TransformLoadVectorCoord(pe,XYEnum); 911 1035 912 1036 /*Clean up and return*/ … … 951 1075 if(control_type!=MaterialsRheologyBbarEnum && 952 1076 control_type!=FrictionCoefficientEnum && 953 control_type!=DamageDbarEnum){ 1077 control_type!=DamageDbarEnum && 1078 control_type!=MaterialsRheologyBEnum){ 954 1079 _error_("Control "<<EnumToStringx(control_type)<<" not supported"); 955 1080 } … … 963 1088 case SurfaceAverageVelMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break; 964 1089 case DragCoefficientAbsGradientEnum: GradientJDragGradient(element,gradient,control_index); break; 965 case RheologyBbarAbsGradientEnum: GradientJBGradient(element,gradient,control_index); break; 1090 case RheologyBbarAbsGradientEnum: GradientJBbarGradient(element,gradient,control_index); break; 1091 case RheologyBAbsGradientEnum: GradientJBGradient(element,gradient,control_index); break; 966 1092 default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); 967 1093 } … … 989 1115 } 990 1116 break; 1117 case MaterialsRheologyBEnum: 1118 switch(approximation){ 1119 case SSAApproximationEnum: GradientJBSSA(element,gradient,control_index); break; 1120 case HOApproximationEnum: GradientJBHO( element,gradient,control_index); break; 1121 case FSApproximationEnum: GradientJBFS( element,gradient,control_index); break; 1122 case NoneApproximationEnum: /*Gradient is 0*/ break; 1123 default: _error_("approximation " << EnumToStringx(approximation) << " not supported yet"); 1124 } 1125 break; 991 1126 case DamageDbarEnum: 992 1127 switch(approximation){ … … 1085 1220 1086 1221 }/*}}}*/ 1087 void AdjointHorizAnalysis::GradientJB Gradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/1222 void AdjointHorizAnalysis::GradientJBbarGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ 1088 1223 1089 1224 /*Intermediaries*/ 1090 int domaintype ,dim;1225 int domaintype; 1091 1226 Element* basalelement; 1092 1227 … … 1096 1231 case Domain2DhorizontalEnum: 1097 1232 basalelement = element; 1098 dim = 2;1099 1233 break; 1100 1234 case Domain2DverticalEnum: 1101 1235 if(!element->IsOnBase()) return; 1102 1236 basalelement = element->SpawnBasalElement(); 1103 dim = 1;1104 1237 break; 1105 1238 case Domain3DEnum: 1106 1239 if(!element->IsOnBase()) return; 1107 1240 basalelement = element->SpawnBasalElement(); 1108 dim = 2;1109 1241 break; 1110 1242 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); … … 1144 1276 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ 1145 1277 for(int i=0;i<numvertices;i++){ 1146 if(d im==2){1278 if(domaintype!=Domain2DverticalEnum){ 1147 1279 ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]); 1148 1280 } … … 1162 1294 delete gauss; 1163 1295 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; 1296 }/*}}}*/ 1297 void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ 1298 1299 /*Intermediaries*/ 1300 int domaintype; 1301 1302 /*Get basal element*/ 1303 element->FindParam(&domaintype,DomainTypeEnum); 1304 switch(domaintype){ 1305 case Domain2DhorizontalEnum: 1306 break; 1307 case Domain2DverticalEnum: 1308 break; 1309 case Domain3DEnum: 1310 break; 1311 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 1312 } 1313 1314 /*Intermediaries*/ 1315 IssmDouble Jdet,weight; 1316 IssmDouble dk[3]; 1317 IssmDouble *xyz_list= NULL; 1318 1319 /*Fetch number of vertices for this finite element*/ 1320 int numvertices = element->GetNumberOfVertices(); 1321 1322 /*Initialize some vectors*/ 1323 IssmDouble* dbasis = xNew<IssmDouble>(3*numvertices); 1324 IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices); 1325 int* vertexpidlist = xNew<int>(numvertices); 1326 1327 /*Retrieve all inputs we will be needing: */ 1328 element->GetVerticesCoordinates(&xyz_list); 1329 element->GradientIndexing(&vertexpidlist[0],control_index); 1330 Input* rheology_input = element->GetInput(MaterialsRheologyBEnum); _assert_(rheology_input); 1331 Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); 1332 /* Start looping on the number of gaussian points: */ 1333 Gauss* gauss=element->NewGauss(2); 1334 for(int ig=gauss->begin();ig<gauss->end();ig++){ 1335 gauss->GaussPoint(ig); 1336 element->JacobianDeterminant(&Jdet,xyz_list,gauss); 1337 element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss); 1338 weights_input->GetInputValue(&weight,gauss,RheologyBAbsGradientEnum); 1339 1340 /*Build alpha_complement_list: */ 1341 rheology_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss); 1342 1343 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ 1344 for(int i=0;i<numvertices;i++){ 1345 if(domaintype!=Domain2DverticalEnum){ 1346 ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]); 1347 } 1348 else{ 1349 ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]+dbasis[2*numvertices+i]*dk[2]); 1350 } 1351 _assert_(!xIsNan<IssmDouble>(ge[i])); 1352 } 1353 } 1354 gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); 1355 1356 /*Clean up and return*/ 1357 xDelete<IssmDouble>(xyz_list); 1358 xDelete<IssmDouble>(dbasis); 1359 xDelete<IssmDouble>(ge); 1360 xDelete<int>(vertexpidlist); 1361 delete gauss; 1164 1362 }/*}}}*/ 1165 1363 void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ … … 1263 1461 1264 1462 /*Intermediaries*/ 1265 int dim=3;1266 1463 IssmDouble Jdet,weight; 1267 1464 IssmDouble drag,dalpha2dk; … … 1269 1466 IssmDouble *xyz_list_base= NULL; 1270 1467 1468 int domaintype,dim; 1469 element->FindParam(&domaintype,DomainTypeEnum); 1271 1470 /*Fetch number of vertices for this finite element*/ 1272 1471 int numvertices = element->GetNumberOfVertices(); … … 1278 1477 1279 1478 /*Build friction element, needed later: */ 1479 if(domaintype!=Domain2DverticalEnum) dim=3; 1480 else dim=2; 1481 Friction* friction=new Friction(element,dim); 1482 1483 /*Retrieve all inputs we will be needing: */ 1484 element->GetVerticesCoordinatesBase(&xyz_list_base); 1485 element->GradientIndexing(&vertexpidlist[0],control_index); 1486 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); 1487 Input* vy_input = NULL; 1488 Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); 1489 Input* adjointy_input = NULL; 1490 Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); 1491 if(domaintype!=Domain2DverticalEnum){ 1492 vy_input = element->GetInput(VyEnum); _assert_(vy_input); 1493 adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); 1494 } 1495 /* Start looping on the number of gaussian points: */ 1496 Gauss* gauss=element->NewGaussBase(4); 1497 for(int ig=gauss->begin();ig<gauss->end();ig++){ 1498 gauss->GaussPoint(ig); 1499 1500 adjointx_input->GetInputValue(&lambda, gauss); 1501 vx_input->GetInputValue(&vx,gauss); 1502 if(domaintype!=Domain2DverticalEnum){ 1503 adjointy_input->GetInputValue(&mu, gauss); 1504 vy_input->GetInputValue(&vy,gauss); 1505 } 1506 dragcoeff_input->GetInputValue(&drag, gauss); 1507 1508 friction->GetAlphaComplement(&dalpha2dk,gauss); 1509 1510 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); 1511 element->NodalFunctionsP1(basis,gauss); 1512 1513 /*Build gradient vector (actually -dJ/dD): */ 1514 for(int i=0;i<numvertices;i++){ 1515 if(domaintype!=Domain2DverticalEnum) ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i]; 1516 else ge[i]+=-2.*drag*dalpha2dk*(lambda*vx)*Jdet*gauss->weight*basis[i]; 1517 _assert_(!xIsNan<IssmDouble>(ge[i])); 1518 } 1519 } 1520 gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); 1521 1522 /*Clean up and return*/ 1523 xDelete<IssmDouble>(xyz_list_base); 1524 xDelete<IssmDouble>(basis); 1525 xDelete<IssmDouble>(ge); 1526 xDelete<int>(vertexpidlist); 1527 delete gauss; 1528 delete friction; 1529 }/*}}}*/ 1530 void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ 1531 1532 /*return if floating or not on bed (gradient is 0)*/ 1533 if(element->IsFloating()) return; 1534 if(!element->IsOnBase()) return; 1535 1536 /*Intermediaries*/ 1537 int domaintype,dim; 1538 IssmDouble Jdet,weight; 1539 IssmDouble drag,dalpha2dk,normal[3]; 1540 IssmDouble vx,vy,vz,lambda,mu,xi; 1541 IssmDouble *xyz_list_base= NULL; 1542 1543 /*Fetch number of vertices for this finite element*/ 1544 int numvertices = element->GetNumberOfVertices(); 1545 1546 /*Initialize some vectors*/ 1547 IssmDouble* basis = xNew<IssmDouble>(numvertices); 1548 IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices); 1549 int* vertexpidlist = xNew<int>(numvertices); 1550 1551 /* get domaintype */ 1552 element->FindParam(&domaintype,DomainTypeEnum); 1553 1554 /*Build friction element, needed later: */ 1555 if(domaintype!=Domain2DverticalEnum) dim=3; 1556 else dim=2; 1280 1557 Friction* friction=new Friction(element,dim); 1281 1558 … … 1287 1564 Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); 1288 1565 Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); 1566 Input* vz_input = NULL; 1567 Input* adjointz_input = NULL; 1568 if(domaintype!=Domain2DverticalEnum){ 1569 Input* vz_input = element->GetInput(VzEnum); _assert_(vy_input); 1570 Input* adjointz_input = element->GetInput(AdjointzEnum); _assert_(adjointz_input); 1571 } 1289 1572 Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); 1290 1573 … … 1298 1581 vx_input->GetInputValue(&vx,gauss); 1299 1582 vy_input->GetInputValue(&vy,gauss); 1583 if(domaintype!=Domain2DverticalEnum){ 1584 adjointz_input->GetInputValue(&xi ,gauss); 1585 vz_input->GetInputValue(&vz,gauss); 1586 } 1300 1587 dragcoeff_input->GetInputValue(&drag, gauss); 1301 1588 1302 1589 friction->GetAlphaComplement(&dalpha2dk,gauss); 1590 element->NormalBase(&normal[0],xyz_list_base); 1303 1591 1304 1592 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); 1305 1593 element->NodalFunctionsP1(basis,gauss); 1306 1594 1307 /*Build gradient vector (actually -dJ/dD): */1308 for(int i=0;i<numvertices;i++){1309 ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];1310 _assert_(!xIsNan<IssmDouble>(ge[i]));1311 }1312 }1313 gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);1314 1315 /*Clean up and return*/1316 xDelete<IssmDouble>(xyz_list_base);1317 xDelete<IssmDouble>(basis);1318 xDelete<IssmDouble>(ge);1319 xDelete<int>(vertexpidlist);1320 delete gauss;1321 delete friction;1322 }/*}}}*/1323 void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/1324 1325 /*return if floating or not on bed (gradient is 0)*/1326 if(element->IsFloating()) return;1327 if(!element->IsOnBase()) return;1328 1329 /*Intermediaries*/1330 int dim=3;1331 IssmDouble Jdet,weight;1332 IssmDouble drag,dalpha2dk,normal[3];1333 IssmDouble vx,vy,vz,lambda,mu,xi;1334 IssmDouble *xyz_list_base= NULL;1335 1336 /*Fetch number of vertices for this finite element*/1337 int numvertices = element->GetNumberOfVertices();1338 1339 /*Initialize some vectors*/1340 IssmDouble* basis = xNew<IssmDouble>(numvertices);1341 IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);1342 int* vertexpidlist = xNew<int>(numvertices);1343 1344 /*Build friction element, needed later: */1345 Friction* friction=new Friction(element,dim);1346 1347 /*Retrieve all inputs we will be needing: */1348 element->GetVerticesCoordinatesBase(&xyz_list_base);1349 element->GradientIndexing(&vertexpidlist[0],control_index);1350 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);1351 Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);1352 Input* vz_input = element->GetInput(VzEnum); _assert_(vy_input);1353 Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input);1354 Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input);1355 Input* adjointz_input = element->GetInput(AdjointzEnum); _assert_(adjointz_input);1356 Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input);1357 1358 /* Start looping on the number of gaussian points: */1359 Gauss* gauss=element->NewGaussBase(4);1360 for(int ig=gauss->begin();ig<gauss->end();ig++){1361 gauss->GaussPoint(ig);1362 1363 adjointx_input->GetInputValue(&lambda, gauss);1364 adjointy_input->GetInputValue(&mu, gauss);1365 adjointz_input->GetInputValue(&xi ,gauss);1366 vx_input->GetInputValue(&vx,gauss);1367 vy_input->GetInputValue(&vy,gauss);1368 vz_input->GetInputValue(&vz,gauss);1369 dragcoeff_input->GetInputValue(&drag, gauss);1370 1371 friction->GetAlphaComplement(&dalpha2dk,gauss);1372 element->NormalBase(&normal[0],xyz_list_base);1373 1374 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);1375 element->NodalFunctionsP1(basis,gauss);1376 1377 1595 /*Build gradient vector (actually -dJ/dk): */ 1378 for(int i=0;i<numvertices;i++){ 1379 ge[i]+=( 1380 -lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2])) 1381 -mu *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2])) 1382 -xi *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2])) 1383 )*Jdet*gauss->weight*basis[i]; 1384 _assert_(!xIsNan<IssmDouble>(ge[i])); 1596 if(domaintype!=Domain2DverticalEnum){ 1597 for(int i=0;i<numvertices;i++){ 1598 ge[i]+=( 1599 -lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2])) 1600 -mu *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2])) 1601 -xi *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2])) 1602 )*Jdet*gauss->weight*basis[i]; 1603 _assert_(!xIsNan<IssmDouble>(ge[i])); 1604 } 1605 } 1606 else{ 1607 for(int i=0;i<numvertices;i++){ 1608 ge[i]+=( 1609 -lambda*2*drag*dalpha2dk*vx 1610 -mu *2*drag*dalpha2dk*vy 1611 )*Jdet*gauss->weight*basis[i]; 1612 _assert_(!xIsNan<IssmDouble>(ge[i])); 1613 } 1385 1614 } 1386 1615 } … … 1493 1722 this->GradientJBbarSSA(element,gradient,control_index); 1494 1723 }/*}}}*/ 1724 void AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ 1725 1726 /*Intermediaries*/ 1727 int domaintype,dim; 1728 Element* basalelement; 1729 1730 /*Get basal element*/ 1731 element->FindParam(&domaintype,DomainTypeEnum); 1732 switch(domaintype){ 1733 case Domain2DhorizontalEnum: 1734 basalelement = element; 1735 dim = 2; 1736 break; 1737 case Domain2DverticalEnum: 1738 if(!element->IsOnBase()) return; 1739 basalelement = element->SpawnBasalElement(); 1740 dim = 1; 1741 break; 1742 case Domain3DEnum: 1743 if(!element->IsOnBase()) return; 1744 basalelement = element->SpawnBasalElement(); 1745 dim = 2; 1746 break; 1747 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 1748 } 1749 1750 /*Intermediaries*/ 1751 IssmDouble Jdet,weight; 1752 IssmDouble thickness,dmudB; 1753 IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; 1754 IssmDouble *xyz_list= NULL; 1755 1756 /*Fetch number of vertices for this finite element*/ 1757 int numvertices = basalelement->GetNumberOfVertices(); 1758 1759 /*Initialize some vectors*/ 1760 IssmDouble* basis = xNew<IssmDouble>(numvertices); 1761 IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices); 1762 int* vertexpidlist = xNew<int>(numvertices); 1763 1764 /*Retrieve all inputs we will be needing: */ 1765 basalelement->GetVerticesCoordinates(&xyz_list); 1766 basalelement->GradientIndexing(&vertexpidlist[0],control_index); 1767 Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); 1768 Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); 1769 Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); 1770 Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input); 1771 Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input); 1772 Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input); 1773 1774 /* Start looping on the number of gaussian points: */ 1775 Gauss* gauss=basalelement->NewGauss(4); 1776 for(int ig=gauss->begin();ig<gauss->end();ig++){ 1777 gauss->GaussPoint(ig); 1778 1779 thickness_input->GetInputValue(&thickness,gauss); 1780 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); 1781 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); 1782 adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss); 1783 adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss); 1784 1785 basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); 1786 1787 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); 1788 basalelement->NodalFunctionsP1(basis,gauss); 1789 1790 /*Build gradient vector (actually -dJ/dB): */ 1791 for(int i=0;i<numvertices;i++){ 1792 ge[i]+=-dmudB*thickness*( 1793 (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1] 1794 )*Jdet*gauss->weight*basis[i]; 1795 _assert_(!xIsNan<IssmDouble>(ge[i])); 1796 } 1797 } 1798 gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); 1799 1800 /*Clean up and return*/ 1801 xDelete<IssmDouble>(xyz_list); 1802 xDelete<IssmDouble>(basis); 1803 xDelete<IssmDouble>(ge); 1804 xDelete<int>(vertexpidlist); 1805 delete gauss; 1806 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; 1807 }/*}}}*/ 1808 void AdjointHorizAnalysis::GradientJBHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ 1809 /*Intermediaries*/ 1810 int domaintype,dim; 1811 1812 /*Get domaintype*/ 1813 element->FindParam(&domaintype,DomainTypeEnum); 1814 1815 /*Intermediaries*/ 1816 IssmDouble Jdet,weight; 1817 IssmDouble thickness,dmudB; 1818 IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; 1819 IssmDouble *xyz_list= NULL; 1820 1821 /*Fetch number of vertices for this finite element*/ 1822 int numvertices = element->GetNumberOfVertices(); 1823 1824 /*Initialize some vectors*/ 1825 IssmDouble* basis = xNew<IssmDouble>(numvertices); 1826 IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices); 1827 int* vertexpidlist = xNew<int>(numvertices); 1828 1829 /*Retrieve all inputs we will be needing: */ 1830 element->GetVerticesCoordinates(&xyz_list); 1831 element->GradientIndexing(&vertexpidlist[0],control_index); 1832 Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); 1833 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); 1834 Input* vy_input = NULL; 1835 Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); 1836 Input* adjointy_input = NULL; 1837 Input* rheologyb_input = element->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input); 1838 if(domaintype!=Domain2DverticalEnum){ 1839 vy_input = element->GetInput(VyEnum); _assert_(vy_input); 1840 adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); 1841 } 1842 /* Start looping on the number of gaussian points: */ 1843 Gauss* gauss=element->NewGauss(4); 1844 for(int ig=gauss->begin();ig<gauss->end();ig++){ 1845 gauss->GaussPoint(ig); 1846 1847 thickness_input->GetInputValue(&thickness,gauss); 1848 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); 1849 adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss); 1850 dim=2; 1851 if(domaintype!=Domain2DverticalEnum){ 1852 adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list, gauss); 1853 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); 1854 dim=3; 1855 } 1856 1857 element->dViscositydBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); 1858 1859 element->JacobianDeterminant(&Jdet,xyz_list,gauss); 1860 element->NodalFunctionsP1(basis,gauss); 1861 1862 /*Build gradient vector (actually -dJ/dB): */ 1863 for(int i=0;i<numvertices;i++){ 1864 if(domaintype!=Domain2DverticalEnum){ 1865 ge[i]+=-dmudB*thickness*( 1866 (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1] 1867 )*Jdet*gauss->weight*basis[i]; 1868 } 1869 else{ 1870 ge[i]+=-dmudB*thickness*4*dvx[0]*dadjx[0]*Jdet*gauss->weight*basis[i]; 1871 } 1872 _assert_(!xIsNan<IssmDouble>(ge[i])); 1873 } 1874 } 1875 gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); 1876 1877 /*Clean up and return*/ 1878 xDelete<IssmDouble>(xyz_list); 1879 xDelete<IssmDouble>(basis); 1880 xDelete<IssmDouble>(ge); 1881 xDelete<int>(vertexpidlist); 1882 delete gauss; 1883 }/*}}}*/ 1884 void AdjointHorizAnalysis::GradientJBFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ 1885 /*WARNING: We use HO as an estimate for now*/ 1886 this->GradientJBHO(element,gradient,control_index); 1887 }/*}}}*/ 1495 1888 void AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/ 1496 1889 … … 1590 1983 int* doflist=NULL; 1591 1984 1985 int domaintype; 1986 element->FindParam(&domaintype,DomainTypeEnum); 1987 1592 1988 /*Fetch number of nodes and dof for this finite element*/ 1593 1989 int numnodes = element->GetNumberOfNodes(); 1594 int numdof = numnodes*2; 1595 1990 int numdof; 1991 if(domaintype!=Domain2DverticalEnum) numdof = numnodes*2; 1992 else numdof = numnodes*1; 1596 1993 /*Fetch dof list and allocate solution vectors*/ 1597 1994 element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); … … 1604 2001 1605 2002 /*Transform solution in Cartesian Space*/ 1606 element->TransformSolutionCoord(&values[0],XYEnum);2003 if(domaintype!=Domain2DverticalEnum) element->TransformSolutionCoord(&values[0],XYEnum); 1607 2004 1608 2005 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ 1609 2006 for(i=0;i<numnodes;i++){ 1610 lambdax[i]=values[i*NDOF2+0]; 1611 lambday[i]=values[i*NDOF2+1]; 1612 2007 if(domaintype!=Domain2DverticalEnum){ 2008 lambdax[i]=values[i*NDOF2+0]; 2009 lambday[i]=values[i*NDOF2+1]; 2010 } 2011 else {lambdax[i]=values[i];lambday[i]=0;} 1613 2012 /*Check solution*/ 1614 2013 if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector"); 1615 if( xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");2014 if(domaintype!=Domain2DverticalEnum && xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector"); 1616 2015 } 1617 2016 … … 1627 2026 }/*}}}*/ 1628 2027 void AdjointHorizAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/ 1629 int i, dim;2028 int i,fe_FS; 1630 2029 int* vdoflist=NULL; 1631 2030 int* pdoflist=NULL; 1632 2031 IssmDouble FSreconditioning; 1633 2032 1634 element->FindParam(&dim,DomainDimensionEnum); 2033 int domaintype,dim; 2034 element->FindParam(&domaintype,DomainTypeEnum); 2035 switch(domaintype){ 2036 case Domain2DhorizontalEnum: 2037 dim = 3; 2038 break; 2039 case Domain2DverticalEnum: 2040 dim = 2; 2041 break; 2042 case Domain3DEnum: 2043 dim = 3; 2044 break; 2045 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2046 } 2047 1635 2048 1636 2049 /*Fetch number of nodes and dof for this finite element*/ … … 1665 2078 /*fill in all arrays: */ 1666 2079 for(i=0;i<vnumnodes;i++){ 1667 lambdax[i] = values[i*NDOF3+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector"); 1668 lambday[i] = values[i*NDOF3+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector"); 1669 lambdaz[i] = values[i*NDOF3+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector"); 2080 lambdax[i] = values[i*dim+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector"); 2081 lambday[i] = values[i*dim+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector"); 2082 if(dim==3){ 2083 lambdaz[i] = values[i*dim+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector"); 2084 } 1670 2085 } 1671 2086 for(i=0;i<pnumnodes;i++){ … … 1680 2095 element->AddInput(AdjointxEnum,lambdax,element->VelocityInterpolation()); 1681 2096 element->AddInput(AdjointyEnum,lambday,element->VelocityInterpolation()); 1682 element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation()); 1683 element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation()); 2097 if(domaintype!=Domain2DverticalEnum) element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation()); 2098 2099 element->FindParam(&fe_FS,FlowequationFeFSEnum); 2100 if(fe_FS!=LATaylorHoodEnum && fe_FS!=LACrouzeixRaviartEnum) 2101 element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation()); 1684 2102 1685 2103 /*Free ressources:*/
Note:
See TracChangeset
for help on using the changeset viewer.
