Changeset 17969

Timestamp:

05/09/14 10:22:13 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: preparing inversions of flowband models with Hongju

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointHorizAnalysis.cpp (modified) (5 diffs)
• modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (modified) (2 diffs)
• modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

@@ -208 +208 @@
         /*Intermediaries */
         bool        incomplete_adjoint;
+        int         dim,epssize;
         IssmDouble  Jdet,mu_prime;
         IssmDouble  eps1dotdphii,eps1dotdphij,eps2dotdphii,eps2dotdphij,eps3dotdphii,eps3dotdphij;
@@ -213 +214 @@
         IssmDouble *xyz_list = NULL;
 
+        /*Get problem dimension*/
+        element->FindParam(&dim,DomainDimensionEnum);
+        if(dim==2) epssize = 3;
+        else       epssize = 6;
+
         /*Fetch number of nodes and dof for this finite element*/
         int vnumnodes = element->NumberofNodesVelocity();
         int pnumnodes = element->NumberofNodesPressure();
-        int numdof    = vnumnodes*3 + pnumnodes;
+        int numdof    = vnumnodes*dim + pnumnodes;
 
         /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
@@ -229 +235 @@
         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_(vz_input);
+        Input* vz_input = NULL;
+        if(dim==3){
+                vz_input = element->GetInput(VzEnum);
+        }
+        else{
+                _error_("Not implemented yet");
+        }
 
         /*Allocate dbasis*/
-        IssmDouble* dbasis = xNew<IssmDouble>(3*vnumnodes);
+        IssmDouble* dbasis = xNew<IssmDouble>(dim*vnumnodes);
 
         /* Start  looping on the number of gaussian points: */
@@ -376 +388 @@
                                         for(i=0;i<vnumnodes;i++){
                                                 dux=vxobs-vx;
-                                                duy=vyobs-vy;
                                                 pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
-                                                pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+                                                if(dim==3){
+                                                        duy=vyobs-vy;
+                                                        pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+                                                }
                                         }
                                         break;
@@ -415 +429 @@
                                          */
                                         for(i=0;i<vnumnodes;i++){
-                                                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*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
-                                                pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+                                                if(dim==3){
+                                                        velocity_mag    =sqrt(vx*vx+vy*vy)+epsvel;
+                                                        obs_velocity_mag=sqrt(vxobs*vxobs+vyobs*vyobs)+epsvel;
+                                                        scale=-8.*meanvel*meanvel/(velocity_mag*velocity_mag)*log(velocity_mag/obs_velocity_mag);
+                                                        dux=scale*vx;
+                                                        duy=scale*vy;
+                                                        pe->values[i*dim+0]+=dux*weight*Jdet*gauss->weight*vbasis[i]; 
+                                                        pe->values[i*dim+1]+=duy*weight*Jdet*gauss->weight*vbasis[i]; 
+                                                }
+                                                else{
+                                                        velocity_mag    =fabs(vx)+epsvel;
+                                                        obs_velocity_mag=fabs(vxobs)+epsvel;
+                                                        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]; 
+                                                }
                                         }
                                         break;

issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp

@@ -21 +21 @@
         for (i=0;i<elements->Size();i++){
                 element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                J+=element->SurfaceAbsVelMisfit();
+                J+=SurfaceAbsVelMisfit(element);
         }
 
@@ -32 +32 @@
         *pJ=J;
 }
+
+IssmDouble SurfaceAbsVelMisfit(Element* element){
+
+        int        domaintype,numcomponents;
+        bool       surfaceintegral;
+        IssmDouble Jelem=0.;
+        IssmDouble misfit,Jdet;
+        IssmDouble vx,vy,vxobs,vyobs,weight;
+        IssmDouble* xyz_list_top = NULL;
+
+        /*Get basal element*/
+        if(!element->IsOnSurface()) return 0.;
+
+        /*If on water, return 0: */
+        if(!element->IsIceInElement()) return 0.;
+
+        /*Get problem dimension*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DverticalEnum:
+                        surfaceintegral = true;
+                        numcomponents   = 1;
+                        break;
+                case Domain3DEnum:
+                        surfaceintegral = true;
+                        numcomponents   = 2;
+                        break;
+                case Domain2DhorizontalEnum:
+                        surfaceintegral = false;
+                        numcomponents   = 2;
+                        break;
+                default: _error_("not supported yet");
+        }
+
+        /* Get node coordinates*/
+        if(surfaceintegral) element->GetVerticesCoordinatesTop(&xyz_list_top);
+        else                element->GetVerticesCoordinates(&xyz_list_top);
+
+        /*Retrieve all inputs we will be needing: */
+        Input* weights_input=element->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
+        Input* vx_input     =element->GetInput(VxEnum);        _assert_(vx_input);
+        Input* vxobs_input  =element->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
+        Input* vy_input     = NULL;
+        Input* vyobs_input  = NULL;
+        if(numcomponents==2){
+                vy_input    =element->GetInput(VyEnum);              _assert_(vy_input);
+                vyobs_input =element->GetInput(InversionVyObsEnum);  _assert_(vyobs_input);
+        }
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=NULL;
+        if(surfaceintegral) gauss=element->NewGaussTop(2);
+        else                gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                /* Get Jacobian determinant: */
+                if(surfaceintegral) element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
+                else                element->JacobianDeterminant(&Jdet,xyz_list_top,gauss);
+
+                /*Get all parameters at gaussian point*/
+                weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
+                vx_input->GetInputValue(&vx,gauss);
+                vxobs_input->GetInputValue(&vxobs,gauss);
+                if(numcomponents==2){
+                        vy_input->GetInputValue(&vy,gauss);
+                        vyobs_input->GetInputValue(&vyobs,gauss);
+                }
+
+                /*Compute SurfaceAbsVelMisfitEnum:
+                 *
+                 *      1  [           2              2 ]
+                 * J = --- | (u - u   )  +  (v - v   )  |
+                 *      2  [       obs            obs   ]
+                 *
+                 */
+                misfit=0.5*(vx-vxobs)*(vx-vxobs);
+                if(numcomponents==2) misfit+=0.5*(vy-vyobs)*(vy-vyobs);
+
+                /*Add to cost function*/
+                Jelem+=misfit*weight*Jdet*gauss->weight;
+        }
+
+        /*clean up and Return: */
+        delete gauss;
+        return Jelem;
+}

issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h

@@ -10 +10 @@
 /* local prototypes: */
 void SurfaceAbsVelMisfitx( IssmDouble* pJ, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters);
+IssmDouble SurfaceAbsVelMisfit(Element* element);
 
 #endif