Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18833) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18834) @@ -337,5 +337,5 @@ /*Intermediaries */ - int num_responses,i,dim; + int num_responses,i,domaintype; IssmDouble Jdet,obs_velocity_mag,velocity_mag; IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; @@ -344,6 +344,6 @@ IssmDouble *xyz_list_top = NULL; - /*Get problem dimension*/ - element->FindParam(&dim,DomainDimensionEnum); + /* Get domaintype*/ + element->FindParam(&domaintype,DomainTypeEnum); /*Fetch number of nodes and dof for this finite element*/ @@ -353,5 +353,5 @@ /*Prepare coordinate system list*/ int* cs_list = xNew(vnumnodes+pnumnodes); - if(dim==2) for(i=0;iGetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); Input* vxobs_input = element->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + Input* vy_input=NULL; + Input* vyobs_input=NULL; + if (domaintype!=Domain2DverticalEnum) { + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + } IssmDouble epsvel = 2.220446049250313e-16; IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ @@ -389,7 +393,9 @@ vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); + if (domaintype!=Domain2DverticalEnum) { + vy_input->GetInputValue(&vy,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + } /*Loop over all requested responses*/ @@ -409,9 +415,13 @@ */ for(i=0;ivalues[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; - if(dim==3){ + if(domaintype!=Domain2DverticalEnum){ + dux=vxobs-vx; + pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; duy=vyobs-vy; - pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + } + else{ + dux=vxobs-vx; + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; } } @@ -430,10 +440,18 @@ */ for(i=0;ivalues[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; - pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + if(domaintype!=Domain2DverticalEnum){ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; + dux=scalex*(vxobs-vx); + duy=scaley*(vyobs-vy); + pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + } + else{ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + dux=scalex*(vxobs-vx); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + } } break; @@ -451,5 +469,5 @@ */ for(i=0;ivalues[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; - pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; } else{ @@ -465,5 +483,5 @@ scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag); dux=scale*vx; - pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; } } @@ -480,9 +498,16 @@ */ for(i=0;ivalues[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; - pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + if (domaintype!=Domain2DverticalEnum){ + scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); + dux=scale*(vxobs-vx); + duy=scale*(vyobs-vy); + pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + } + else{ + scale=1./(S*2*fabs(vx-vxobs)+epsvel); + dux=scale*(vxobs-vx); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + } } break; @@ -498,8 +523,14 @@ */ for(i=0;ivalues[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; - pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + if(domaintype!=Domain2DverticalEnum){ + dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); + duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); + pe->values[i*3+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + pe->values[i*3+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; + } + else{ + dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; + } } break; @@ -542,5 +573,5 @@ /*Intermediaries */ - int num_responses,i; + int num_responses,i,domaintype; IssmDouble Jdet,obs_velocity_mag,velocity_mag; IssmDouble vx,vy,vxobs,vyobs,dux,duy,weight; @@ -548,4 +579,7 @@ int *responses = NULL; IssmDouble *xyz_list_top = NULL; + + /* Get domaintype*/ + element->FindParam(&domaintype,DomainTypeEnum); /*Fetch number of nodes and dof for this finite element*/ @@ -562,7 +596,11 @@ Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); Input* vxobs_input = element->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + Input* vy_input=NULL; + Input* vyobs_input=NULL; + if(domaintype!=Domain2DverticalEnum){ + vy_input = element->GetInput(VyEnum); _assert_(vy_input); + vyobs_input = element->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + } IssmDouble epsvel = 2.220446049250313e-16; IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ @@ -584,8 +622,9 @@ vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - + if(domaintype!=Domain2DverticalEnum){ + vy_input->GetInputValue(&vy,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + } /*Loop over all requested responses*/ for(int resp=0;respvalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + dux=vxobs-vx; + duy=vyobs-vy; + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + dux=vxobs-vx; + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -623,10 +668,17 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; + dux=scalex*(vxobs-vx); + duy=scaley*(vyobs-vy); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + dux=scalex*(vxobs-vx); + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -644,11 +696,20 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; + obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; + scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); + dux=scale*vx; + duy=scale*vy; + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + velocity_mag =fabs(vx)+epsvel; + obs_velocity_mag=fabs(vxobs)+epsvel; + scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); + dux=scale*vx; + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -664,9 +725,16 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); + dux=scale*(vxobs-vx); + duy=scale*(vyobs-vy); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + scale=1./(S*2*fabs(vx-vxobs)+epsvel); + dux=scale*(vxobs-vx); + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -682,8 +750,14 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); + duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -701,4 +775,7 @@ break; case RheologyBbarAbsGradientEnum: + /*Nothing in P vector*/ + break; + case RheologyBAbsGradientEnum: /*Nothing in P vector*/ break; @@ -710,5 +787,5 @@ /*Transform coordinate system*/ - element->TransformLoadVectorCoord(pe,XYEnum); + if(domaintype!=Domain2DverticalEnum) element->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ @@ -731,4 +808,8 @@ case Domain2DhorizontalEnum: basalelement = element; + break; + case Domain2DverticalEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(); break; case Domain3DEnum: @@ -760,7 +841,11 @@ Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); Input* vxobs_input = basalelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); - Input* vyobs_input = basalelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + Input* vy_input=NULL; + Input* vyobs_input=NULL; + if(domaintype!=Domain2DverticalEnum){ + vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + vyobs_input = basalelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + } IssmDouble epsvel = 2.220446049250313e-16; IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ @@ -782,8 +867,9 @@ vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); vxobs_input->GetInputValue(&vxobs,gauss); - vyobs_input->GetInputValue(&vyobs,gauss); - + if(domaintype!=Domain2DverticalEnum){ + vy_input->GetInputValue(&vy,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + } /*Loop over all requested responses*/ for(int resp=0;respvalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + dux=vxobs-vx; + duy=vyobs-vy; + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else { + dux=vxobs-vx; + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -821,10 +913,17 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0; + dux=scalex*(vxobs-vx); + duy=scaley*(vyobs-vy); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0; + dux=scalex*(vxobs-vx); + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -842,11 +941,20 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + velocity_mag =sqrt(pow(vx, 2)+pow(vy, 2))+epsvel; + obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel; + scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); + dux=scale*vx; + duy=scale*vy; + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + velocity_mag =fabs(vx)+epsvel; + obs_velocity_mag=fabs(vxobs)+epsvel; + scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); + dux=scale*vx; + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -862,9 +970,16 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel); + dux=scale*(vxobs-vx); + duy=scale*(vyobs-vy); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + scale=1./(S*2*fabs(vx-vxobs)+epsvel); + dux=scale*(vxobs-vx); + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -880,8 +995,14 @@ */ for(i=0;ivalues[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; - pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + if(domaintype!=Domain2DverticalEnum){ + dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); + duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel); + pe->values[i*2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel); + pe->values[i]+=dux*weight*Jdet*gauss->weight*basis[i]; + } } break; @@ -899,4 +1020,7 @@ break; case RheologyBbarAbsGradientEnum: + /*Nothing in P vector*/ + break; + case RheologyBAbsGradientEnum: /*Nothing in P vector*/ break; @@ -908,5 +1032,5 @@ /*Transform coordinate system*/ - basalelement->TransformLoadVectorCoord(pe,XYEnum); + if(domaintype!=Domain2DverticalEnum) basalelement->TransformLoadVectorCoord(pe,XYEnum); /*Clean up and return*/ @@ -951,5 +1075,6 @@ if(control_type!=MaterialsRheologyBbarEnum && control_type!=FrictionCoefficientEnum && - control_type!=DamageDbarEnum){ + control_type!=DamageDbarEnum && + control_type!=MaterialsRheologyBEnum){ _error_("Control "<* gradient,int control_index){/*{{{*/ +void AdjointHorizAnalysis::GradientJBbarGradient(Element* element,Vector* gradient,int control_index){/*{{{*/ /*Intermediaries*/ - int domaintype,dim; + int domaintype; Element* basalelement; @@ -1096,15 +1231,12 @@ case Domain2DhorizontalEnum: basalelement = element; - dim = 2; break; case Domain2DverticalEnum: if(!element->IsOnBase()) return; basalelement = element->SpawnBasalElement(); - dim = 1; break; case Domain3DEnum: if(!element->IsOnBase()) return; basalelement = element->SpawnBasalElement(); - dim = 2; break; default: _error_("mesh "<weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]); } @@ -1162,4 +1294,70 @@ delete gauss; if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; +}/*}}}*/ +void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + int domaintype; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + break; + case Domain2DverticalEnum: + break; + case Domain3DEnum: + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* dbasis = xNew(3*numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinates(&xyz_list); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* rheology_input = element->GetInput(MaterialsRheologyBEnum); _assert_(rheology_input); + Input* weights_input = element->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss); + weights_input->GetInputValue(&weight,gauss,RheologyBAbsGradientEnum); + + /*Build alpha_complement_list: */ + rheology_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for(int i=0;iweight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]); + } + else{ + ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]+dbasis[2*numvertices+i]*dk[2]); + } + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(dbasis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; }/*}}}*/ void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ @@ -1263,5 +1461,4 @@ /*Intermediaries*/ - int dim=3; IssmDouble Jdet,weight; IssmDouble drag,dalpha2dk; @@ -1269,4 +1466,6 @@ IssmDouble *xyz_list_base= NULL; + int domaintype,dim; + element->FindParam(&domaintype,DomainTypeEnum); /*Fetch number of vertices for this finite element*/ int numvertices = element->GetNumberOfVertices(); @@ -1278,4 +1477,82 @@ /*Build friction element, needed later: */ + if(domaintype!=Domain2DverticalEnum) dim=3; + else dim=2; + Friction* friction=new Friction(element,dim); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = NULL; + Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = NULL; + Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); + if(domaintype!=Domain2DverticalEnum){ + vy_input = element->GetInput(VyEnum); _assert_(vy_input); + adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); + } + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + adjointx_input->GetInputValue(&lambda, gauss); + vx_input->GetInputValue(&vx,gauss); + if(domaintype!=Domain2DverticalEnum){ + adjointy_input->GetInputValue(&mu, gauss); + vy_input->GetInputValue(&vy,gauss); + } + dragcoeff_input->GetInputValue(&drag, gauss); + + friction->GetAlphaComplement(&dalpha2dk,gauss); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dD): */ + for(int i=0;iweight*basis[i]; + else ge[i]+=-2.*drag*dalpha2dk*(lambda*vx)*Jdet*gauss->weight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list_base); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + delete friction; +}/*}}}*/ +void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*return if floating or not on bed (gradient is 0)*/ + if(element->IsFloating()) return; + if(!element->IsOnBase()) return; + + /*Intermediaries*/ + int domaintype,dim; + IssmDouble Jdet,weight; + IssmDouble drag,dalpha2dk,normal[3]; + IssmDouble vx,vy,vz,lambda,mu,xi; + IssmDouble *xyz_list_base= NULL; + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /* get domaintype */ + element->FindParam(&domaintype,DomainTypeEnum); + + /*Build friction element, needed later: */ + if(domaintype!=Domain2DverticalEnum) dim=3; + else dim=2; Friction* friction=new Friction(element,dim); @@ -1287,4 +1564,10 @@ Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* vz_input = NULL; + Input* adjointz_input = NULL; + if(domaintype!=Domain2DverticalEnum){ + Input* vz_input = element->GetInput(VzEnum); _assert_(vy_input); + Input* adjointz_input = element->GetInput(AdjointzEnum); _assert_(adjointz_input); + } Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); @@ -1298,89 +1581,35 @@ vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); + if(domaintype!=Domain2DverticalEnum){ + adjointz_input->GetInputValue(&xi ,gauss); + vz_input->GetInputValue(&vz,gauss); + } dragcoeff_input->GetInputValue(&drag, gauss); friction->GetAlphaComplement(&dalpha2dk,gauss); + element->NormalBase(&normal[0],xyz_list_base); element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); element->NodalFunctionsP1(basis,gauss); - /*Build gradient vector (actually -dJ/dD): */ - for(int i=0;iweight*basis[i]; - _assert_(!xIsNan(ge[i])); - } - } - gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); - - /*Clean up and return*/ - xDelete(xyz_list_base); - xDelete(basis); - xDelete(ge); - xDelete(vertexpidlist); - delete gauss; - delete friction; -}/*}}}*/ -void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector* gradient,int control_index){/*{{{*/ - - /*return if floating or not on bed (gradient is 0)*/ - if(element->IsFloating()) return; - if(!element->IsOnBase()) return; - - /*Intermediaries*/ - int dim=3; - IssmDouble Jdet,weight; - IssmDouble drag,dalpha2dk,normal[3]; - IssmDouble vx,vy,vz,lambda,mu,xi; - IssmDouble *xyz_list_base= NULL; - - /*Fetch number of vertices for this finite element*/ - int numvertices = element->GetNumberOfVertices(); - - /*Initialize some vectors*/ - IssmDouble* basis = xNew(numvertices); - IssmDouble* ge = xNewZeroInit(numvertices); - int* vertexpidlist = xNew(numvertices); - - /*Build friction element, needed later: */ - Friction* friction=new Friction(element,dim); - - /*Retrieve all inputs we will be needing: */ - element->GetVerticesCoordinatesBase(&xyz_list_base); - element->GradientIndexing(&vertexpidlist[0],control_index); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = element->GetInput(VzEnum); _assert_(vy_input); - Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); - Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); - Input* adjointz_input = element->GetInput(AdjointzEnum); _assert_(adjointz_input); - Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGaussBase(4); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); - - adjointx_input->GetInputValue(&lambda, gauss); - adjointy_input->GetInputValue(&mu, gauss); - adjointz_input->GetInputValue(&xi ,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); - dragcoeff_input->GetInputValue(&drag, gauss); - - friction->GetAlphaComplement(&dalpha2dk,gauss); - element->NormalBase(&normal[0],xyz_list_base); - - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctionsP1(basis,gauss); - /*Build gradient vector (actually -dJ/dk): */ - for(int i=0;iweight*basis[i]; - _assert_(!xIsNan(ge[i])); + if(domaintype!=Domain2DverticalEnum){ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + else{ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } } } @@ -1493,4 +1722,168 @@ this->GradientJBbarSSA(element,gradient,control_index); }/*}}}*/ +void AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ + + /*Intermediaries*/ + int domaintype,dim; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + dim = 2; + break; + case Domain2DverticalEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 1; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + dim = 2; + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->GradientIndexing(&vertexpidlist[0],control_index); + Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + thickness_input->GetInputValue(&thickness,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss); + adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss); + + basalelement->dViscositydBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dB): */ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; +}/*}}}*/ +void AdjointHorizAnalysis::GradientJBHO(Element* element,Vector* gradient,int control_index){/*{{{*/ + /*Intermediaries*/ + int domaintype,dim; + + /*Get domaintype*/ + element->FindParam(&domaintype,DomainTypeEnum); + + /*Intermediaries*/ + IssmDouble Jdet,weight; + IssmDouble thickness,dmudB; + IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3]; + IssmDouble *xyz_list= NULL; + + /*Fetch number of vertices for this finite element*/ + int numvertices = element->GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + + /*Retrieve all inputs we will be needing: */ + element->GetVerticesCoordinates(&xyz_list); + element->GradientIndexing(&vertexpidlist[0],control_index); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = NULL; + Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input = NULL; + Input* rheologyb_input = element->GetInput(MaterialsRheologyBEnum); _assert_(rheologyb_input); + if(domaintype!=Domain2DverticalEnum){ + vy_input = element->GetInput(VyEnum); _assert_(vy_input); + adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input); + } + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(4); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + thickness_input->GetInputValue(&thickness,gauss); + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss); + dim=2; + if(domaintype!=Domain2DverticalEnum){ + adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list, gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + dim=3; + } + + element->dViscositydBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dB): */ + for(int i=0;iweight*basis[i]; + } + else{ + ge[i]+=-dmudB*thickness*4*dvx[0]*dadjx[0]*Jdet*gauss->weight*basis[i]; + } + _assert_(!xIsNan(ge[i])); + } + } + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; +}/*}}}*/ +void AdjointHorizAnalysis::GradientJBFS(Element* element,Vector* gradient,int control_index){/*{{{*/ + /*WARNING: We use HO as an estimate for now*/ + this->GradientJBHO(element,gradient,control_index); +}/*}}}*/ void AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector* gradient,int control_index){/*{{{*/ @@ -1590,8 +1983,12 @@ int* doflist=NULL; + int domaintype; + element->FindParam(&domaintype,DomainTypeEnum); + /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); - int numdof = numnodes*2; - + int numdof; + if(domaintype!=Domain2DverticalEnum) numdof = numnodes*2; + else numdof = numnodes*1; /*Fetch dof list and allocate solution vectors*/ element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum); @@ -1604,14 +2001,16 @@ /*Transform solution in Cartesian Space*/ - element->TransformSolutionCoord(&values[0],XYEnum); + if(domaintype!=Domain2DverticalEnum) element->TransformSolutionCoord(&values[0],XYEnum); /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ for(i=0;i(lambdax[i])) _error_("NaN found in solution vector"); - if(xIsNan(lambday[i])) _error_("NaN found in solution vector"); + if(domaintype!=Domain2DverticalEnum && xIsNan(lambday[i])) _error_("NaN found in solution vector"); } @@ -1627,10 +2026,24 @@ }/*}}}*/ void AdjointHorizAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/ - int i,dim; + int i,fe_FS; int* vdoflist=NULL; int* pdoflist=NULL; IssmDouble FSreconditioning; - element->FindParam(&dim,DomainDimensionEnum); + int domaintype,dim; + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + dim = 3; + break; + case Domain2DverticalEnum: + dim = 2; + break; + case Domain3DEnum: + dim = 3; + break; + default: _error_("mesh "<(lambdax[i])) _error_("NaN found in solution vector"); - lambday[i] = values[i*NDOF3+1]; if(xIsNan(lambday[i])) _error_("NaN found in solution vector"); - lambdaz[i] = values[i*NDOF3+2]; if(xIsNan(lambdaz[i])) _error_("NaN found in solution vector"); + lambdax[i] = values[i*dim+0]; if(xIsNan(lambdax[i])) _error_("NaN found in solution vector"); + lambday[i] = values[i*dim+1]; if(xIsNan(lambday[i])) _error_("NaN found in solution vector"); + if(dim==3){ + lambdaz[i] = values[i*dim+2]; if(xIsNan(lambdaz[i])) _error_("NaN found in solution vector"); + } } for(i=0;iAddInput(AdjointxEnum,lambdax,element->VelocityInterpolation()); element->AddInput(AdjointyEnum,lambday,element->VelocityInterpolation()); - element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation()); - element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation()); + if(domaintype!=Domain2DverticalEnum) element->AddInput(AdjointzEnum,lambdaz,element->VelocityInterpolation()); + + element->FindParam(&fe_FS,FlowequationFeFSEnum); + if(fe_FS!=LATaylorHoodEnum && fe_FS!=LACrouzeixRaviartEnum) + element->AddInput(AdjointpEnum,lambdap,element->PressureInterpolation()); /*Free ressources:*/