Changeset 27721

Timestamp:

05/03/23 13:41:37 (23 months ago)

Author:

Mathieu Morlighem

Message:

CHG: done with transient cost function

Location:

issm/trunk-jpl/src/c/classes

Files:

• Cfsurfacesquaretransient.cpp (modified) (6 diffs)
• Cfsurfacesquaretransient.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.cpp

@@ -32 +32 @@
         this->datatimes      = NULL;
         this->passedflags    = NULL;
-}
-/*}}}*/
-Cfsurfacesquaretransient::Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum, int in_num_datatimes, IssmDouble* in_datatimes, bool* in_passedflags){/*{{{*/
+        this->J              = 0.;
+}
+/*}}}*/
+Cfsurfacesquaretransient::Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum, int in_num_datatimes, IssmDouble* in_datatimes, bool* in_passedflags, IssmDouble in_J){/*{{{*/
 
         this->definitionenum=in_definitionenum;
@@ -50 +51 @@
         xMemCpy<IssmDouble>(this->datatimes,in_datatimes,this->num_datatimes);
         xMemCpy<bool>(this->passedflags,in_passedflags,this->num_datatimes);
+
+        this->J = in_J;
 }
 /*}}}*/
@@ -70 +73 @@
         /*initialize passedtimes to false*/
         for(int i=0;i<this->num_datatimes;i++) this->passedflags[i]= false;
+        this->J = 0;
 }
 /*}}}*/
@@ -81 +85 @@
 /*Object virtual function resolutoin: */
 Object* Cfsurfacesquaretransient::copy() {/*{{{*/
-        Cfsurfacesquaretransient* output = new Cfsurfacesquaretransient(this->name,this->definitionenum, this->model_enum, this->num_datatimes, this->datatimes,this->passedflags);
+        Cfsurfacesquaretransient* output = new Cfsurfacesquaretransient(this->name,this->definitionenum, this->model_enum, this->num_datatimes, this->datatimes,this->passedflags, this->J);
         return (Object*)output;
 }
@@ -110 +114 @@
         marshallhandle->call(this->datatimes,this->num_datatimes);
         marshallhandle->call(this->passedflags,this->num_datatimes);
+        marshallhandle->call(this->J);
 } 
 /*}}}*/
@@ -130 +135 @@
 /*}}}*/
 IssmDouble Cfsurfacesquaretransient::Response(FemModel* femmodel){/*{{{*/
-         _error_("Not implemented yet");
-
-}
-/*}}}*/
+
+        /*recover model time parameters: */
+        IssmDouble time;
+        femmodel->parameters->FindParam(&time,TimeEnum);
+
+        /*Find closest datatime that is less than time*/
+        int pos=-1;
+        for(int i=0;i<this->num_datatimes;i++){
+                if(this->datatimes[i]<=time){
+                        pos = i;
+                }
+                else{
+                        break;
+                }
+        }
+
+        /*if pos=-1, time is earlier than the first data observation in this dataset*/
+        if(pos==-1){
+                _assert_(this->J==0.);
+                return 0.;
+        }
+
+        /*Check that we have not yet calculated this cost function*/
+        if(this->passedflags[pos]){
+                return 0.;
+                //return this->J; //FIXME
+        }
+
+        /*Calculate cost function for this time slice*/
+        IssmDouble J_part=0.;
+        for(Object* & object : femmodel->elements->objects){
+                Element* element=xDynamicCast<Element*>(object);
+                J_part+=this->Cfsurfacesquaretransient_Calculation(element,model_enum);
+        }
+
+        /*Sum across partition*/
+        IssmDouble J_sum;
+        ISSM_MPI_Allreduce( (void*)&J_part,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+        /*Record this cost function so that we do not recalculate it later*/
+        this->passedflags[pos]= true;
+        this->J += J_sum;
+
+        /*Return full cost function this far*/
+        //return this->J; //FIXME
+        return J_sum;
+}/*}}}*/
 IssmDouble Cfsurfacesquaretransient::Cfsurfacesquaretransient_Calculation(Element* element, int model_enum){/*{{{*/
 
-        _error_("not implemented");
+        IssmDouble Jelem=0.;
+        IssmDouble misfit,Jdet;
+        IssmDouble model,obs,weight;
+        IssmDouble* xyz_list = NULL;
+
+        /*Get basal element*/
+        if(!element->IsOnSurface()) return 0.;
+
+        /*If on water, return 0: */
+        if(!element->IsIceInElement()) return 0.;
+
+        /*Spawn surface element*/
+        Element* topelement = element->SpawnTopElement();
+
+        /* Get node coordinates*/
+        topelement->GetVerticesCoordinates(&xyz_list);
+
+        /*Retrieve all inputs we will be needing: */
+        DatasetInput *datasetinput = topelement->GetDatasetInput(definitionenum); _assert_(datasetinput);
+        Input        *model_input  = topelement->GetInput(model_enum);            _assert_(model_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=topelement->NewGauss(2);
+        while(gauss->next()){
+
+                /* Get Jacobian determinant: */
+                topelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+
+                /*Get all parameters at gaussian point*/
+                datasetinput->GetInputValue(&weight,gauss,WeightsSurfaceObservationEnum);
+                datasetinput->GetInputValue(&obs,gauss,SurfaceObservationEnum);
+                model_input->GetInputValue(&model,gauss);
+
+                /*Compute Misfit
+                 *     
+                 *       1  [           2 ]
+                 *  J = --- | (x - x   )  |
+                 *       2  [       obs   ]
+                 **/
+                misfit=0.5*(model-obs)*(model-obs);
+
+                /*Add to cost function*/
+                Jelem+=misfit*weight*Jdet*gauss->weight;
+        }
+
+        /*clean up and Return: */
+        if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;};
+        xDelete<IssmDouble>(xyz_list);
+        delete gauss;
+        return Jelem;
 }/*}}}*/

issm/trunk-jpl/src/c/classes/Cfsurfacesquaretransient.h

@@ -22 +22 @@
                 IssmDouble *datatimes;
                 bool       *passedflags;
+                IssmDouble  J;
 
                 /*Cfsurfacesquaretransient constructors, destructors :*/
                 Cfsurfacesquaretransient();
                 Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum,int num_datatimes, IssmDouble* in_datatime);
-                Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum,int num_datatimes, IssmDouble* in_datatime, bool* in_timepassedflag);
+                Cfsurfacesquaretransient(char* in_name, int in_definitionenum, int in_model_enum,int num_datatimes, IssmDouble* in_datatime, bool* in_timepassedflag, IssmDouble in_J);
                 ~Cfsurfacesquaretransient();