Changeset 4839


Ignore:
Timestamp:
07/27/10 15:58:59 (15 years ago)
Author:
Mathieu Morlighem
Message:

Nicer way to compute adjoint

Location:
issm/trunk/src
Files:
2 added
2 deleted
23 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r4785 r4839  
    409409                                        ./modules/OutputResultsx/ElementResultsToPatch.cpp\
    410410                                        ./modules/OutputResultsx/MatlabWriteResults.cpp\
    411                                         ./modules/Dux/Dux.h\
    412                                         ./modules/Dux/Dux.cpp\
    413411                                        ./modules/MinVelx/MinVelx.h\
    414412                                        ./modules/MinVelx/MinVelx.cpp\
    … …  
    970968                                        ./modules/OutputResultsx/ElementResultsToPatch.cpp\
    971969                                        ./modules/OutputResultsx/FileWriteResults.cpp\
    972                                         ./modules/Dux/Dux.h\
    973                                         ./modules/Dux/Dux.cpp\
    974970                                        ./modules/MinVelx/MinVelx.h\
    975971                                        ./modules/MinVelx/MinVelx.cpp\
    … …  
    11491145                                        ./solutions/SolutionConfiguration.cpp\
    11501146                                        ./solvers/solver_linear.cpp\
     1147                                        ./solvers/solver_adjoint_linear.cpp\
    11511148                                        ./solvers/solver_diagnostic_nonlinear.cpp\
    11521149                                        ./solvers/solver_thermal_nonlinear.cpp\

  • issm/trunk/src/c/modules/modules.h

    r4773 r4839  
    2121#include "./CostFunctionx/CostFunctionx.h"
    2222#include "./DakotaResponsesx/DakotaResponsesx.h"
    23 #include "./Dux/Dux.h"
    2423#include "./ElementConnectivityx/ElementConnectivityx.h"
    2524#include "./FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.h"

  • issm/trunk/src/c/objects/Elements/Beam.cpp

    r4780 r4839  
    285285        }
    286286
    287 }
    288 /*}}}*/
    289 /*FUNCTION Beam::Du{{{1*/
    290 void  Beam::Du(Vec){
    291         ISSMERROR(" not supported yet!");
    292287}
    293288/*}}}*/

  • issm/trunk/src/c/objects/Elements/Element.h

    r4780 r4839  
    3737                virtual void   GetThicknessList(double* thickness_list)=0;
    3838                virtual void   GetBedList(double* bed_list)=0;
    39                 virtual void   Du(Vec du_g)=0;
    4039                virtual void   Gradj(Vec gradient,int control_type)=0;
    4140                virtual void   GradjDrag(Vec gradient)=0;

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4838 r4839  
    356356
    357357        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    358         if (analysis_type==ControlAnalysisEnum){
    359                 InputUpdateFromSolutionDiagnosticHoriz( solution);
    360         }
    361         else if (analysis_type==DiagnosticHorizAnalysisEnum){
     358        if (analysis_type==DiagnosticHorizAnalysisEnum){
    362359                InputUpdateFromSolutionDiagnosticHoriz( solution);
    363360        }
    … …  
    730727
    731728        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    732         if (analysis_type==ControlAnalysisEnum){
     729        if (analysis_type==DiagnosticHorizAnalysisEnum){
    733730                CreateKMatrixDiagnosticHoriz( Kgg);
    734731        }
    735         else if (analysis_type==DiagnosticHorizAnalysisEnum){
     732        else if (analysis_type==AdjointHorizAnalysisEnum){
     733                /*Same as diagnostic*/
    736734                CreateKMatrixDiagnosticHoriz( Kgg);
    737735        }
    … …  
    745743                CreateKMatrixDiagnosticStokes( Kgg);
    746744        }
     745        else if (analysis_type==AdjointStokesAnalysisEnum){
     746                /*Same as diagnostic*/
     747                CreateKMatrixDiagnosticStokes( Kgg);
     748        }
    747749        else if (analysis_type==BedSlopeXAnalysisEnum || analysis_type==SurfaceSlopeXAnalysisEnum || analysis_type==BedSlopeYAnalysisEnum || analysis_type==SurfaceSlopeYAnalysisEnum){
    748750                CreateKMatrixSlope( Kgg);
    … …  
    779781
    780782        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    781         if (analysis_type==ControlAnalysisEnum){
     783        if (analysis_type==DiagnosticHorizAnalysisEnum){
    782784                CreatePVectorDiagnosticHoriz( pg);
    783785        }
    784         else if (analysis_type==DiagnosticHorizAnalysisEnum){
    785                 CreatePVectorDiagnosticHoriz( pg);
     786        else if (analysis_type==AdjointHorizAnalysisEnum){
     787                CreatePVectorAdjointHoriz( pg);
    786788        }
    787789        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    794796                CreatePVectorDiagnosticStokes( pg);
    795797        }
     798        else if (analysis_type==AdjointStokesAnalysisEnum){
     799                CreatePVectorAdjointStokes( pg);
     800        }
    796801        else if (analysis_type==BedSlopeXAnalysisEnum || analysis_type==SurfaceSlopeXAnalysisEnum || analysis_type==BedSlopeYAnalysisEnum || analysis_type==SurfaceSlopeYAnalysisEnum){
    797802                CreatePVectorSlope( pg);
    … …  
    814819        else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
    815820
    816 }
    817 /*}}}*/
    818 /*FUNCTION Penta::Du {{{1*/
    819 void  Penta::Du(Vec du_g){
    820 
    821         int i;
    822         Tria* tria=NULL;
    823 
    824         /*inputs: */
    825         bool onwater;
    826         bool collapse;
    827         bool onsurface;
    828         bool onbed;
    829 
    830         /*retrieve inputs :*/
    831         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    832         inputs->GetParameterValue(&collapse,CollapseEnum);
    833         inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    834         inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
    835 
    836         /*If on water, skip: */
    837         if(onwater)return;
    838 
    839         /*Bail out if this element if:
    840          * -> Non collapsed and not on the surface
    841          * -> collapsed (2d model) and not on bed) */
    842         if ((!collapse && !onsurface) || (collapse && !onbed)){
    843                 return;
    844         }
    845         else if (collapse){
    846 
    847                 /*This element should be collapsed into a tria element at its base. Create this tria element,
    848                  * and compute Du*/
    849                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
    850                 tria->Du(du_g);
    851                 delete tria;
    852                 return;
    853         }
    854         else{
    855 
    856                 tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    857                 tria->Du(du_g);
    858                 delete tria;
    859                 return;
    860         }
    861821}
    862822/*}}}*/
    … …  
    42164176}
    42174177/*}}}*/
     4178/*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/
     4179void  Penta::CreatePVectorAdjointHoriz(Vec p_g){
     4180
     4181        int i;
     4182        Tria* tria=NULL;
     4183
     4184        /*inputs: */
     4185        bool onwater;
     4186        bool collapse;
     4187        bool onsurface;
     4188        bool onbed;
     4189
     4190        /*retrieve inputs :*/
     4191        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     4192        inputs->GetParameterValue(&collapse,CollapseEnum);
     4193        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     4194        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     4195
     4196        /*If on water, skip: */
     4197        if(onwater) return;
     4198
     4199        /*Bail out if this element if:
     4200         * -> Non collapsed and not on the surface
     4201         * -> collapsed (2d model) and not on bed) */
     4202        if ((!collapse && !onsurface) || (collapse && !onbed)){
     4203                return;
     4204        }
     4205        else if (collapse){
     4206
     4207                /*This element should be collapsed into a tria element at its base. Create this tria element,
     4208                 * and compute pe_g*/
     4209                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
     4210                tria->CreatePVectorAdjointHoriz(p_g);
     4211                delete tria;
     4212                return;
     4213        }
     4214        else{
     4215
     4216                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     4217                tria->CreatePVectorAdjointHoriz(p_g);
     4218                delete tria;
     4219                return;
     4220        }
     4221}
     4222/*}}}*/
    42184223/*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
    42194224void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
    … …  
    44974502        xfree((void**)&vert_gauss_weights);
    44984503
     4504}
     4505/*}}}*/
     4506/*FUNCTION Penta::CreatePVectorAdjointStokes{{{1*/
     4507void  Penta::CreatePVectorAdjointStokes(Vec p_g){
     4508
     4509        int i;
     4510        Tria* tria=NULL;
     4511
     4512        /*inputs: */
     4513        bool onwater;
     4514        bool onsurface;
     4515
     4516        /*retrieve inputs :*/
     4517        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     4518        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     4519
     4520        /*If on water, skip: */
     4521        if(onwater || !onsurface)return;
     4522
     4523        /*Call Tria's function*/
     4524        tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     4525        tria->CreatePVectorAdjointStokes(p_g);
     4526        delete tria;
     4527        return;
    44994528}
    45004529/*}}}*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r4838 r4839  
    7373                void   CreateKMatrix(Mat Kgg);
    7474                void   CreatePVector(Vec pg);
    75                 void   Du(Vec du_g);
    7675                void   GetBedList(double* bed_list);
    7776                void*  GetMatPar();
    … …  
    127126                void      CreatePVectorBalancedvelocities( Vec pg);
    128127                void      CreatePVectorDiagnosticHoriz( Vec pg);
     128                void      CreatePVectorAdjointHoriz( Vec pg);
    129129                void      CreatePVectorDiagnosticHutter( Vec pg);
    130130                void      CreatePVectorDiagnosticStokes( Vec pg);
     131                void      CreatePVectorAdjointStokes( Vec pg);
    131132                void      CreatePVectorDiagnosticVert( Vec pg);
    132133                void      CreatePVectorMelting( Vec pg);

  • issm/trunk/src/c/objects/Elements/Sing.cpp

    r4780 r4839  
    256256        }
    257257
    258 }
    259 /*}}}*/
    260 /*FUNCTION Sing::Du {{{1*/
    261 void  Sing::Du(Vec){
    262         ISSMERROR(" not supported yet!");
    263258}
    264259/*}}}*/

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r4836 r4839  
    364364
    365365        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    366         if (analysis_type==ControlAnalysisEnum){
    367                 InputUpdateFromSolutionDiagnosticHoriz( solution);
    368         }
    369         else if (analysis_type==DiagnosticHorizAnalysisEnum){
     366        if (analysis_type==DiagnosticHorizAnalysisEnum){
    370367                InputUpdateFromSolutionDiagnosticHoriz( solution);
    371368        }
    … …  
    664661
    665662        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    666         if (analysis_type==ControlAnalysisEnum){
     663        if (analysis_type==DiagnosticHorizAnalysisEnum){
    667664                CreateKMatrixDiagnosticHoriz( Kgg);
    668665        }
    669         else if (analysis_type==DiagnosticHorizAnalysisEnum){
     666        else if (analysis_type==AdjointHorizAnalysisEnum){
     667                /*Same as diagnostic*/
    670668                CreateKMatrixDiagnosticHoriz( Kgg);
    671                 }
     669        }
    672670        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    673671                CreateKMatrixDiagnosticHutter( Kgg);
    … …  
    715713
    716714        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    717         if (analysis_type==ControlAnalysisEnum){
     715        if (analysis_type==DiagnosticHorizAnalysisEnum){
    718716                CreatePVectorDiagnosticHoriz( pg);
    719717        }
    720         else if (analysis_type==DiagnosticHorizAnalysisEnum){
    721                 CreatePVectorDiagnosticHoriz( pg);
     718        else if (analysis_type==AdjointHorizAnalysisEnum){
     719                CreatePVectorAdjointHoriz( pg);
    722720        }
    723721        else if (analysis_type==DiagnosticHutterAnalysisEnum){
    … …  
    749747                ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
    750748        }
    751 
    752 }
    753 /*}}}*/
    754 /*FUNCTION Tria::Du {{{1*/
    755 void Tria::Du(Vec du_g){
    756 
    757         int i;
    758 
    759         /* node data: */
    760         const int    numgrids=3;
    761         const int    numdof=2*numgrids;
    762         const int    NDOF2=2;
    763         double       xyz_list[numgrids][3];
    764         int          doflist[numdof];
    765         int          numberofdofspernode;
    766         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    767 
    768         /* grid data: */
    769         double vx_list[numgrids];
    770         double vy_list[numgrids];
    771         double obs_vx_list[numgrids];
    772         double obs_vy_list[numgrids];
    773         double dux_list[numgrids];
    774         double duy_list[numgrids];
    775         double weights_list[numgrids];
    776 
    777         /* gaussian points: */
    778         int     num_gauss,ig;
    779         double* first_gauss_area_coord  =  NULL;
    780         double* second_gauss_area_coord =  NULL;
    781         double* third_gauss_area_coord  =  NULL;
    782         double* gauss_weights           =  NULL;
    783         double  gauss_weight;
    784         double  gauss_l1l2l3[3];
    785 
    786         /* parameters: */
    787         double  obs_velocity_mag,velocity_mag;
    788         double  dux,duy;
    789         double  meanvel, epsvel;
    790 
    791         /*element vector : */
    792         double  due_g[numdof]={0.0};
    793         double  due_g_gaussian[numdof];
    794 
    795         /* Jacobian: */
    796         double Jdet;
    797 
    798         /*nodal functions: */
    799         double l1l2l3[3];
    800 
    801         /*relative and algorithmic fitting: */
    802         double scalex=0;
    803         double scaley=0;
    804         double scale=0;
    805         double S=0;
    806         int    fit=-1;
    807 
    808         /*retrieve some parameters: */
    809         this->parameters->FindParam(&meanvel,MeanVelEnum);
    810         this->parameters->FindParam(&epsvel,EpsVelEnum);
    811 
    812         /* Get node coordinates and dof list: */
    813         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    814         GetDofList(&doflist[0],&numberofdofspernode);
    815 
    816         /* Recover input data: */
    817         inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
    818         inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
    819 
    820         inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
    821         inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
    822        
    823         inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
    824        
    825         inputs->GetParameterValue(&fit,FitEnum);
    826         if(fit==3){
    827                 inputs->GetParameterValue(&S,SurfaceAreaEnum);
    828         }
    829 
    830         /*Get Du at the 3 nodes (integration of the linearized function)
    831          * Here we integrate linearized functions:
    832          *               
    833          * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
    834          *
    835          *       d J                  dJ_i
    836          * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
    837          *       d u                  du_i
    838          *
    839          * where J_i are the misfits at the 3 nodes of the triangle
    840          *       Phi_i is the nodal function (P1) with respect to
    841          *       the vertex i
    842          */
    843         if(fit==0){
    844                 /*We are using an absolute misfit:
    845                  *
    846                  *      1  [           2              2 ]
    847                  * J = --- | (u - u   )  +  (v - v   )  |
    848                  *      2  [       obs            obs   ]
    849                  *
    850                  *        dJ             2
    851                  * DU = - -- = (u   - u )
    852                  *        du     obs
    853                  */
    854                 for (i=0;i<numgrids;i++){
    855                         dux_list[i]=obs_vx_list[i]-vx_list[i];
    856                         duy_list[i]=obs_vy_list[i]-vy_list[i];
    857                 }
    858         }
    859         else if(fit==1){
    860                 /*We are using a relative misfit:
    861                  *                       
    862                  *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    863                  * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    864                  *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    865                  *              obs                        obs                     
    866                  *
    867                  *        dJ     \bar{v}^2
    868                  * DU = - -- = ------------- (u   - u )
    869                  *        du   (u   + eps)^2    obs
    870                  *               obs
    871                  */
    872                 for (i=0;i<numgrids;i++){
    873                         scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
    874                         scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
    875                         if(obs_vx_list[i]==0)scalex=0;
    876                         if(obs_vy_list[i]==0)scaley=0;
    877                         dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
    878                         duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
    879                 }
    880         }
    881         else if(fit==2){
    882                 /*We are using a logarithmic misfit:
    883                  *                       
    884                  *                 [        vel + eps     ] 2
    885                  * J = 4 \bar{v}^2 | log ( -----------  ) | 
    886                  *                 [       vel   + eps    ]
    887                  *                            obs
    888                  *
    889                  *        dJ                 2 * log(...)
    890                  * DU = - -- = - 4 \bar{v}^2 -------------  u
    891                  *        du                 vel^2 + eps
    892                  *           
    893                  */
    894                 for (i=0;i<numgrids;i++){
    895                         velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
    896                         obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
    897                         scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    898                         dux_list[i]=scale*vx_list[i];
    899                         duy_list[i]=scale*vy_list[i];
    900                 }
    901         }
    902         else if(fit==3){
    903                 /*We are using an spacially average absolute misfit:
    904                  *
    905                  *      1                    2              2
    906                  * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    907                  *      S                obs            obs
    908                  *
    909                  *        dJ      1       1
    910                  * DU = - -- = - --- ----------- * 2 (u - u   )
    911                  *        du      S  2 sqrt(...)           obs
    912                  */
    913                 for (i=0;i<numgrids;i++){
    914                         scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
    915                         dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
    916                         duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
    917                 }
    918         }
    919         else if(fit==4){
    920                 /*We are using an logarithmic 2 misfit:
    921                  *
    922                  *      1            [        |u| + eps     2          |v| + eps     2  ]
    923                  * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    924                  *      2            [       |u    |+ eps              |v    |+ eps     ]
    925                  *                              obs                       obs
    926                  *        dJ                              1      u                             1
    927                  * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    928                  *        du                         |u| + eps  |u|                           u + eps
    929                  */
    930                 for (i=0;i<numgrids;i++){
    931                         dux_list[i] = - pow(meanvel,(double)2)*(
    932                                                 log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
    933                         duy_list[i] = - pow(meanvel,(double)2)*(
    934                                                 log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
    935                 }
    936         }
    937         else{
    938                 /*Not supported yet! : */
    939                 ISSMERROR("%s%g","unsupported type of fit: ",fit);
    940         }
    941 
    942         /*Apply weights to DU*/
    943         for (i=0;i<numgrids;i++){
    944                 dux_list[i]=weights_list[i]*dux_list[i];
    945                 duy_list[i]=weights_list[i]*duy_list[i];
    946         }
    947 
    948         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    949         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    950 
    951         /* Start  looping on the number of gaussian points: */
    952         for (ig=0; ig<num_gauss; ig++){
    953                 /*Pick up the gaussian point: */
    954                 gauss_weight=*(gauss_weights+ig);
    955                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    956                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    957                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    958 
    959                 /* Get Jacobian determinant: */
    960                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    961 
    962                 /* Get nodal functions value at gaussian point:*/
    963                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    964 
    965                 /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
    966 
    967                 /*Compute absolute(x/y) at gaussian point: */
    968                 GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
    969                 GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
    970 
    971                 /*compute Du*/
    972                 for (i=0;i<numgrids;i++){
    973                         due_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i];
    974                         due_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i];
    975                 }
    976 
    977                 /*Add due_g_gaussian vector to due_g: */
    978                 for( i=0; i<numdof; i++){
    979                         due_g[i]+=due_g_gaussian[i];
    980                 }
    981         }
    982 
    983         /*Add due_g to global vector du_g: */
    984         VecSetValues(du_g,numdof,doflist,(const double*)due_g,ADD_VALUES);
    985 
    986         xfree((void**)&first_gauss_area_coord);
    987         xfree((void**)&second_gauss_area_coord);
    988         xfree((void**)&third_gauss_area_coord);
    989         xfree((void**)&gauss_weights);
    990749
    991750}
    … …  
    43894148}
    43904149/*}}}*/
     4150/*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/
     4151void Tria::CreatePVectorAdjointHoriz(Vec p_g){
     4152
     4153        int i;
     4154
     4155        /* node data: */
     4156        const int    numgrids=3;
     4157        const int    numdof=2*numgrids;
     4158        const int    NDOF2=2;
     4159        double       xyz_list[numgrids][3];
     4160        int          doflist[numdof];
     4161        int          numberofdofspernode;
     4162        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     4163
     4164        /* grid data: */
     4165        double vx_list[numgrids];
     4166        double vy_list[numgrids];
     4167        double obs_vx_list[numgrids];
     4168        double obs_vy_list[numgrids];
     4169        double dux_list[numgrids];
     4170        double duy_list[numgrids];
     4171        double weights_list[numgrids];
     4172
     4173        /* gaussian points: */
     4174        int     num_gauss,ig;
     4175        double* first_gauss_area_coord  =  NULL;
     4176        double* second_gauss_area_coord =  NULL;
     4177        double* third_gauss_area_coord  =  NULL;
     4178        double* gauss_weights           =  NULL;
     4179        double  gauss_weight;
     4180        double  gauss_l1l2l3[3];
     4181
     4182        /* parameters: */
     4183        double  obs_velocity_mag,velocity_mag;
     4184        double  dux,duy;
     4185        double  meanvel, epsvel;
     4186
     4187        /*element vector : */
     4188        double  pe_g[numdof]={0.0};
     4189        double  pe_g_gaussian[numdof];
     4190
     4191        /* Jacobian: */
     4192        double Jdet;
     4193
     4194        /*nodal functions: */
     4195        double l1l2l3[3];
     4196
     4197        /*relative and algorithmic fitting: */
     4198        double scalex=0;
     4199        double scaley=0;
     4200        double scale=0;
     4201        double S=0;
     4202        int    fit=-1;
     4203
     4204        /*retrieve some parameters: */
     4205        this->parameters->FindParam(&meanvel,MeanVelEnum);
     4206        this->parameters->FindParam(&epsvel,EpsVelEnum);
     4207
     4208        /* Get node coordinates and dof list: */
     4209        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4210        GetDofList(&doflist[0],&numberofdofspernode);
     4211
     4212        /* Recover input data: */
     4213        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
     4214        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
     4215
     4216        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
     4217        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
     4218
     4219        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
     4220
     4221        inputs->GetParameterValue(&fit,FitEnum);
     4222        if(fit==3){
     4223                inputs->GetParameterValue(&S,SurfaceAreaEnum);
     4224        }
     4225
     4226        /*Get Du at the 3 nodes (integration of the linearized function)
     4227         * Here we integrate linearized functions:
     4228         *               
     4229         * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
     4230         *
     4231         *       d J                  dJ_i
     4232         * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
     4233         *       d u                  du_i
     4234         *
     4235         * where J_i are the misfits at the 3 nodes of the triangle
     4236         *       Phi_i is the nodal function (P1) with respect to
     4237         *       the vertex i
     4238         */
     4239        if(fit==0){
     4240                /*We are using an absolute misfit:
     4241                 *
     4242                 *      1  [           2              2 ]
     4243                 * J = --- | (u - u   )  +  (v - v   )  |
     4244                 *      2  [       obs            obs   ]
     4245                 *
     4246                 *        dJ             2
     4247                 * DU = - -- = (u   - u )
     4248                 *        du     obs
     4249                 */
     4250                for (i=0;i<numgrids;i++){
     4251                        dux_list[i]=obs_vx_list[i]-vx_list[i];
     4252                        duy_list[i]=obs_vy_list[i]-vy_list[i];
     4253                }
     4254        }
     4255        else if(fit==1){
     4256                /*We are using a relative misfit:
     4257                 *                       
     4258                 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     4259                 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     4260                 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     4261                 *              obs                        obs                     
     4262                 *
     4263                 *        dJ     \bar{v}^2
     4264                 * DU = - -- = ------------- (u   - u )
     4265                 *        du   (u   + eps)^2    obs
     4266                 *               obs
     4267                 */
     4268                for (i=0;i<numgrids;i++){
     4269                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
     4270                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
     4271                        if(obs_vx_list[i]==0)scalex=0;
     4272                        if(obs_vy_list[i]==0)scaley=0;
     4273                        dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
     4274                        duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
     4275                }
     4276        }
     4277        else if(fit==2){
     4278                /*We are using a logarithmic misfit:
     4279                 *                       
     4280                 *                 [        vel + eps     ] 2
     4281                 * J = 4 \bar{v}^2 | log ( -----------  ) | 
     4282                 *                 [       vel   + eps    ]
     4283                 *                            obs
     4284                 *
     4285                 *        dJ                 2 * log(...)
     4286                 * DU = - -- = - 4 \bar{v}^2 -------------  u
     4287                 *        du                 vel^2 + eps
     4288                 *           
     4289                 */
     4290                for (i=0;i<numgrids;i++){
     4291                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
     4292                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
     4293                        scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
     4294                        dux_list[i]=scale*vx_list[i];
     4295                        duy_list[i]=scale*vy_list[i];
     4296                }
     4297        }
     4298        else if(fit==3){
     4299                /*We are using an spacially average absolute misfit:
     4300                 *
     4301                 *      1                    2              2
     4302                 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     4303                 *      S                obs            obs
     4304                 *
     4305                 *        dJ      1       1
     4306                 * DU = - -- = - --- ----------- * 2 (u - u   )
     4307                 *        du      S  2 sqrt(...)           obs
     4308                 */
     4309                for (i=0;i<numgrids;i++){
     4310                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
     4311                        dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
     4312                        duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
     4313                }
     4314        }
     4315        else if(fit==4){
     4316                /*We are using an logarithmic 2 misfit:
     4317                 *
     4318                 *      1            [        |u| + eps     2          |v| + eps     2  ]
     4319                 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     4320                 *      2            [       |u    |+ eps              |v    |+ eps     ]
     4321                 *                              obs                       obs
     4322                 *        dJ                              1      u                             1
     4323                 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
     4324                 *        du                         |u| + eps  |u|                           u + eps
     4325                 */
     4326                for (i=0;i<numgrids;i++){
     4327                        dux_list[i] = - pow(meanvel,(double)2)*(
     4328                                                log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
     4329                        duy_list[i] = - pow(meanvel,(double)2)*(
     4330                                                log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
     4331                }
     4332        }
     4333        else{
     4334                /*Not supported yet! : */
     4335                ISSMERROR("%s%g","unsupported type of fit: ",fit);
     4336        }
     4337
     4338        /*Apply weights to DU*/
     4339        for (i=0;i<numgrids;i++){
     4340                dux_list[i]=weights_list[i]*dux_list[i];
     4341                duy_list[i]=weights_list[i]*duy_list[i];
     4342        }
     4343
     4344        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     4345        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4346
     4347        /* Start  looping on the number of gaussian points: */
     4348        for (ig=0; ig<num_gauss; ig++){
     4349                /*Pick up the gaussian point: */
     4350                gauss_weight=*(gauss_weights+ig);
     4351                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     4352                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     4353                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4354
     4355                /* Get Jacobian determinant: */
     4356                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4357
     4358                /* Get nodal functions value at gaussian point:*/
     4359                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4360
     4361                /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
     4362
     4363                /*Compute absolute(x/y) at gaussian point: */
     4364                GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
     4365                GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
     4366
     4367                /*compute Du*/
     4368                for (i=0;i<numgrids;i++){
     4369                        pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i];
     4370                        pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i];
     4371                }
     4372
     4373                /*Add pe_g_gaussian vector to pe_g: */
     4374                for( i=0; i<numdof; i++){
     4375                        pe_g[i]+=pe_g_gaussian[i];
     4376                }
     4377        }
     4378
     4379        /*Add pe_g to global vector p_g: */
     4380        VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     4381
     4382        /*Clean up*/
     4383        xfree((void**)&first_gauss_area_coord);
     4384        xfree((void**)&second_gauss_area_coord);
     4385        xfree((void**)&third_gauss_area_coord);
     4386        xfree((void**)&gauss_weights);
     4387
     4388}
     4389/*}}}*/
     4390/*FUNCTION Tria::CreatePVectorAdjointStokes{{{1*/
     4391void Tria::CreatePVectorAdjointStokes(Vec p_g){
     4392
     4393        int i;
     4394
     4395        /* node data: */
     4396        const int    numgrids=3;
     4397        const int    numdof=2*numgrids;
     4398        const int    NDOF2=2;
     4399        double       xyz_list[numgrids][3];
     4400        int          doflist[numdof];
     4401        int          numberofdofspernode;
     4402        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     4403
     4404        /* grid data: */
     4405        double vx_list[numgrids];
     4406        double vy_list[numgrids];
     4407        double obs_vx_list[numgrids];
     4408        double obs_vy_list[numgrids];
     4409        double dux_list[numgrids];
     4410        double duy_list[numgrids];
     4411        double weights_list[numgrids];
     4412
     4413        /* gaussian points: */
     4414        int     num_gauss,ig;
     4415        double* first_gauss_area_coord  =  NULL;
     4416        double* second_gauss_area_coord =  NULL;
     4417        double* third_gauss_area_coord  =  NULL;
     4418        double* gauss_weights           =  NULL;
     4419        double  gauss_weight;
     4420        double  gauss_l1l2l3[3];
     4421
     4422        /* parameters: */
     4423        double  obs_velocity_mag,velocity_mag;
     4424        double  dux,duy;
     4425        double  meanvel, epsvel;
     4426
     4427        /*element vector : */
     4428        double  pe_g[numdof]={0.0};
     4429        double  pe_g_gaussian[numdof];
     4430
     4431        /* Jacobian: */
     4432        double Jdet;
     4433
     4434        /*nodal functions: */
     4435        double l1l2l3[3];
     4436
     4437        /*relative and algorithmic fitting: */
     4438        double scalex=0;
     4439        double scaley=0;
     4440        double scale=0;
     4441        double S=0;
     4442        int    fit=-1;
     4443
     4444        /*retrieve some parameters: */
     4445        this->parameters->FindParam(&meanvel,MeanVelEnum);
     4446        this->parameters->FindParam(&epsvel,EpsVelEnum);
     4447
     4448        /* Get node coordinates and dof list: */
     4449        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4450        GetDofList(&doflist[0],&numberofdofspernode);
     4451
     4452        /* Recover input data: */
     4453        inputs->GetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum);
     4454        inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum);
     4455
     4456        inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum);
     4457        inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum);
     4458
     4459        inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum);
     4460
     4461        inputs->GetParameterValue(&fit,FitEnum);
     4462        if(fit==3){
     4463                inputs->GetParameterValue(&S,SurfaceAreaEnum);
     4464        }
     4465
     4466        /*Get Du at the 3 nodes (integration of the linearized function)
     4467         * Here we integrate linearized functions:
     4468         *               
     4469         * J(E) = int_E   sum_{i=1}^3  J_i Phi_i
     4470         *
     4471         *       d J                  dJ_i
     4472         * DU= - --- = sum_{i=1}^3  - ---  Phi_i = sum_{i=1}^3 DU_i Phi_i
     4473         *       d u                  du_i
     4474         *
     4475         * where J_i are the misfits at the 3 nodes of the triangle
     4476         *       Phi_i is the nodal function (P1) with respect to
     4477         *       the vertex i
     4478         */
     4479        if(fit==0){
     4480                /*We are using an absolute misfit:
     4481                 *
     4482                 *      1  [           2              2 ]
     4483                 * J = --- | (u - u   )  +  (v - v   )  |
     4484                 *      2  [       obs            obs   ]
     4485                 *
     4486                 *        dJ             2
     4487                 * DU = - -- = (u   - u )
     4488                 *        du     obs
     4489                 */
     4490                for (i=0;i<numgrids;i++){
     4491                        dux_list[i]=obs_vx_list[i]-vx_list[i];
     4492                        duy_list[i]=obs_vy_list[i]-vy_list[i];
     4493                }
     4494        }
     4495        else if(fit==1){
     4496                /*We are using a relative misfit:
     4497                 *                       
     4498                 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
     4499                 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
     4500                 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
     4501                 *              obs                        obs                     
     4502                 *
     4503                 *        dJ     \bar{v}^2
     4504                 * DU = - -- = ------------- (u   - u )
     4505                 *        du   (u   + eps)^2    obs
     4506                 *               obs
     4507                 */
     4508                for (i=0;i<numgrids;i++){
     4509                        scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2);
     4510                        scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2);
     4511                        if(obs_vx_list[i]==0)scalex=0;
     4512                        if(obs_vy_list[i]==0)scaley=0;
     4513                        dux_list[i]=scalex*(obs_vx_list[i]-vx_list[i]);
     4514                        duy_list[i]=scaley*(obs_vy_list[i]-vy_list[i]);
     4515                }
     4516        }
     4517        else if(fit==2){
     4518                /*We are using a logarithmic misfit:
     4519                 *                       
     4520                 *                 [        vel + eps     ] 2
     4521                 * J = 4 \bar{v}^2 | log ( -----------  ) | 
     4522                 *                 [       vel   + eps    ]
     4523                 *                            obs
     4524                 *
     4525                 *        dJ                 2 * log(...)
     4526                 * DU = - -- = - 4 \bar{v}^2 -------------  u
     4527                 *        du                 vel^2 + eps
     4528                 *           
     4529                 */
     4530                for (i=0;i<numgrids;i++){
     4531                        velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil.
     4532                        obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil.
     4533                        scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
     4534                        dux_list[i]=scale*vx_list[i];
     4535                        duy_list[i]=scale*vy_list[i];
     4536                }
     4537        }
     4538        else if(fit==3){
     4539                /*We are using an spacially average absolute misfit:
     4540                 *
     4541                 *      1                    2              2
     4542                 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
     4543                 *      S                obs            obs
     4544                 *
     4545                 *        dJ      1       1
     4546                 * DU = - -- = - --- ----------- * 2 (u - u   )
     4547                 *        du      S  2 sqrt(...)           obs
     4548                 */
     4549                for (i=0;i<numgrids;i++){
     4550                        scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel);
     4551                        dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]);
     4552                        duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]);
     4553                }
     4554        }
     4555        else if(fit==4){
     4556                /*We are using an logarithmic 2 misfit:
     4557                 *
     4558                 *      1            [        |u| + eps     2          |v| + eps     2  ]
     4559                 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
     4560                 *      2            [       |u    |+ eps              |v    |+ eps     ]
     4561                 *                              obs                       obs
     4562                 *        dJ                              1      u                             1
     4563                 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
     4564                 *        du                         |u| + eps  |u|                           u + eps
     4565                 */
     4566                for (i=0;i<numgrids;i++){
     4567                        dux_list[i] = - pow(meanvel,(double)2)*(
     4568                                                log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel));
     4569                        duy_list[i] = - pow(meanvel,(double)2)*(
     4570                                                log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel));
     4571                }
     4572        }
     4573        else{
     4574                /*Not supported yet! : */
     4575                ISSMERROR("%s%g","unsupported type of fit: ",fit);
     4576        }
     4577
     4578        /*Apply weights to DU*/
     4579        for (i=0;i<numgrids;i++){
     4580                dux_list[i]=weights_list[i]*dux_list[i];
     4581                duy_list[i]=weights_list[i]*duy_list[i];
     4582        }
     4583
     4584        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     4585        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     4586
     4587        /* Start  looping on the number of gaussian points: */
     4588        for (ig=0; ig<num_gauss; ig++){
     4589                /*Pick up the gaussian point: */
     4590                gauss_weight=*(gauss_weights+ig);
     4591                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     4592                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     4593                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4594
     4595                /* Get Jacobian determinant: */
     4596                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4597
     4598                /* Get nodal functions value at gaussian point:*/
     4599                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4600
     4601                /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
     4602
     4603                /*Compute absolute(x/y) at gaussian point: */
     4604                GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
     4605                GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
     4606
     4607                /*compute Du*/
     4608                for (i=0;i<numgrids;i++){
     4609                        pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i];
     4610                        pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i];
     4611                }
     4612
     4613                /*Add pe_g_gaussian vector to pe_g: */
     4614                for( i=0; i<numdof; i++){
     4615                        pe_g[i]+=pe_g_gaussian[i];
     4616                }
     4617        }
     4618
     4619        /*Add pe_g to global vector p_g: */
     4620        VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     4621
     4622        /*Clean up*/
     4623        xfree((void**)&first_gauss_area_coord);
     4624        xfree((void**)&second_gauss_area_coord);
     4625        xfree((void**)&third_gauss_area_coord);
     4626        xfree((void**)&gauss_weights);
     4627
     4628}
     4629/*}}}*/
    43914630/*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
    43924631void  Tria::CreatePVectorDiagnosticHutter(Vec pg){

  • issm/trunk/src/c/objects/Elements/Tria.h

    r4835 r4839  
    7070                void   CreateKMatrix(Mat Kgg);
    7171                void   CreatePVector(Vec pg);
    72                 void   Du(Vec du_g);
    7372                void   GetBedList(double* bed_list);
    7473                void*  GetMatPar();
    … …  
    128127                void      CreatePVectorDiagnosticBaseVert(Vec pg);
    129128                void      CreatePVectorDiagnosticHoriz(Vec pg);
     129                void      CreatePVectorAdjointHoriz(Vec pg);
     130                void      CreatePVectorAdjointStokes(Vec pg);
    130131                void      CreatePVectorDiagnosticHutter(Vec pg);
    131132                void      CreatePVectorPrognostic_CG(Vec pg);

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r4684 r4839  
    312312
    313313        /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
    314         if (analysis_type==ControlAnalysisEnum){
    315                 CreatePVectorDiagnosticHoriz( pg);
    316         }
    317         else if (analysis_type==DiagnosticHorizAnalysisEnum){
     314        if (analysis_type==DiagnosticHorizAnalysisEnum){
    318315                CreatePVectorDiagnosticHoriz( pg);
    319316        }

  • issm/trunk/src/c/solutions/adjoint_core.cpp

    r4475 r4839  
    1414void adjoint_core(FemModel* femmodel){
    1515       
    16        
    1716        /*parameters: */
    18         int   verbose=0;
    19         char* solverstring=NULL;
    20         bool  isstokes=false;
    21         bool  conserve_loads=true;
    22         int   dim;
    23         int   solution_type;
    24        
    25         /*intermediary: */
    26         Vec u_g=NULL;
    27         Mat K_ff0=NULL;
    28         Mat K_fs0=NULL;
    29 
    30         Vec du_g=NULL;
    31         Vec du_f=NULL;
    32 
    33         Vec adjoint_f=NULL;
    34         Vec adjoint_g=NULL;
     17        int  verbose        = 0;
     18        bool isstokes;
     19        bool conserve_loads = true;
     20        int  solution_type;
    3521
    3622        /*retrieve parameters:*/
    37         femmodel->parameters->FindParam(&solverstring,SolverStringEnum);
    3823        femmodel->parameters->FindParam(&verbose,VerboseEnum);
    3924        femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
    40         femmodel->parameters->FindParam(&dim,DimEnum);
    4125        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    4226
    4327        /*set analysis type to compute velocity: */
     28        _printf_("%s\n","      computing velocities");
    4429        if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
    4530        else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
    46        
    47         _printf_("%s\n","      recover solution for this stiffness and right hand side:");
    48         solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0, femmodel,conserve_loads);
    49 
    50         _printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
    51         Dux(&du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    52 
    53         _printf_("%s\n","      reduce adjoint load from g-set to f-set:");
    54         Reduceloadfromgtofx(&du_f, du_g, femmodel->Gmn, K_fs0, femmodel->ys, femmodel->nodesets,true); //true means that ys0 flag is activated: all spcs show 0 displacement
    55        
    56         /*free some ressources: */
    57         VecFree(&du_g);MatFree(&K_fs0);
    58 
    59         _printf_("%s\n","      solve for adjoint vector:");
    60         Solverx(&adjoint_f, K_ff0, du_f, NULL, solverstring);
    61        
    62         /*free some ressources: */
    63         VecFree(&du_f); MatFree(&K_ff0);
    64        
    65         _printf_("%s\n","      merge back to g set:");
    66         Mergesolutionfromftogx(&adjoint_g, adjoint_f,femmodel->Gmn,femmodel->ys,femmodel->nodesets,true);//true means that ys0 flag is activated: all spc are 0
    67        
    68         /*free some ressources: */
    69         VecFree(&adjoint_f);
     31        solver_diagnostic_nonlinear(femmodel,conserve_loads);
    7032
    7133        /*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */
     34        _printf_("%s\n","      computing adjoint");
    7235        if(isstokes)femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum);
    7336        else femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum);
     37        solver_adjoint_linear(femmodel);
    7438       
    75         InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g);
    76 
    7739        if(verbose)_printf_("saving results:\n");
    7840        if(solution_type==AdjointSolutionEnum){
    7941                InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointxEnum);
    8042                InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointyEnum);
    81                 if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointzEnum);
     43                if (isstokes){
     44                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointzEnum);
     45                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointpEnum);
     46                }
    8247        }
    83        
    84         /*Free ressources:*/
    85         xfree((void**)&solverstring);
    86         VecFree(&u_g);
    87         MatFree(&K_ff0);
    88         MatFree(&K_fs0);
    89         VecFree(&du_g);
    90         VecFree(&du_f);
    91         VecFree(&adjoint_f);
    92         VecFree(&adjoint_g);
    93 
    9448}

  • issm/trunk/src/c/solutions/diagnostic_core.cpp

    r4837 r4839  
    1515
    1616        /*parameters: */
    17         int  verbose=0;
    18         bool qmu_analysis=false;
    19         int  dim=-1;
    20         bool ishutter=false;
    21         bool ismacayealpattyn=false;
    22         bool isstokes=false;
     17        int    verbose              = 0;
     18        bool   qmu_analysis         = false;
     19        int    dim                  = -1;
     20        bool   ishutter             = false;
     21        bool   ismacayealpattyn     = false;
     22        bool   isstokes             = false;
    2323        double stokesreconditioning;
    24         bool conserve_loads=true;
    25         bool modify_loads=true;
    26         int solution_type;
     24        bool   conserve_loads       = true;
     25        bool   modify_loads         = true;
     26        int    solution_type;
    2727
    2828        /* recover parameters:*/
    … …  
    6565                if(verbose)_printf_("%s\n"," computing horizontal velocities...");
    6666                femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
    67                 solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,modify_loads);
     67                solver_diagnostic_nonlinear(femmodel,modify_loads);
    6868        }
    6969       
    … …  
    8282                        if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
    8383                        femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
    84                         solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
     84                        solver_diagnostic_nonlinear(femmodel,conserve_loads);
    8585                }
    8686        }

  • issm/trunk/src/c/solutions/objectivefunctionC.cpp

    r4560 r4839  
    6969
    7070        /*Run diagnostic with updated inputs: */
    71         if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);  //true means we conserve loads at each diagnostic run
     71        if(!control_steady) solver_diagnostic_nonlinear(femmodel,conserve_loads);  //true means we conserve loads at each diagnostic run
    7272        else                diagnostic_core(femmodel);  //We need a 3D velocity!! (vz is required for the next thermal run)
    7373

  • issm/trunk/src/c/solutions/stokescontrolinit.cpp

    r4837 r4839  
    4141        /*Run a complete diagnostic to update the Stokes spcs: */
    4242        femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
    43         solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
     43        solver_diagnostic_nonlinear(femmodel,conserve_loads);
    4444        femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum);
    4545        solver_linear(femmodel);
    … …  
    5050        if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
    5151        femmodel->SetCurrentConfiguration(DiagnosticStokesAnalysisEnum);
    52         solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
     52        solver_diagnostic_nonlinear(femmodel,conserve_loads);
    5353}

  • issm/trunk/src/c/solutions/thermal_core_step.cpp

    r4837 r4839  
    2121        if(verbose)_printf_("computing temperatures:\n");
    2222        femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
    23         solver_thermal_nonlinear(NULL,NULL,femmodel);
     23        solver_thermal_nonlinear(femmodel);
    2424
    2525        if(verbose)_printf_("computing melting:\n");

  • issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp

    r4832 r4839  
    1010#include "./solvers.h"
    1111
    12 void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* femmodel,bool conserve_loads){
     12void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads){
    1313
    1414
    … …  
    127127                }
    128128        }
    129        
    130         //more output might be needed, when running in control
    131         if(pKff0){
    132 
    133                 kflag=1; pflag=0; //stiffness generation only
    134        
    135                 SystemMatricesx(&Kgg, &pg,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters,kflag,pflag);
    136                 Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->Gmn,femmodel->nodesets);
    137                 MatFree(&Kgg);VecFree(&pg);
    138 
    139         }
    140129
    141130        /*Delete loads only if no ouput was requested: */
    … …  
    144133        /*clean up*/
    145134        VecFree(&uf);
     135        VecFree(&ug);
    146136        VecFree(&old_uf);
    147137        VecFree(&old_ug);
    148138        xfree((void**)&solver_string);
    149139       
    150         /*Assign output pointers: */
    151         if(pug)*pug=ug;
    152         else VecFree(&ug);
    153 
    154         if(pKff0)*pKff0=Kff;
    155         if(pKfs0)*pKfs0=Kfs;
    156 
    157140}

  • issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp

    r4524 r4839  
    88#include "../modules/modules.h"
    99
    10 void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem){
     10void solver_thermal_nonlinear(FemModel* fem){
    1111
    1212        /*solution : */
    … …  
    125125        MatFree(&Kgg_nopenalty);
    126126        VecFree(&pg_nopenalty);
     127        VecFree(&tg);
    127128        VecFree(&tf);
    128129        VecFree(&tf_old);
    129130        xfree((void**)&solver_string);
    130 
    131         /*Assign output pointers: */
    132         if(ptg)*ptg=tg;
    133         else VecFree(&tg);
    134         if(pmelting_offset)*pmelting_offset=melting_offset;
    135131}

  • issm/trunk/src/c/solvers/solvers.h

    r4837 r4839  
    1212class FemModel;
    1313
    14 void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* femmodel);
    15 void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* femmodel,bool conserve_loads);
     14void solver_thermal_nonlinear(FemModel* femmodel);
     15void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads);
    1616void solver_linear(FemModel* femmodel);
     17void solver_adjoint_linear(FemModel* femmodel);
    1718
    1819#endif

  • issm/trunk/src/m/solutions/adjoint_core.m

    r4687 r4839  
    1010
    1111        %set analysis type to compute velocity:
     12        displaystring('\n%s',['      computing velocities']);
    1213        if(isstokes),
    1314                femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum);
    … …  
    1516                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
    1617        end
     18        femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads);
    1719
    18         displaystring('\n%s',['      recover solution for this stiffness and right hand side:']);
    19         [femmodel,ug,K_ff0,K_fs0]=solver_diagnostic_nonlinear(femmodel,conserve_loads);
    20 
    21         displaystring('\n%s',['      Build Du, difference between observed and modeled velocities:']);
    22         du_g=Du(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    23 
    24         displaystring('\n%s',['      reduce adjoint load from g_set to f_set:']);
    25         du_f=Reduceloadfromgtof(du_g,femmodel.Gmn,K_fs0,femmodel.ys,femmodel.nodesets,true); %true means that ys0 flag is activated: all spcs show 0 displacement
    26 
    27         displaystring('\n%s',['      solve for adjoint vector:']);
    28         adjoint_f=Solver(K_ff0,du_f,[],femmodel.parameters);
    29 
    30         displaystring('\n%s',['      merge back to g set:']);
    31         adjoint_g=Mergesolutionfromftog(adjoint_f,femmodel.Gmn,femmodel.ys,femmodel.nodesets,true);%true means that ys0 flag is activated: all spc are 0
    32 
    33         %Update inputs using adjoint solution, and same type of setup as diagnostic solution:
     20        displaystring('\n%s',['      computing asjoint']);
    3421        if(isstokes),
    3522                femmodel=SetCurrentConfiguration(femmodel,DiagnosticStokesAnalysisEnum,AdjointStokesAnalysisEnum);
    … …  
    3724                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum);
    3825        end
    39 
    40         [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,adjoint_g);
     26        femmodel=solver_adjoint_linear(femmodel);
    4127
    4228        displaystring(verbose,'\n%s',['      saving results...']);
    … …  
    4430                femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointxEnum);
    4531                femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointyEnum);
    46                 if(dim==3),
     32                if(isstokes),
    4733                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointzEnum);
     34                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,AdjointpEnum);
    4835                end
    4936        end

  • issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m

    r4833 r4839  
    1 function [femmodel ug varargout]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
     1function femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads)
    22%SOLVER_DIAGNOSTIC_NONLINEAR - core solver of diagnostic run
    33%
    44%   Usage:
    5 %      [femmodel ug varargout]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
     5%      [femmodel]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
    66       
    77
    … …  
    8282        end
    8383
    84         %more output might be needed, when running in control_core.m
    85         nout=max(nargout,1)-2;
    86         if nout==2,
    87                 femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=0;
    88                 [K_gg, p_g]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
    89                 [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets);
    90                 varargout(1)={K_ff};
    91                 varargout(2)={K_fs};
    92         end
    93 
    94 
    9584        %deal with loads:
    9685        if conserve_loads==false,

  • issm/trunk/src/m/solvers/solver_linear.m

    r4687 r4839  
    1 function [femmodel u_g]=solver_linear(femmodel)
     1function femmodel=solver_linear(femmodel)
    22%SOLVER_LINEAR - core solver of any linear solution sequence
    33%
    44%   Usage:
    5 %      [femmodel, u_g]=solver_linear(femmodel)
     5%      femmodel =solver_linear(femmodel)
    66
    77        %stiffness and load generation only:

  • issm/trunk/src/m/solvers/solver_thermal_nonlinear.m

    r4687 r4839  
    1 function [femmodel t_g ,melting_offset]=solver_thermal_nonlinear(femmodel)
     1function femmodel=solver_thermal_nonlinear(femmodel)
    22%SOLVER_THERMAL_NONLINEAR - core of thermal solution sequence.
    33%   femmodel is returned together with temperature and melting_offset, in case loads have been modified
    44%
    55%   Usage:
    6 %      [femmodel t_g melting_offset]=solver_thermal_nonlinear(femmodel)
     6%      [femmodel]=solver_thermal_nonlinear(femmodel)
    77
    88        count=1;

  • issm/trunk/src/mex/Makefile.am

    r4773 r4839  
    1919                                CostFunction \
    2020                                DakotaResponses\
    21                                 Du\
    2221                                Echo\
    2322                                ElementConnectivity\
    … …  
    157156                          DakotaResponses/DakotaResponses.h
    158157
    159 Du_SOURCES = Du/Du.cpp\
    160                           Du/Du.h
    161 
    162158Echo_SOURCES = Echo/Echo.cpp\
    163159                          Echo/Echo.h
Note: See TracChangeset for help on using the changeset viewer.