Ignore:
Timestamp:
05/09/14 15:37:18 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: pushing misfits away from elements

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r17907 r17976  
    11341134void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/
    11351135
    1136         IssmDouble J = 0;
     1136        /*output: */
     1137        IssmDouble J=0.;
    11371138        IssmDouble J_sum;
    11381139
    1139         for(int i=0;i<this->elements->Size();i++){
    1140                 Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
    1141                 J+=element->BalancethicknessMisfit();
    1142         }
     1140        IssmDouble  weight,vx,vy,H,dvx[2],dvy[2],dH[2];
     1141        IssmDouble  temp,Jdet,dhdt,basal_melting,surface_mass_balance;
     1142        IssmDouble* xyz_list = NULL;
     1143        IssmDouble  dp[3];
     1144
     1145        /*Compute Misfit: */
     1146        for(int i=0;i<elements->Size();i++){
     1147                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1148
     1149                /*If on water, return 0: */
     1150                if(!element->IsIceInElement()) continue;
     1151
     1152                /* Get node coordinates*/
     1153                element->GetVerticesCoordinates(&xyz_list);
     1154                Input* weights_input              = element->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
     1155                Input* thickness_input            = element->GetInput(ThicknessEnum); _assert_(thickness_input);
     1156                Input* vx_input                   = element->GetInput(VxEnum);                                 _assert_(vx_input);
     1157                Input* vy_input                   = element->GetInput(VyEnum);                                 _assert_(vy_input);
     1158                Input* surface_mass_balance_input = element->GetInput(SurfaceforcingsMassBalanceEnum);         _assert_(surface_mass_balance_input);
     1159                Input* basal_melting_input        = element->GetInput(BasalforcingsMeltingRateEnum);           _assert_(basal_melting_input);
     1160                Input* dhdt_input                 = element->GetInput(BalancethicknessThickeningRateEnum);     _assert_(dhdt_input);
     1161
     1162                /* Start  looping on the number of gaussian points: */
     1163                Gauss* gauss=element->NewGauss(2);
     1164                for(int ig=gauss->begin();ig<gauss->end();ig++){
     1165
     1166                        gauss->GaussPoint(ig);
     1167
     1168                        /* Get Jacobian determinant: */
     1169                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1170
     1171                        /*Get all parameters at gaussian point*/
     1172                        weights_input->GetInputValue(&weight,gauss,BalancethicknessMisfitEnum);
     1173                        thickness_input->GetInputValue(&H, gauss);
     1174                        thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss);
     1175                        surface_mass_balance_input->GetInputValue(&surface_mass_balance,gauss);
     1176                        basal_melting_input->GetInputValue(&basal_melting,gauss);
     1177                        dhdt_input->GetInputValue(&dhdt,gauss);
     1178                        vx_input->GetInputValue(&vx,gauss);
     1179                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     1180                        vy_input->GetInputValue(&vy,gauss);
     1181                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     1182
     1183                        /*Balance thickness soft constraint J = 1/2 (div(Hv)-a)^2*/
     1184                        temp  = vx*dH[0]+vy*dH[1]+H*(dvx[0]+dvy[1]) - (surface_mass_balance-basal_melting-dhdt);
     1185                        J    +=weight*1/2*temp*temp*Jdet*gauss->weight;
     1186                }
     1187
     1188                /*clean up and Return: */
     1189                xDelete<IssmDouble>(xyz_list);
     1190                delete gauss;
     1191        }
     1192
     1193        /*Sum all J from all cpus of the cluster:*/
    11431194        ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    11441195        ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     
    11511202void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/
    11521203
    1153         /*Intermediary*/
    1154         int i;
    1155         Element* element=NULL;
    1156 
    11571204        /*output: */
    1158         IssmDouble J=0;
     1205        IssmDouble J=0.;
    11591206        IssmDouble J_sum;
    11601207
     1208        IssmDouble  thickness,thicknessobs,weight;
     1209        IssmDouble  Jdet;
     1210        IssmDouble* xyz_list = NULL;
     1211        IssmDouble  dp[3];
     1212
    11611213        /*Compute Misfit: */
    1162         for (i=0;i<elements->Size();i++){
    1163                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
    1164                 J+=element->ThicknessAbsGradient();
     1214        for(int i=0;i<elements->Size();i++){
     1215                Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     1216
     1217                /*If on water, return 0: */
     1218                if(!element->IsIceInElement()) continue;
     1219
     1220                /* Get node coordinates*/
     1221                element->GetVerticesCoordinates(&xyz_list);
     1222
     1223                /*Retrieve all inputs we will be needing: */
     1224                Input* weights_input   =element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     1225                Input* thickness_input =element->GetInput(ThicknessEnum);                          _assert_(thickness_input);
     1226
     1227                /* Start  looping on the number of gaussian points: */
     1228                Gauss* gauss=element->NewGauss(2);
     1229                for(int ig=gauss->begin();ig<gauss->end();ig++){
     1230
     1231                        gauss->GaussPoint(ig);
     1232
     1233                        /* Get Jacobian determinant: */
     1234                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     1235
     1236                        /*Get all parameters at gaussian point*/
     1237                        weights_input->GetInputValue(&weight,gauss,ThicknessAcrossGradientEnum);
     1238                        thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
     1239
     1240                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
     1241                        J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
     1242                }
     1243
     1244                /*clean up and Return: */
     1245                xDelete<IssmDouble>(xyz_list);
     1246                delete gauss;
    11651247        }
    11661248
Note: See TracChangeset for help on using the changeset viewer.