Changeset 22178

Timestamp:

10/20/17 14:03:10 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on balancethickness2 with Jerome, new model, new implementation, new adjoint

Location:

issm/trunk-jpl/src/c

Files:

• analyses/AdjointBalancethickness2Analysis.cpp (modified) (3 diffs)
• analyses/Balancethickness2Analysis.cpp (modified) (5 diffs)
• analyses/Balancethickness2Analysis.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (1 diff)
• cores/adjointbalancethickness2_core.cpp (modified) (1 diff)
• cores/balancethickness2_core.cpp (modified) (3 diffs)
• cores/controlvalidation_core.cpp (modified) (1 diff)

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

@@ -192 +192 @@
 
         /*Intermediaries*/
-        IssmDouble dlambda[2],ds[2],D0,omega,Jdet; 
+        int         n=3;
+        IssmDouble  dlambda[2],ds[2],slopex,slopey,slope,omega,Jdet; 
         IssmDouble *xyz_list= NULL;
 
@@ -207 +208 @@
         element->GradientIndexing(&vertexpidlist[0],control_index);
         Input* adjoint_input = element->GetInput(AdjointEnum);            _assert_(adjoint_input);
-        Input* s_input       = element->GetInput(SurfaceEnum);            _assert_(s_input);
-        Input* D0_input      = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input);
         Input* omega_input   = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
+        Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); 
+        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); 
+        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); 
+        IssmDouble rhog    = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
         Gauss* gauss=element->NewGauss(2);
@@ -218 +221 @@
                 element->NodalFunctionsP1(basis,gauss);
 
-                D0_input->GetInputValue(&D0,gauss);
                 omega_input->GetInputValue(&omega,gauss);
                 adjoint_input->GetInputDerivativeValue(&dlambda[0],xyz_list,gauss);
-                s_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+                surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
+                surfaceslopex_input->GetInputValue(&slopex,gauss);
+                surfaceslopey_input->GetInputValue(&slopey,gauss);
+
+                slope = sqrt(slopex*slopex + slopey*slopey);
 
                 /*Build gradient vector (actually -dJ/da): */
                 for(int i=0;i<numvertices;i++){
-                        ge[i]+= -Jdet*gauss->weight*basis[i]*exp(omega)*D0*(ds[0]*dlambda[0] + ds[1]*dlambda[1]);
+                        ge[i]+= -Jdet*gauss->weight*basis[i]*pow(rhog,n)*exp(omega)*pow(slope,n-1)*(ds[0]*dlambda[0] + ds[1]*dlambda[1]);
                         _assert_(!xIsNan<IssmDouble>(ge[i]));
                 }

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

@@ -39 +39 @@
         iomodel->FetchDataToInput(elements,"md.balancethickness.thickening_rate",BalancethicknessThickeningRateEnum);
         iomodel->FetchDataToInput(elements,"md.balancethickness.omega",BalancethicknessOmegaEnum);
+        iomodel->FetchDataToInput(elements,"md.balancethickness.slopex",SurfaceSlopeXEnum);
+        iomodel->FetchDataToInput(elements,"md.balancethickness.slopey",SurfaceSlopeYEnum);
 
         /*Update elements: */
@@ -63 +65 @@
         return NULL;
 }/*}}}*/
-void           Balancethickness2Analysis::CreateD0(Element* element){/*{{{*/
-
-        /*Intermediaries */
-        IssmDouble       Gamma,h,mu0,ds[2],Cmu,B,k,s,b,normds;
-        const int        n = 3;
-        const IssmDouble Hstar = 500.;
-        const IssmDouble Lstar = 500.e+3;
-        IssmDouble *xyz_list  = NULL;
-
-        /*Fetch number of vertices and allocate output*/
-        int  numnodes = element->GetNumberOfNodes();
-        IssmDouble* D0     = xNew<IssmDouble>(numnodes);
-
-        /*retrieve what we need: */
-        element->GetVerticesCoordinates(&xyz_list);
-        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum);
-        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum);
-        Input* surface_input       = element->GetInput(SurfaceEnum);            _assert_(surface_input);
-        Input* B_input             = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
-        IssmDouble rhog            = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
-
-        /*Calculate damage evolution source term: */
-        Gauss* gauss=element->NewGauss();
-        for (int i=0;i<numnodes;i++){
-                gauss->GaussNode(element->GetElementType(),i);
-                
-                B_input->GetInputValue(&B,gauss);
-                if(surfaceslopex_input && surfaceslopey_input){
-                        surfaceslopex_input->GetInputValue(&ds[0],gauss);
-                        surfaceslopey_input->GetInputValue(&ds[1],gauss);
-                }
-                else{
-                        surface_input->GetInputDerivativeValue(&ds[0],xyz_list,gauss);
-                }
-
-                /*check slopes*/
-                normds = sqrt(ds[0]*ds[0]+ds[1]*ds[1]);
-                if (normds==0.){
-                        _error_("surface slope is zero");
-                }
-                if(normds<1.e-5){
-                        ds[0] = ds[0]/normds*1.e+5;
-                        ds[1] = ds[1]/normds*1.e+5;
-                        normds = 1.e-5;
-                }
-
-                mu0   = pow(2.,(1-3*n)/(2.*n)) * B;
-                Gamma = pow(rhog,n) * pow(Hstar,2*(n+1)) * pow(Hstar/Lstar,2*(n+1)) * 1./pow(mu0,n);
-
-                D0[i] = Gamma*pow(ds[0]*ds[0]+ds[1]*ds[1],(n-1)/2)/(n+2);
-        }
-
-        /*Add input*/
-        element->AddInput(BalancethicknessD0Enum,D0,element->GetElementType());
-        //if(surfaceslopex_input && surfaceslopey_input){
-        //      element->DeleteInput(SurfaceSlopeXEnum);
-        //      element->DeleteInput(SurfaceSlopeYEnum);
-        //}
-        
-        /*Clean up and return*/
-        xDelete<IssmDouble>(D0);
-        xDelete<IssmDouble>(xyz_list);
-        delete gauss;
-}/*}}}*/
 ElementMatrix* Balancethickness2Analysis::CreateJacobianMatrix(Element* element){/*{{{*/
 _error_("Not implemented");
@@ -133 +71 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,D0,omega;
+        int         n=3;
+        IssmDouble  Jdet,omega,ds[2],slope;
         IssmDouble* xyz_list = NULL;
 
@@ -145 +84 @@
         /*Retrieve all inputs and parameters*/
         element->GetVerticesCoordinates(&xyz_list);
-        Input* omega_input = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
-        Input* D0_input    = element->GetInput(BalancethicknessD0Enum);
-        if(!D0_input){
-                this->CreateD0(element);
-                D0_input = element->GetInput(BalancethicknessD0Enum); _assert_(D0_input);
-        }
+        Input* omega_input         = element->GetInput(BalancethicknessOmegaEnum); _assert_(omega_input);
+        Input* surfaceslopex_input = element->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input); 
+        Input* surfaceslopey_input = element->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input); 
+        IssmDouble rhog    = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
 
         /* Start  looping on the number of gaussian points: */
