![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 15718
- Timestamp:
- 08/06/13 08:54:01 (12 years ago)
- Location:
- issm/trunk-jpl/src/c/classes/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r15717 r15718 5187 5187 } 5188 5188 /*}}}*/ 5189 /*FUNCTION Tria::CreatePVectorAdjointFS{{{*/5190 ElementVector* Tria::CreatePVectorAdjointFS(void){5191 5192 /*Intermediaries */5193 int i,resp;5194 int *responses=NULL;5195 int num_responses;5196 IssmDouble Jdet;5197 IssmDouble obs_velocity_mag,velocity_mag;5198 IssmDouble dux,duy;5199 IssmDouble epsvel=2.220446049250313e-16;5200 IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/5201 IssmDouble scalex=0,scaley=0,scale=0,S=0;5202 IssmDouble vx,vy,vxobs,vyobs,weight;5203 IssmDouble xyz_list[NUMVERTICES][3];5204 IssmDouble basis[3];5205 GaussTria* gauss=NULL;5206 5207 /*Initialize Element vector*/5208 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);5209 5210 /*Retrieve all inputs and parameters*/5211 GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);5212 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);5213 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);5214 Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);5215 Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);5216 Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);5217 Input* vxobs_input = inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);5218 Input* vyobs_input = inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);5219 5220 /*Get Surface if required by one response*/5221 for(resp=0;resp<num_responses;resp++){5222 if(responses[resp]==SurfaceAverageVelMisfitEnum){5223 inputs->GetInputValue(&S,SurfaceAreaEnum); break;5224 }5225 }5226 5227 /* Start looping on the number of gaussian points: */5228 gauss=new GaussTria(4);5229 for(int ig=gauss->begin();ig<gauss->end();ig++){5230 5231 gauss->GaussPoint(ig);5232 5233 /* Get Jacobian determinant: */5234 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);5235 5236 /*Get all parameters at gaussian point*/5237 vx_input->GetInputValue(&vx,gauss);5238 vy_input->GetInputValue(&vy,gauss);5239 vxobs_input->GetInputValue(&vxobs,gauss);5240 vyobs_input->GetInputValue(&vyobs,gauss);5241 GetNodalFunctions(basis, gauss);5242 5243 /*Loop over all requested responses*/5244 for(resp=0;resp<num_responses;resp++){5245 5246 weights_input->GetInputValue(&weight,gauss,resp);5247 5248 switch(responses[resp]){5249 5250 case SurfaceAbsVelMisfitEnum:5251 /*5252 * 1 [ 2 2 ]5253 * J = --- | (u - u ) + (v - v ) |5254 * 2 [ obs obs ]5255 *5256 * dJ5257 * DU = - -- = (u - u )5258 * du obs5259 */5260 for (i=0;i<NUMVERTICES;i++){5261 dux=vxobs-vx;5262 duy=vyobs-vy;5263 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];5264 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];5265 }5266 break;5267 case SurfaceRelVelMisfitEnum:5268 /*5269 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]5270 * J = --- | ------------- (u - u ) + ------------- (v - v ) |5271 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]5272 * obs obs5273 *5274 * dJ \bar{v}^25275 * DU = - -- = ------------- (u - u )5276 * du (u + eps)^2 obs5277 * obs5278 */5279 for (i=0;i<NUMVERTICES;i++){5280 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;5281 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;5282 dux=scalex*(vxobs-vx);5283 duy=scaley*(vyobs-vy);5284 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];5285 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];5286 }5287 break;5288 case SurfaceLogVelMisfitEnum:5289 /*5290 * [ vel + eps ] 25291 * J = 4 \bar{v}^2 | log ( ----------- ) |5292 * [ vel + eps ]5293 * obs5294 *5295 * dJ 2 * log(...)5296 * DU = - -- = - 4 \bar{v}^2 ------------- u5297 * du vel^2 + eps5298 *5299 */5300 for (i=0;i<NUMVERTICES;i++){5301 velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel;5302 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;5303 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);5304 dux=scale*vx;5305 duy=scale*vy;5306 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];5307 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];5308 }5309 break;5310 case SurfaceAverageVelMisfitEnum:5311 /*5312 * 1 2 25313 * J = --- sqrt( (u - u ) + (v - v ) )5314 * S obs obs5315 *5316 * dJ 1 15317 * DU = - -- = - --- ----------- * 2 (u - u )5318 * du S 2 sqrt(...) obs5319 */5320 for (i=0;i<NUMVERTICES;i++){5321 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);5322 dux=scale*(vxobs-vx);5323 duy=scale*(vyobs-vy);5324 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];5325 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];5326 }5327 break;5328 case SurfaceLogVxVyMisfitEnum:5329 /*5330 * 1 [ |u| + eps 2 |v| + eps 2 ]5331 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |5332 * 2 [ |u |+ eps |v |+ eps ]5333 * obs obs5334 * dJ 1 u 15335 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------5336 * du |u| + eps |u| u + eps5337 */5338 for (i=0;i<NUMVERTICES;i++){5339 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);5340 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);5341 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];5342 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];5343 }5344 break;5345 case DragCoefficientAbsGradientEnum:5346 /*Nothing in P vector*/5347 break;5348 case ThicknessAbsGradientEnum:5349 /*Nothing in P vector*/5350 break;5351 case ThicknessAcrossGradientEnum:5352 /*Nothing in P vector*/5353 break;5354 case ThicknessAlongGradientEnum:5355 /*Nothing in P vector*/5356 break;5357 case RheologyBbarAbsGradientEnum:5358 /*Nothing in P vector*/5359 break;5360 default:5361 _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");5362 }5363 }5364 }5365 5366 /*Clean up and return*/5367 delete gauss;5368 xDelete<int>(responses);5369 return pe;5370 }5371 /*}}}*/5372 5189 /*FUNCTION Tria::DragCoefficientAbsGradient{{{*/ 5373 5190 IssmDouble Tria::DragCoefficientAbsGradient(int weight_index){ … … 6050 5867 ElementVector* Tria::CreatePVectorHydrologyDCInefficient(void){ 6051 5868 6052 /*Constants*/6053 const int numdof=NDOF1*NUMVERTICES;6054 6055 5869 /*Intermediaries */ 6056 5870 IssmDouble Jdet; … … 6059 5873 IssmDouble water_load,transfer; 6060 5874 IssmDouble sediment_storing; 6061 IssmDouble basis[numdof]; 6062 GaussTria* gauss=NULL; 6063 6064 /*Initialize Element vector*/ 6065 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); 5875 5876 /*Fetch number of nodes and dof for this finite element*/ 5877 int numnodes = this->NumberofNodes(); 5878 5879 /*Initialize Element vector and other vectors*/ 5880 ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); 5881 IssmDouble* basis = xNew<IssmDouble>(numnodes); 6066 5882 6067 5883 /*Retrieve all inputs and parameters*/ … … 6078 5894 6079 5895 /* Start looping on the number of gaussian points: */ 6080 gauss=new GaussTria(2);5896 GaussTria* gauss=new GaussTria(2); 6081 5897 for(int ig=gauss->begin();ig<gauss->end();ig++){ 6082 5898 … … 6090 5906 scalar = Jdet*gauss->weight*(water_load+transfer); 6091 5907 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt; 6092 for(int i=0;i<num dof;i++) pe->values[i]+=scalar*basis[i];5908 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i]; 6093 5909 6094 5910 /*Transient term*/ … … 6096 5912 old_wh_input->GetInputValue(&water_head,gauss); 6097 5913 scalar = Jdet*gauss->weight*water_head*sediment_storing; 6098 for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i]; 6099 } 6100 } 5914 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i]; 5915 } 5916 } 5917 6101 5918 /*Clean up and return*/ 5919 xDelete<IssmDouble>(basis); 6102 5920 delete gauss; 6103 5921 return pe; … … 6106 5924 /*FUNCTION Tria::CreatePVectorHydrologyDCEfficient {{{*/ 6107 5925 ElementVector* Tria::CreatePVectorHydrologyDCEfficient(void){ 6108 6109 /*Constants*/6110 const int numdof=NDOF1*NUMVERTICES;6111 5926 6112 5927 /*Intermediaries */ … … 6117 5932 IssmDouble transfer,residual; 6118 5933 IssmDouble epl_storing; 6119 IssmDouble basis[numdof]; 6120 GaussTria* gauss=NULL; 5934 GaussTria* gauss = NULL; 6121 5935 6122 5936 /*Check that all nodes are active, else return empty matrix*/ … … 6124 5938 return NULL; 6125 5939 } 6126 /*Initialize Element vector*/ 6127 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); 5940 5941 /*Fetch number of nodes and dof for this finite element*/ 5942 int numnodes = this->NumberofNodes(); 5943 5944 /*Initialize Element vector and other vectors*/ 5945 ElementVector* pe = new ElementVector(nodes,numnodes,this->parameters); 5946 IssmDouble* basis = xNew<IssmDouble>(numnodes); 6128 5947 6129 5948 /*Retrieve all inputs and parameters*/ … … 6145 5964 gauss->GaussPoint(ig); 6146 5965 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 6147 GetNodalFunctions(basis, 5966 GetNodalFunctions(basis,gauss); 6148 5967 6149 5968 /*Loading term*/ … … 6151 5970 scalar = Jdet*gauss->weight*(-transfer); 6152 5971 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt; 6153 for(int i=0;i<num dof;i++) pe->values[i]+=scalar*basis[i];5972 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i]; 6154 5973 6155 5974 /*Transient term*/ … … 6157 5976 old_wh_input->GetInputValue(&water_head,gauss); 6158 5977 scalar = Jdet*gauss->weight*water_head*epl_storing; 6159 for(int i=0;i<num dof;i++) pe->values[i]+=scalar*basis[i];5978 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i]; 6160 5979 } 6161 5980 } … … 6171 5990 pe->values[iv]+=residual/connectivity; 6172 5991 } 5992 6173 5993 /*Clean up and return*/ 5994 xDelete<IssmDouble>(basis); 6174 5995 delete gauss; 6175 5996 return pe; -
issm/trunk-jpl/src/c/classes/Elements/Tria.h
r15695 r15718 241 241 ElementMatrix* CreateKMatrixAdjointSSA(void); 242 242 ElementVector* CreatePVectorAdjointHoriz(void); 243 ElementVector* CreatePVectorAdjointFS(void);244 243 ElementVector* CreatePVectorAdjointBalancethickness(void); 245 244 void InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution);
Note:
See TracChangeset
for help on using the changeset viewer.
