Changeset 17976 for issm/trunk-jpl/src/c/classes/FemModel.cpp
- Timestamp:
- 05/09/14 15:37:18 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/FemModel.cpp
r17907 r17976 1134 1134 void FemModel::BalancethicknessMisfitx(IssmDouble* presponse){/*{{{*/ 1135 1135 1136 IssmDouble J = 0; 1136 /*output: */ 1137 IssmDouble J=0.; 1137 1138 IssmDouble J_sum; 1138 1139 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:*/ 1143 1194 ISSM_MPI_Reduce (&J,&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); 1144 1195 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); … … 1151 1202 void FemModel::ThicknessAbsGradientx( IssmDouble* pJ){/*{{{*/ 1152 1203 1153 /*Intermediary*/1154 int i;1155 Element* element=NULL;1156 1157 1204 /*output: */ 1158 IssmDouble J=0 ;1205 IssmDouble J=0.; 1159 1206 IssmDouble J_sum; 1160 1207 1208 IssmDouble thickness,thicknessobs,weight; 1209 IssmDouble Jdet; 1210 IssmDouble* xyz_list = NULL; 1211 IssmDouble dp[3]; 1212 1161 1213 /*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; 1165 1247 } 1166 1248
Note:
See TracChangeset
for help on using the changeset viewer.