@@ -158 +95 @@
                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-                D0_input->GetInputValue(&D0,gauss);
                 omega_input->GetInputValue(&omega,gauss);
+                surfaceslopex_input->GetInputValue(&ds[0],gauss);
+                surfaceslopey_input->GetInputValue(&ds[1],gauss);
+
+                slope = sqrt(ds[0]*ds[0] + ds[1]*ds[1]);
 
                 for(int i=0;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
-                                Ke->values[i*numnodes+j] += D0*exp(omega)*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+                                Ke->values[i*numnodes+j] += pow(rhog,n)*exp(omega)*pow(slope,n-1)*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
                         }
                 }

issm/trunk-jpl/src/c/analyses/Balancethickness2Analysis.h

@@ -23 +23 @@
                 void           Core(FemModel* femmodel);
                 ElementVector* CreateDVector(Element* element);
-                void           CreateD0(Element* element);
                 ElementMatrix* CreateJacobianMatrix(Element* element);
                 ElementMatrix* CreateKMatrix(Element* element);

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -1759 +1759 @@
                                                 for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
                                                 this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                        }
+                                        break;
+                                case BalancethicknessOmegaEnum:
+                                        if(iomodel->Data("md.balancethickness.omega")){
+                                                for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data("md.balancethickness.omega")[tria_vertex_ids[j]-1];
+                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
+                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(tria_vertex_ids[j]-1)*num_control_type+i];
+                                                this->inputs->AddInput(new ControlInput(BalancethicknessOmegaEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
                                         }
                                         break;

issm/trunk-jpl/src/c/cores/adjointbalancethickness2_core.cpp

@@ -24 +24 @@
 
         /*Call SurfaceAreax, because some it might be needed by PVector*/
-        SurfaceAreax(NULL,femmodel);
+        //SurfaceAreax(NULL,femmodel);
 
         /*compute adjoint*/

issm/trunk-jpl/src/c/cores/balancethickness2_core.cpp

@@ -13 +13 @@
 
         /*parameters: */
-        bool        save_results;
-        IssmDouble  l = 3.;
+        bool save_results;
+        //IssmDouble  l = 3.;
 
         /*recover parameters: */
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
 
-        if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n");
+        //if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n");
         //femmodel->parameters->SetParam(l,SmoothThicknessMultiplierEnum);
         //femmodel->SetCurrentConfiguration(SmoothAnalysisEnum);
@@ -26 +26 @@
         //femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToSmoothEnum);
         //solutionsequence_linear(femmodel);
-        surfaceslope_core(femmodel);
+        //surfaceslope_core(femmodel);
 
         if(VerboseSolution()) _printf0_("call computational core:\n");
@@ -35 +35 @@
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
-                const int numoutputs = 6;
-                int outputs[numoutputs] = {SurfaceEnum,SurfaceSlopeXEnum,SurfaceSlopeYEnum,VxEnum,VyEnum,VelEnum};
-                //const int numoutputs = 4;
-                //int outputs[numoutputs] = {SurfaceEnum,VxEnum,VyEnum,VelEnum};
+                const int numoutputs = 1;
+                int outputs[numoutputs] = {SurfaceEnum};
                 femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs);
         }

issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp

@@ -104 +104 @@
         xDelete<IssmDouble>(X);
         xDelete<IssmDouble>(X0);
+        xDelete<IssmDouble>(scaling_factors);
 }