Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 26477) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 26478) @@ -52,4 +52,8 @@ case FSApproximationEnum: return CreateKMatrixFS(element); + case MLHOApproximationEnum: + // a more accurate option, but integrate in the vertical direction numerically. + // return CreateKMatrixMLHOVerticalIntergrated(element); + return CreateKMatrixMLHO(element); case NoneApproximationEnum: return NULL; @@ -215,4 +219,277 @@ xDelete(xyz_list); return Ke; +}/*}}}*/ +ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/ + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*Intermediaries */ + bool incomplete_adjoint; + IssmDouble Jdet,mu_prime,n,thickness,mu,effmu; + IssmDouble *xyz_list = NULL; + IssmDouble viscosity[9]; //9 mu for different integrand + int domaintype; + int dim=2; + + IssmDouble eb1i,eb1j,esh1i,esh1j,eb2i,eb2j,esh2i,esh2j; + IssmDouble epsilon[5],epsilonbase[5],epsilonshear[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble e1b[2], e2b[2], e1sh[2], e2sh[2]; + IssmDouble vxshear, vyshear; + + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: case Domain2DverticalEnum: + _error_("mesh "<GetNumberOfNodes(); + IssmDouble* dbasis = xNew(2*numnodes); // like SSA + IssmDouble* basis = xNew(numnodes); // like SSA + + /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ + element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); + ElementMatrix* Ke=analysis->CreateKMatrix(basalelement); + delete analysis; + if(incomplete_adjoint) return Ke; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); //vertically integrated vx + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); //vertically integrated vy + Input* vxbase_input = element->GetInput(VxBaseEnum); _assert_(vxbase_input); + Input* vybase_input = element->GetInput(VyBaseEnum); _assert_(vybase_input); + Input* vxshear_input = element->GetInput(VxShearEnum); _assert_(vxshear_input); + Input* vyshear_input = element->GetInput(VyShearEnum); _assert_(vyshear_input); + Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(5); + Gauss* gauss_base = basalelement->NewGauss(); + while(gauss->next()){ + gauss->SynchronizeGaussBase(gauss_base); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); + basalelement->NodalFunctions(basis, gauss_base); + + thickness_input->GetInputValue(&thickness, gauss); + n_input->GetInputValue(&n,gauss); + + vxshear_input->GetInputValue(&vxshear,gauss); + vyshear_input->GetInputValue(&vyshear,gauss); + + element->material->ViscosityMLHOAdjoint(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input); + + effmu = 2.0*(1-n)/2.0/n; + + element->StrainRateHO(&epsilonbase[0],xyz_list,gauss,vxbase_input, vybase_input); + element->StrainRateHO(&epsilonshear[0],xyz_list,gauss,vxshear_input, vyshear_input); + + e1b[0] = 2*epsilonbase[0]+epsilonbase[1]; e1b[1] = epsilonbase[2]; + e2b[1] = epsilonbase[0]+2*epsilonbase[1]; e2b[0] = epsilonbase[2]; + + e1sh[0] = 2*epsilonshear[0]+epsilonshear[1]; e1sh[1] = epsilonshear[2]; + e2sh[1] = epsilonshear[0]+2*epsilonshear[1]; e2sh[0] = epsilonshear[2]; + + for(int i=0;ivalues[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb1i+viscosity[1]*(esh1j*eb1i+eb1j*esh1i)+viscosity[2]*esh1j*esh1i); + Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1j*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb1j*eb2i+viscosity[1]*(esh1j*eb2i+eb1j*esh2i)+viscosity[2]*esh1j*esh2i); + Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1j*basis[i]); + + Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb1i+viscosity[2]*(esh1j*eb1i+eb1j*esh1i)+viscosity[4]*esh1j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh1i*basis[j]); + Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb1i+viscosity[4]*(esh1j*eb1i+eb1j*esh1i)+viscosity[5]*esh1j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(eb1i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*(esh1i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb1j*eb2i+viscosity[2]*(esh1j*eb2i+eb1j*esh2i)+viscosity[4]*esh1j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2i*basis[j]); + Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb1j*eb2i+viscosity[4]*(esh1j*eb2i+eb1j*esh2i)+viscosity[5]*esh1j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*eb2i*basis[j]+vyshear*eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2i*basis[j]+vyshear*esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vyshear*basis[j]*basis[i]); + + + Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb1i+viscosity[1]*(esh2j*eb1i+eb2j*esh1i)+viscosity[2]*esh2j*esh1i); + Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vxshear*esh2j*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[0]*eb2j*eb2i+viscosity[1]*(esh2j*eb2i+eb2j*esh2i)+viscosity[2]*esh2j*esh2i); + Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2j*basis[i]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2j*basis[i]); + + + Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb1i+viscosity[2]*(esh2j*eb1i+eb2j*esh1i)+viscosity[4]*esh2j*esh1i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb1i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh1i*basis[j]); + Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb1i+viscosity[4]*(esh2j*eb1i+eb2j*esh1i)+viscosity[5]*esh2j*esh1i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*(vyshear*eb1i*basis[j]+vxshear*eb2j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*(vxshear*esh2j*basis[i]+vyshear*esh1i*basis[j])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vxshear*vxshear*basis[j]*basis[i]); + Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[1]*eb2j*eb2i+viscosity[2]*(esh2j*eb2i+eb2j*esh2i)+viscosity[4]*esh2j*esh2i+viscosity[3]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*eb2i*basis[j]+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*esh2i*basis[j]); + Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*(viscosity[2]*eb2j*eb2i+viscosity[4]*(esh2j*eb2i+eb2j*esh2i)+viscosity[5]*esh2j*esh2i+viscosity[6]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(eb2i*basis[j]+eb1j*basis[i])+viscosity[7]*(n+1)*(n+1)/2.0/thickness/thickness*vyshear*(esh2i*basis[j]+esh1j*basis[i])+viscosity[8]*(n+1)*(n+1)*(n+1)*(n+1)/4.0/thickness/thickness/thickness/thickness*vyshear*vyshear*basis[j]*basis[i]); + } + } + } + + /*Transform Coordinate System*/ + /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/ + + /*Clean up and return*/ + delete gauss; + delete gauss_base; + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; + xDelete(xyz_list); + xDelete(dbasis); + xDelete(basis); + return Ke; +}/*}}}*/ +ElementMatrix* AdjointHorizAnalysis::CreateKMatrixMLHOVerticalIntergrated(Element* element){/*{{{*/ + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*Intermediaries */ + bool incomplete_adjoint; + IssmDouble Jdet,mu_prime,n,thickness,mu,effmu; + IssmDouble *xyz_list = NULL; + IssmDouble zeta, epsilon_eff; + int domaintype; + int dim=2; + + IssmDouble e1phi1i, e1phi1j, e2phi1i, e2phi1j, e1phi2i, e1phi2j, e2phi2i, e2phi2j; + IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble e1[3],e2[3]; + + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: case Domain2DverticalEnum: + _error_("mesh "<GetNumberOfNodes(); + IssmDouble* dbasis = xNew(2*numnodes); // like SSA + IssmDouble* basis = xNew(numnodes); // like SSA + + /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/ + element->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum); + StressbalanceAnalysis* analysis = new StressbalanceAnalysis(); + ElementMatrix* Ke=analysis->CreateKMatrix(basalelement); + delete analysis; + if(incomplete_adjoint) return Ke; + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vxbase_input = element->GetInput(VxBaseEnum); _assert_(vxbase_input); + Input* vybase_input = element->GetInput(VyBaseEnum); _assert_(vybase_input); + Input* vxshear_input = element->GetInput(VxShearEnum); _assert_(vxshear_input); + Input* vyshear_input = element->GetInput(VyShearEnum); _assert_(vyshear_input); + Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss = element->NewGauss(10); + Gauss* gauss_base = basalelement->NewGauss(); + + GaussSeg* gauss_seg=new GaussSeg(5); + + while(gauss->next()){ + gauss->SynchronizeGaussBase(gauss_base); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); + basalelement->NodalFunctions(basis, gauss_base); + + thickness_input->GetInputValue(&thickness, gauss); + n_input->GetInputValue(&n,gauss); + + effmu = 2*(1-n)/2.0/n; + + /* Get the integration in the vertical direction */ + gauss_seg->Reset(); + while(gauss_seg->next()){ + /*Compute zeta for gauss_seg point (0=surface, 1=base)*/ + zeta=0.5*(gauss_seg->coord1+1); + + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy (for a given zeta)*/ + element->StrainRateMLHO(&epsilon[0],xyz_list,gauss, + vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); + epsilon_eff=sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] + + epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]); + + /*Get viscosity at zeta*/ + element->material->GetViscosity(&mu,epsilon_eff,gauss); + /*the adjoint viscosity with zeta dependent term*/ + mu = mu /epsilon_eff/epsilon_eff; + + e1[0] = 2.0*epsilon[0]+epsilon[1]; + e1[1] = epsilon[2]; + e1[2] = epsilon[3]; + + e2[0] = epsilon[2]; + e2[1] = epsilon[0]+2.0*epsilon[1]; + e2[2] = epsilon[4]; + + for(int i=0;ivalues[4*numnodes*(4*i+0)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e1phi1j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+1)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e1phi1j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+2)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e1phi1j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+3)+4*j+0]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e1phi1j*thickness*0.5*gauss_seg->weight; + + Ke->values[4*numnodes*(4*i+0)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e1phi2j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+1)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e1phi2j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+2)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e1phi2j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+3)+4*j+1]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e1phi2j*thickness*0.5*gauss_seg->weight; + + Ke->values[4*numnodes*(4*i+0)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e2phi1j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+1)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e2phi1j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+2)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e2phi1j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+3)+4*j+2]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e2phi1j*thickness*0.5*gauss_seg->weight; + + Ke->values[4*numnodes*(4*i+0)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e1phi1i*e2phi2j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+1)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e1phi2i*e2phi2j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+2)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e2phi1i*e2phi2j*thickness*0.5*gauss_seg->weight; + Ke->values[4*numnodes*(4*i+3)+4*j+3]+=gauss->weight*Jdet*effmu*mu*e2phi2i*e2phi2j*thickness*0.5*gauss_seg->weight; + } + } + } + } + /*Transform Coordinate System*/ + /*element->TransformStiffnessMatrixCoord(Ke,XYEnum);*/ + + /*Clean up and return*/ + delete gauss; + delete gauss_base; + delete gauss_seg; + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; + xDelete(xyz_list); + xDelete(dbasis); + xDelete(basis); + return Ke; }/*}}}*/ ElementMatrix* AdjointHorizAnalysis::CreateKMatrixL1L2(Element* element){/*{{{*/ @@ -335,4 +612,6 @@ case FSApproximationEnum: return CreatePVectorFS(element); + case MLHOApproximationEnum: + return CreatePVectorMLHO(element); case NoneApproximationEnum: return NULL; @@ -820,4 +1099,259 @@ }/*}}}*/ +ElementVector* AdjointHorizAnalysis::CreatePVectorMLHO(Element* element){/*{{{*/ + + /*Intermediaries*/ + int domaintype; + Element* basalelement; + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain2DverticalEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + case Domain3DEnum: + if(!element->IsOnBase()) return NULL; + basalelement = element->SpawnBasalElement(true); + break; + default: _error_("mesh "<GetNumberOfNodes(); + + /*Initialize Element vector and vectors*/ + ElementVector* pe = basalelement->NewElementVector(MLHOApproximationEnum); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum); + basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum); + DatasetInput* weights_input = basalelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input = basalelement->GetInput(VxSurfaceEnum); _assert_(vx_input); + Input* vxobs_input = basalelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); + Input* thickness_input=basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* n_input =element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* vy_input=NULL; + Input* vyobs_input=NULL; + if(domaintype!=Domain2DverticalEnum){ + vy_input = basalelement->GetInput(VySurfaceEnum); _assert_(vy_input); + vyobs_input = basalelement->GetInput(InversionVyObsEnum); _assert_(vyobs_input); + } + IssmDouble epsvel = 2.220446049250313e-16; + IssmDouble meanvel = 3.170979198376458e-05; /*1000 m/yr*/ + + /*Get Surface if required by one response*/ + Input* S_input = NULL; + for(int resp=0;respGetInput(SurfaceAreaEnum); _assert_(S_input); break; + } + } + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(4); + while(gauss->next()){ + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctions(basis, gauss); + + thickness_input->GetInputValue(&thickness,gauss); + n_input->GetInputValue(&n,gauss); + vx_input->GetInputValue(&vx,gauss); + vxobs_input->GetInputValue(&vxobs,gauss); + if(domaintype!=Domain2DverticalEnum){ + vy_input->GetInputValue(&vy,gauss); + vyobs_input->GetInputValue(&vyobs,gauss); + } + /*Loop over all requested responses*/ + for(int resp=0;respGetInputValue(&weight,gauss,responses[resp]); + + switch(responses[resp]){ + case SurfaceAbsVelMisfitEnum: + /* + * 1 [ 2 2 ] + * J = --- | (u - u ) + (v - v ) | + * 2 [ obs obs ] + * + * dJ + * DU = - -- = (u - u ) + * du obs + */ + for(i=0;ivalues[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else { + _error_("2D vertical is not implemented for MLHO"); + } + } + break; + case SurfaceRelVelMisfitEnum: + /* + * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] + * J = --- | ------------- (u - u ) + ------------- (v - v ) | + * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] + * obs obs + * + * dJ \bar{v}^2 + * DU = - -- = ------------- (u - u ) + * du (u + eps)^2 obs + * obs + */ + for(i=0;ivalues[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + _error_("2D vertical is not implemented for MLHO"); + } + } + break; + case SurfaceLogVelMisfitEnum: + /* + * [ vel + eps ] 2 + * J = 4 \bar{v}^2 | log ( ----------- ) | + * [ vel + eps ] + * obs + * + * dJ 2 * log(...) + * DU = - -- = - 4 \bar{v}^2 ------------- u + * du vel^2 + eps + * + */ + for(i=0;ivalues[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + _error_("2D vertical is not implemented for MLHO"); + } + } + break; + case SurfaceAverageVelMisfitEnum: + /* + * 1 2 2 + * J = --- sqrt( (u - u ) + (v - v ) ) + * S obs obs + * + * dJ 1 1 + * DU = - -- = - --- ----------- * 2 (u - u ) + * du S 2 sqrt(...) obs + */ + S_input->GetInputValue(&S,gauss); + for(i=0;ivalues[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + _error_("2D vertical is not implemented for MLHO"); + } + } + break; + case SurfaceLogVxVyMisfitEnum: + /* + * 1 [ |u| + eps 2 |v| + eps 2 ] + * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | + * 2 [ |u |+ eps |v |+ eps ] + * obs obs + * dJ 1 u 1 + * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------ + * du |u| + eps |u| u + eps + */ + for(i=0;ivalues[i*4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+1]+=dux*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+2]+=duy*weight*Jdet*gauss->weight*basis[i]; + pe->values[i*4+3]+=duy*weight*Jdet*gauss->weight*basis[i]; + } + else{ + _error_("2D vertical is not implemented for MLHO"); + } + } + break; + case DragCoefficientAbsGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAbsGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAlongGradientEnum: + /*Nothing in P vector*/ + break; + case ThicknessAcrossGradientEnum: + /*Nothing in P vector*/ + break; + case RheologyBbarAbsGradientEnum: + /*Nothing in P vector*/ + break; + case RheologyBAbsGradientEnum: + /*Nothing in P vector*/ + break; + case RheologyBInitialguessMisfitEnum: + /*Nothing in P vector*/ + break; + default: + _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); + } + } + } + + /*Transform coordinate system*/ + //if(domaintype!=Domain2DverticalEnum) basalelement->TransformLoadVectorCoord(pe,XYEnum); + /*Clean up and return*/ + xDelete(responses); + xDelete(xyz_list); + xDelete(basis); + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; + delete gauss; + return pe; +}/*}}}*/ ElementVector* AdjointHorizAnalysis::CreatePVectorSSA(Element* element){/*{{{*/ @@ -1129,4 +1663,5 @@ case L1L2ApproximationEnum:GradientJDragL1L2(element,gradient,control_interp,control_index); break; case HOApproximationEnum: GradientJDragHO( element,gradient,control_interp,control_index); break; + case MLHOApproximationEnum: GradientJDragMLHO( element,gradient,control_interp,control_index); break; case FSApproximationEnum: GradientJDragFS( element,gradient,control_interp,control_index); break; case NoneApproximationEnum: /*Gradient is 0*/ break; @@ -1149,4 +1684,5 @@ case L1L2ApproximationEnum:GradientJBbarL1L2(element,gradient,control_interp,control_index); break; case HOApproximationEnum: GradientJBbarHO( element,gradient,control_interp,control_index); break; + case MLHOApproximationEnum: GradientJBbarMLHO( element,gradient,control_interp,control_index); break; case FSApproximationEnum: GradientJBbarFS( element,gradient,control_interp,control_index); break; case NoneApproximationEnum: /*Gradient is 0*/ break; @@ -1158,4 +1694,5 @@ case SSAApproximationEnum: GradientJBSSA(element,gradient,control_interp,control_index); break; case HOApproximationEnum: GradientJBHO( element,gradient,control_interp,control_index); break; + // case MLHOApproximationEnum: GradientJBMLHO( element,gradient,control_interp,control_index); break; case FSApproximationEnum: GradientJBFS( element,gradient,control_interp,control_index); break; case NoneApproximationEnum: /*Gradient is 0*/ break; @@ -1266,4 +1803,132 @@ /*WARNING: We use SSA as an estimate for now*/ this->GradientJBbarSSA(element,gradient,control_interp,control_index); +}/*}}}*/ +void AdjointHorizAnalysis::GradientJBbarMLHO(Element* element,Vector* gradient,int control_interp,int control_index){/*{{{*/ + + if(control_interp!=P1Enum) _error_("not implemented yet..."); + /*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(true); + dim = 2; + break; + default: _error_("mesh "<GetNumberOfVertices(); + + /*Initialize some vectors*/ + IssmDouble* basis = xNew(numvertices); + IssmDouble* ge = xNewZeroInit(numvertices); + int* vertexpidlist = xNew(numvertices); + IssmDouble zeta; + IssmDouble e1[3],e2[3], phishx[3], phishy[3]; + IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble adjshx, adjshy; + + /*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* n_input = basalelement->GetInput(MaterialsRheologyNEnum); _assert_(n_input); + Input* vxbase_input = basalelement->GetInput(VxBaseEnum); _assert_(vxbase_input); + Input* vybase_input = basalelement->GetInput(VyBaseEnum); _assert_(vybase_input); + Input* vxshear_input = basalelement->GetInput(VxShearEnum); _assert_(vxshear_input); + Input* vyshear_input = basalelement->GetInput(VyShearEnum); _assert_(vyshear_input); + + Input* adjointbx_input = basalelement->GetInput(AdjointxBaseEnum); _assert_(adjointbx_input); + Input* adjointby_input = basalelement->GetInput(AdjointyBaseEnum); _assert_(adjointby_input); + Input* adjointshx_input = basalelement->GetInput(AdjointxShearEnum); _assert_(adjointshx_input); + Input* adjointshy_input = basalelement->GetInput(AdjointyShearEnum); _assert_(adjointshy_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(5); + GaussSeg* gauss_seg=new GaussSeg(2); + + while(gauss->next()){ + + thickness_input->GetInputValue(&thickness,gauss); + n_input->GetInputValue(&n,gauss); + adjointbx_input->GetInputDerivativeValue(&dadjbx[0],xyz_list,gauss); + adjointby_input->GetInputDerivativeValue(&dadjby[0],xyz_list,gauss); + adjointshx_input->GetInputDerivativeValue(&dadjshx[0],xyz_list,gauss); + adjointshy_input->GetInputDerivativeValue(&dadjshy[0],xyz_list,gauss); + adjointshx_input->GetInputValue(&adjshx, gauss); + adjointshy_input->GetInputValue(&adjshy, gauss); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsP1(basis,gauss); + + /* Get the integration in the vertical direction */ + gauss_seg->Reset(); + while(gauss_seg->next()){ + zeta=0.5*(gauss_seg->coord1+1); + + basalelement->StrainRateMLHO(&epsilon[0],xyz_list,gauss, + vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); + + basalelement->dViscositydBMLHO(&dmudB,dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); + + e1[0] = 2.0*epsilon[0]+epsilon[1]; + e1[1] = epsilon[2]; + e1[2] = epsilon[3]; + + e2[0] = epsilon[2]; + e2[1] = epsilon[0]+2.0*epsilon[1]; + e2[2] = epsilon[4]; + + phishx[0] = dadjbx[0] + (1-pow(zeta, n+1))*dadjshx[0]; + phishx[1] = dadjbx[1] + (1-pow(zeta, n+1))*dadjshx[1]; + phishx[2] = (n+1)/thickness*pow(zeta, n)*adjshx; + phishy[0] = dadjby[0] + (1-pow(zeta, n+1))*dadjshy[0]; + phishy[1] = dadjby[1] + (1-pow(zeta, n+1))*dadjshy[1]; + phishy[2] = (n+1)/thickness*pow(zeta, n)*adjshy; + + /*Build gradient vector (actually -dJ/dB): */ + for(int i=0;iweight*basis[i]*thickness*0.5*gauss_seg->weight; + _assert_(!xIsNan(ge[i])); + } + } + } + } + if(control_interp==P1Enum){ + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + } + else{ + _error_("not supported"); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + delete gauss_seg; + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; }/*}}}*/ void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector* gradient,int control_interp,int control_index){/*{{{*/ @@ -1580,4 +2245,7 @@ delete gauss; }/*}}}*/ +void AdjointHorizAnalysis::GradientJBMLHO(Element* element,Vector* gradient,int control_interp,int control_index){/*{{{*/ + _error_("not implemented yet..."); +}/*}}}*/ void AdjointHorizAnalysis::GradientJBSSA(Element* element,Vector* gradient,int control_interp,int control_index){/*{{{*/ @@ -1977,4 +2645,122 @@ delete friction; }/*}}}*/ +void AdjointHorizAnalysis::GradientJDragMLHO(Element* element,Vector* gradient,int control_interp,int control_index){/*{{{*/ + + /*return if floating (gradient is 0)*/ + if(element->IsAllFloating()) return; + + /*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); + + /*Build friction element, needed later: */ + Friction* friction=new Friction(basalelement,dim); + + /*Retrieve all inputs we will be needing: */ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->GradientIndexing(&vertexpidlist[0],control_index); + Input* vx_input = basalelement->GetInput(VxBaseEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyBaseEnum); _assert_(vy_input); + Input* adjointx_input = basalelement->GetInput(AdjointxBaseEnum); _assert_(adjointx_input); + Input* adjointy_input = basalelement->GetInput(AdjointyBaseEnum); _assert_(adjointy_input); + + /* get the friction law: 1- Budd, 11-Schoof*/ + int frictionlaw;element->FindParam(&frictionlaw, FrictionLawEnum); + Input* dragcoeff_input = NULL; + switch(frictionlaw) { + case 1: + dragcoeff_input = basalelement->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input); + break; + case 2: + case 11: + dragcoeff_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoeff_input); + break; + default: + _error_("Friction law "<< frictionlaw <<" not supported in the inversion."); + } + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(4); + while(gauss->next()){ + + adjointx_input->GetInputValue(&lambda, gauss); + adjointy_input->GetInputValue(&mu, gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + dragcoeff_input->GetInputValue(&drag, gauss); + + friction->GetAlphaComplement(&dalpha2dk,gauss); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctionsP1(basis,gauss); + + /*Build gradient vector (actually -dJ/dD): */ + if(control_interp==P1Enum){ + for(int i=0;iweight*basis[i]; + _assert_(!xIsNan(ge[i])); + } + } + else if(control_interp==P0Enum){ + ge[0]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight; + _assert_(!xIsNan(ge[0])); + } + else{ + _error_("not supported"); + } + } + if(control_interp==P1Enum){ + gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL); + } + else if(control_interp==P0Enum){ + gradient->SetValue(vertexpidlist[0],ge[0],ADD_VAL); + } + else{ + _error_("not supported"); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(ge); + xDelete(vertexpidlist); + delete gauss; + delete friction; + if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; +}/*}}}*/ void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector* gradient,int control_interp,int control_index){/*{{{*/ @@ -2480,4 +3266,7 @@ InputUpdateFromSolutionFS(solution,element); } + else if (approximation==MLHOApproximationEnum) { + InputUpdateFromSolutionMLHO(solution, element); + } else{ InputUpdateFromSolutionHoriz(solution,element); @@ -2625,4 +3414,77 @@ xDelete(doflist); }/*}}}*/ +void AdjointHorizAnalysis::InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element){/*{{{*/ + int i; + int* doflist=NULL; + + int domaintype; + element->FindParam(&domaintype,DomainTypeEnum); + if (domaintype!=Domain2DhorizontalEnum) _error_("mesh "<GetNumberOfNodes(); + int numdof = numnodes * 4; + + /*Fetch dof list and allocate solution vectors*/ + element->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); + IssmDouble* values = xNew(numdof); + IssmDouble* lambdax = xNew(numnodes); + IssmDouble* lambday = xNew(numnodes); + IssmDouble* lambdabx = xNew(numnodes); + IssmDouble* lambdaby = xNew(numnodes); + IssmDouble* lambdashx = xNew(numnodes); + IssmDouble* lambdashy = xNew(numnodes); + IssmDouble* n = xNew(numnodes); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iTransformSolutionCoord(&values[0],XYEnum); + + element->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum,0.); + + /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ + for(i=0;i(lambdabx[i])) _error_("NaN found in solution vector"); + if(xIsInf(lambdabx[i])) _error_("Inf found in solution vector"); + if(xIsNan(lambdashx[i])) _error_("NaN found in solution vector"); + if(xIsInf(lambdashx[i])) _error_("Inf found in solution vector"); + /* adjoint for the surface velocity */ + lambdax[i] = lambdabx[i] + lambdashx[i]*(n[i]+1)/(n[i]+2); + + lambdaby[i] =values[i*4+2]; + lambdashy[i]=values[i*4+3]; + /*Check solution*/ + if(xIsNan(lambdaby[i])) _error_("NaN found in solution vector"); + if(xIsInf(lambdaby[i])) _error_("Inf found in solution vector"); + if(xIsNan(lambdashy[i])) _error_("NaN found in solution vector"); + if(xIsInf(lambdashy[i])) _error_("Inf found in solution vector"); + /* adjoint for the surface velocity */ + lambday[i] = lambdaby[i] + lambdashy[i]*(n[i]+1)/(n[i]+2); + } + + /*Add vx and vy as inputs to the tria element: */ + element->AddInput(AdjointxBaseEnum,lambdabx,element->GetElementType()); + element->AddInput(AdjointyBaseEnum,lambdaby,element->GetElementType()); + element->AddInput(AdjointxShearEnum,lambdashx,element->GetElementType()); + element->AddInput(AdjointyShearEnum,lambdashy,element->GetElementType()); + + element->AddInput(AdjointxEnum,lambdax,element->GetElementType()); + element->AddInput(AdjointyEnum,lambday,element->GetElementType()); + + /*Free ressources:*/ + xDelete(values); + xDelete(lambdax); + xDelete(lambday); + xDelete(lambdabx); + xDelete(lambdaby); + xDelete(lambdashx); + xDelete(lambdashy); + xDelete(doflist); +}/*}}}*/ void AdjointHorizAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ /*Default, do nothing*/ Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 26477) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.h (revision 26478) @@ -28,4 +28,6 @@ ElementMatrix* CreateKMatrixFS(Element* element); ElementMatrix* CreateKMatrixHO(Element* element); + ElementMatrix* CreateKMatrixMLHO(Element* element); + ElementMatrix* CreateKMatrixMLHOVerticalIntergrated(Element* element); ElementMatrix* CreateKMatrixL1L2(Element* element); ElementMatrix* CreateKMatrixSSA(Element* element); @@ -34,4 +36,5 @@ ElementVector* CreatePVectorL1L2(Element* element); ElementVector* CreatePVectorHO(Element* element); + ElementVector* CreatePVectorMLHO(Element* element); ElementVector* CreatePVectorSSA(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); @@ -42,8 +45,10 @@ void GradientJBbarL1L2(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJBbarHO(Element* element,Vector* gradient,int control_interp,int control_index); + void GradientJBbarMLHO(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJBbarSSA(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJBFS(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJBGradient(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJBHO(Element* element,Vector* gradient,int control_interp,int control_index); + void GradientJBMLHO(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJBSSA(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJDragFS(Element* element,Vector* gradient,int control_interp,int control_index); @@ -51,4 +56,5 @@ void GradientJDragL1L2(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJDragHO(Element* element,Vector* gradient,int control_interp,int control_index); + void GradientJDragMLHO(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJDragSSA(Element* element,Vector* gradient,int control_interp,int control_index); void GradientJDragHydroFS(Element* element,Vector* gradient,int control_interp,int control_index); @@ -60,4 +66,5 @@ void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element); void InputUpdateFromSolutionHoriz(IssmDouble* solution,Element* element); + void InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); }; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26478) @@ -807,4 +807,33 @@ case MatestarEnum: material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff); + break; + default: _error_("not supported"); + } + + /*Assign output pointer*/ + *pdmudB=dmudB; + +} +/*}}}*/ +void Element::dViscositydBMLHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vxbase_input,Input* vybase_input, Input* vxshear_input ,Input* vyshear_input,Input* thickness_input,Input* n_input, IssmDouble zeta){/*{{{*/ + + /*Intermediaries*/ + int materialstype; + IssmDouble dmudB; + IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz]; */ + IssmDouble eps_eff; + IssmDouble eps0=1.e-27; + + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */ + this->StrainRateMLHO(&epsilon[0],xyz_list,gauss, + vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); + eps_eff=sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] + + epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1] + eps0*eps0); + + /*Get viscosity*/ + materialstype=this->material->ObjectEnum(); + switch(materialstype){ + case MaticeEnum: + material->GetViscosity_B(&dmudB,eps_eff,gauss); break; default: _error_("not supported"); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26478) @@ -81,4 +81,5 @@ void dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void dViscositydBMLHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vxbase_input,Input* vybase_input, Input* vxshear_input ,Input* vyshear_input,Input* thickness_input,Input* n_input, IssmDouble zeta); void dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void Echo(); Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Materials/Material.h (revision 26478) @@ -52,4 +52,5 @@ virtual void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; virtual void ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input)=0; + virtual void ViscosityMLHOAdjoint(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input)=0; virtual void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0; virtual void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0; Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.cpp (revision 26478) @@ -570,4 +570,7 @@ _error_("not implemented yet"); }/*}}}*/ +void Matestar::ViscosityMLHOAdjoint(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input){/*{{{*/ + _error_("not implemented yet"); +}/*}}}*/ void Matestar::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/ _error_("not implemented yet"); Index: /issm/trunk-jpl/src/c/classes/Materials/Matestar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Materials/Matestar.h (revision 26478) @@ -77,4 +77,5 @@ void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input); + void ViscosityMLHOAdjoint(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input); void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 26478) @@ -762,7 +762,8 @@ + epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]); - /*Get viscosity at zeta - FIXME for now Bbar*/ - this->GetViscosityBar(&mu,epsilon_eff,gauss); - + /*Get viscosity at zeta */ +//?? need to use Bar for the current inversion +// this->GetViscosity(&mu,epsilon_eff,gauss); + this->GetViscosityBar(&mu, epsilon_eff,gauss); thickness_input->GetInputValue(&H, gauss); n_input->GetInputValue(&n,gauss); @@ -792,4 +793,72 @@ pviscosity[2]=mubar[2]; pviscosity[3]=mubar[3]; + + /*Clean up*/ + delete gauss_seg; +}/*}}}*/ +void Matice::ViscosityMLHOAdjoint(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vxbase_input,Input* vybase_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input){/*{{{*/ + + /* To compute the additional 5 terms in the viscosity appear in the adjoint equation*/ + /*Intermediaries*/ + IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz]; */ + IssmDouble epsilon_eff; + IssmDouble zeta,H,n; + IssmDouble f[9],F[9]; + IssmDouble mubar[9]; + IssmDouble mu; + int order=10; + + for(int i=0;i<9;++i) mubar[i]=0; + + GaussSeg* gauss_seg=new GaussSeg(order); + //IssmDouble eps_eff_averaged=0; + while(gauss_seg->next()){ + + /*Compute zeta for gauss_seg point (0=surface, 1=base)*/ + zeta=0.5*(gauss_seg->coord1+1); + + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy (for a given zeta)*/ + element->StrainRateMLHO(&epsilon[0],xyz_list,gauss, + vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input,zeta); + epsilon_eff=sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] + + epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]); + + this->GetViscosity(&mu,epsilon_eff,gauss); + /*the adjoint viscosity with zeta dependent term*/ + mu = mu /epsilon_eff/epsilon_eff; + + thickness_input->GetInputValue(&H, gauss); + n_input->GetInputValue(&n,gauss); + + /*Compute fi and Fi at zeta*/ + f[0]=1; + f[1]=(1-pow(zeta,n+1)); + f[2]=(1-pow(zeta,n+1))*(1-pow(zeta,n+1)); + f[3]=pow(zeta,2*n); // NOTE: this is different from the forward formulation + f[4]=(1-pow(zeta,n+1))*(1-pow(zeta,n+1))*(1-pow(zeta,n+1)); + f[5]=(1-pow(zeta,n+1))*(1-pow(zeta,n+1))*(1-pow(zeta,n+1))*(1-pow(zeta,n+1)); + f[6]=(1-pow(zeta,n+1))*pow(zeta,2*n); + f[7]=(1-pow(zeta,n+1))*(1-pow(zeta,n+1))*pow(zeta,2*n); + f[8]=pow(zeta,4*n); + + + F[0]=H; + F[1]=H*(n+1)/(n+2); + F[2]=2*H*(n+1)*(n+1)/( (2*n+3)*(n+2) ); + F[3]=H/(2*n+1); + F[4]=6*H*(n+1)*(n+1)*(n+1)/( (n+2)*(2*n+3)*(3*n+4) ); + F[5]=24*H*(n+1)*(n+1)*(n+1)*(n+1)/( (n+2)*(2*n+3)*(3*n+4)*(4*n+5) ); + F[6]=H*(n+1)/( (2*n+1)*(3*n+2) ); + F[7]=2*H*(n+1)*(n+1)/( (2*n+1)*(3*n+2)*(4*n+3) ); + F[8]=H/(4*n+1); + + /*Sum the viscosity*/ + for(int i=0;i<9;i++) { + mubar[i]+=(H/2)*gauss_seg->weight*mu*f[i]; + } + }//while + + /*Assign output pointer*/ + for(int i=0;i<9;i++) pviscosity[i]=mubar[i]; /*Clean up*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 26478) @@ -79,4 +79,5 @@ void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input); + void ViscosityMLHOAdjoint(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input); void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); Index: /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h (revision 26477) +++ /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h (revision 26478) @@ -73,4 +73,5 @@ void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; void ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input){_error_("not supported");}; + void ViscosityMLHOAdjoint(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input,Input* n_input){_error_("not supported");}; void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");}; void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; Index: /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 26477) +++ /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 26478) @@ -188,4 +188,5 @@ /*}}}*/ #elif defined(_HAVE_CODIPACK_) + /*{{{*/ simul_starttrace2(femmodel); IssmDouble* aX=xNew(n); @@ -252,5 +253,5 @@ /*Clear tape*/ tape_codi.reset(); - +/*}}}*/ #else /*{{{*/ Index: /issm/trunk-jpl/src/m/classes/inversion.m =================================================================== --- /issm/trunk-jpl/src/m/classes/inversion.m (revision 26477) +++ /issm/trunk-jpl/src/m/classes/inversion.m (revision 26478) @@ -104,6 +104,6 @@ %Only SSA, HO and FS are supported right now if strcmp(solution,'StressbalanceSolution') - if ~(md.flowequation.isSSA || md.flowequation.isHO || md.flowequation.isFS || md.flowequation.isL1L2), - md = checkmessage(md,['inversion can only be performed for SSA, HO or FS ice flow models']); + if ~(md.flowequation.isSSA || md.flowequation.isMLHO || md.flowequation.isHO || md.flowequation.isFS || md.flowequation.isL1L2), + md = checkmessage(md,['inversion can only be performed for SSA, MLHO, HO or FS ice flow models']); end end Index: /issm/trunk-jpl/test/NightlyRun/test333.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test333.m (revision 26478) +++ /issm/trunk-jpl/test/NightlyRun/test333.m (revision 26478) @@ -0,0 +1,34 @@ +%Test Name: SquareSheetConstrainedCMDragMLHO +md=triangle(model(),'../Exp/Square.exp',200000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/SquareSheetConstrained.par'); +md=setflowequation(md,'MLHO','all'); + +%control parameters +md.inversion.iscontrol=1; +md.inversion.control_parameters={'FrictionCoefficient'}; +md.inversion.min_parameters=1.*ones(md.mesh.numberofvertices,1); +md.inversion.max_parameters=200.*ones(md.mesh.numberofvertices,1); +md.inversion.nsteps=2; +md.inversion.cost_functions=[102 501]; +md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,2); md.inversion.cost_functions_coefficients(:,2)=2.*10^-7; +md.inversion.gradient_scaling=3.*ones(md.inversion.nsteps,1); +md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1); +md.inversion.step_threshold=0.3*ones(md.inversion.nsteps,1); +md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy; + +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,'Stressbalance'); + +%Fields and tolerances to track changes +field_names ={'Gradient','Misfits','FrictionCoefficient','Pressure','Vel','Vx','Vy'}; +field_tolerances={1e-12,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.StressbalanceSolution.Gradient1),... + (md.results.StressbalanceSolution.J),... + (md.results.StressbalanceSolution.FrictionCoefficient),... + (md.results.StressbalanceSolution.Pressure),... + (md.results.StressbalanceSolution.Vel),... + (md.results.StressbalanceSolution.Vx),... + (md.results.StressbalanceSolution.Vy) +}; Index: /issm/trunk-jpl/test/NightlyRun/test333.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test333.py (revision 26478) +++ /issm/trunk-jpl/test/NightlyRun/test333.py (revision 26478) @@ -0,0 +1,45 @@ +#Test Name: SquareSheetConstrainedCMDragMLHO +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * + + +md = triangle(model(), '../Exp/Square.exp', 200000.) +md = setmask(md, '', '') +md = parameterize(md, '../Par/SquareSheetConstrained.py') +md = setflowequation(md, 'MLHO', 'all') + +#control parameters + +md.inversion.iscontrol = 1 +md.inversion.control_parameters = ['FrictionCoefficient'] +md.inversion.min_parameters = 1. * np.ones((md.mesh.numberofvertices, len(md.inversion.control_parameters))) +md.inversion.max_parameters = 200. * np.ones((md.mesh.numberofvertices, len(md.inversion.control_parameters))) +md.inversion.nsteps = 2 +md.inversion.cost_functions = [102, 501] +md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, 2)) +md.inversion.cost_functions_coefficients[:, 1] = 2. * 10**-7 +md.inversion.gradient_scaling = 3. * np.ones((md.inversion.nsteps, len(md.inversion.control_parameters))) +md.inversion.maxiter_per_step = 2 * np.ones(md.inversion.nsteps) +md.inversion.step_threshold = 0.3 * np.ones(md.inversion.nsteps) +md.inversion.vx_obs = md.initialization.vx +md.inversion.vy_obs = md.initialization.vy + +md.cluster = generic('name', gethostname(), 'np', 3) +md = solve(md, 'Stressbalance') + +#Fields and tolerances to track changes +field_names = ['Gradient', 'Misfits', 'FrictionCoefficient', 'Pressure', 'Vel', 'Vx', 'Vy'] +field_tolerances = [1e-12, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13] +field_values = [md.results.StressbalanceSolution.Gradient1, + md.results.StressbalanceSolution.J, + md.results.StressbalanceSolution.FrictionCoefficient, + md.results.StressbalanceSolution.Pressure, + md.results.StressbalanceSolution.Vel, + md.results.StressbalanceSolution.Vx, + md.results.StressbalanceSolution.Vy] Index: /issm/trunk-jpl/test/NightlyRun/test334.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test334.m (revision 26478) +++ /issm/trunk-jpl/test/NightlyRun/test334.m (revision 26478) @@ -0,0 +1,35 @@ +%Test Name: SquareShelfCMBMLHO +md=triangle(model(),'../Exp/Square.exp',200000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelf.par'); +md=setflowequation(md,'MLHO','all'); + +%control parameters +md.inversion.iscontrol=1; +md.inversion.control_parameters={'MaterialsRheologyBbar'}; +md.inversion.min_parameters=10^6*ones(md.mesh.numberofvertices,1); +md.inversion.max_parameters=2*10^9*ones(md.mesh.numberofvertices,1); +md.inversion.nsteps=2; +md.inversion.cost_functions=101; +md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); +md.inversion.gradient_scaling=10^7*ones(md.inversion.nsteps,1); +md.inversion.maxiter_per_step=2*ones(md.inversion.nsteps,1); +md.inversion.step_threshold=0.3*ones(md.inversion.nsteps,1); +md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy; +md.verbose.control=true; + +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,'Stressbalance'); + +%Fields and tolerances to track changes +field_names ={'Gradient','Misfits','MaterialsRheologyBbar','Pressure','Vel','Vx','Vy'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.StressbalanceSolution.Gradient1),... + (md.results.StressbalanceSolution.J),... + (md.results.StressbalanceSolution.MaterialsRheologyBbar),... + (md.results.StressbalanceSolution.Pressure),... + (md.results.StressbalanceSolution.Vel),... + (md.results.StressbalanceSolution.Vx),... + (md.results.StressbalanceSolution.Vy) +}; Index: /issm/trunk-jpl/test/NightlyRun/test334.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test334.py (revision 26478) +++ /issm/trunk-jpl/test/NightlyRun/test334.py (revision 26478) @@ -0,0 +1,47 @@ +#Test Name: SquareShelfCMBMLHO +from model import * +from socket import gethostname +import numpy as np +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from generic import generic + +md = triangle(model(), '../Exp/Square.exp', 200000) +md = setmask(md, 'all', '') +md = parameterize(md, '../Par/SquareShelf.py') +md = setflowequation(md, 'MLHO', 'all') + +# control parameters + +md.inversion.iscontrol = 1 +md.inversion.control_parameters = ['MaterialsRheologyBbar'] +md.inversion.min_parameters = 1.0e6 * np.ones((md.mesh.numberofvertices, len(md.inversion.control_parameters))) +md.inversion.max_parameters = 2.0e9 * np.ones((md.mesh.numberofvertices, len(md.inversion.control_parameters))) +md.inversion.nsteps = 2 +md.inversion.cost_functions = [101] +md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, len(md.inversion.cost_functions))) +md.inversion.gradient_scaling = 1.0e7 * np.ones((md.inversion.nsteps, len(md.inversion.control_parameters))) +md.inversion.maxiter_per_step = 2. * np.ones((md.inversion.nsteps)) +md.inversion.step_threshold = 0.3 * np.ones((md.inversion.nsteps)) +md.inversion.vx_obs = md.initialization.vx +md.inversion.vy_obs = md.initialization.vy + + +md.cluster = generic('name', gethostname(), 'np', 3) +md = solve(md, 'Stressbalance') + + +# Fields and tolerances to track changes + +field_names = ['Gradient', 'Misfits', 'MaterialsRheologyBbar', 'Pressure', 'Vel', 'Vx', 'Vy'] +field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13] +field_values = [md.results.StressbalanceSolution.Gradient1, + md.results.StressbalanceSolution.J, + md.results.StressbalanceSolution.MaterialsRheologyBbar, + md.results.StressbalanceSolution.Pressure, + md.results.StressbalanceSolution.Vel, + md.results.StressbalanceSolution.Vx, + md.results.StressbalanceSolution.Vy]