Changeset 25836 for issm/trunk/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp

Timestamp:

12/08/20 08:45:53 (4 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 25834

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c (modified) (1 prop)
• src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp (modified) (4 diffs)

issm/trunk/src/c

@@ -23 +23 @@
 issm_ocean
 issm_dakota
+issm_post

@@ -23 +23 @@
 issm_ocean
 issm_dakota
+issm_post

issm/trunk/src/c/modules/InputUpdateFromDakotax/InputUpdateFromDakotax.cpp

@@ -6 +6 @@
 #include "../../shared/shared.h"
 #include "../../toolkits/toolkits.h"
+#include "../../classes/Inputs/TransientInput.h"
+#include "../../classes/Inputs/DatasetInput.h"
+#include "../../classes/Inputs/TriaInput.h"
 #include "../InputUpdateFromMatrixDakotax/InputUpdateFromMatrixDakotax.h"
 #include "../InputUpdateFromConstantx/InputUpdateFromConstantx.h"
 #include "../InputUpdateFromVectorDakotax/InputUpdateFromVectorDakotax.h"
 
-void InputUpdateFromDakotax(FemModel* femmodel,double* variables,char* *variables_descriptors,int numvariables){
+void  InputUpdateFromDakotax(FemModel* femmodel,double* variables,char* *variables_descriptors,int numdakotavariables){ /*{{{*/
 
         int     i,j,k,l;
-        int     dummy;
-
-        int     numberofvertices;
-        int     numberofelements;
-        int     nrows;
-        int     ncols;
-        int     npart;
-        double *qmu_vpart  = NULL;
-        double *qmu_epart  = NULL;
+
+        IssmDouble **variable_partitions         = NULL;
+        IssmDouble * variable_partition         = NULL;
+        int * variable_partitions_npart         = NULL;
+        int * variable_partitions_nt         = NULL;
+        int          variable_partitions_num;
+        int          npart;
+        int          nt;
+        int variablecount=0;
 
         double *distributed_values = NULL;
@@ -28 +31 @@
         char    root[50]; //root name of variable, ex: DragCoefficent, RhoIce, etc ...
 
+        if (VerboseQmu())_printf0_("dakota variables updates\n");
+
         /*retrieve parameters: */
-        femmodel->parameters->FindParam(&npart,QmuNumberofpartitionsEnum);
-        femmodel->parameters->FindParam(&qmu_vpart,NULL,QmuVpartitionEnum);
-        femmodel->parameters->FindParam(&qmu_epart,NULL,QmuEpartitionEnum);
-        numberofvertices=femmodel->vertices->NumberOfVertices();
-        numberofelements=femmodel->elements->NumberOfElements();
-
-        /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and 
-         * for each descriptor, take the variable value and plug it into the inputs: */
-
-        for(i=0;i<numvariables;i++){
+        femmodel->parameters->FindParam(&variable_partitions,&variable_partitions_num,NULL,NULL,QmuVariablePartitionsEnum);
+        femmodel->parameters->FindParam(&variable_partitions_npart,NULL,NULL,QmuVariablePartitionsNpartEnum);
+        femmodel->parameters->FindParam(&variable_partitions_nt,NULL,NULL,QmuVariablePartitionsNtEnum);
+
+
+        /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and
+         * for each descriptor, take the variable value and plug it into the inputs (more or less :)):
+         * We also start with distributed and standard values , as they tend to be used to pluck data from a multi-modle ensemble (mme)
+         * which can then be scaled. Doing the scaling first would be impractical, as the entire mme would have to be scaled,
+         * which is a waste of time:*/
+
+        variablecount=0;
+        for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions.
 
                 descriptor=variables_descriptors[i];
 
-                /*From descriptor, figure out if the variable is scaled, indexed, nodal, or just a simple variable: */
+
+                /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */
                 if (strncmp(descriptor,"scaled_",7)==0){
-
-                        /*Variable is scaled. Determine root name of variable (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the 
-                         * distributed_values with the next npart variables: */
-
-                        //strcpy(root,strstr(descriptor,"_")+1); *strstr(root,"_")='\0';
+                        /*we are skipping these for now.*/
+                        npart=variable_partitions_npart[variablecount];
+                        nt=variable_partitions_nt[variablecount];
+
+                        /*increment i to skip the distributed values just collected: */
+                        i+=npart*nt-1; //careful, the for loop will add 1.
+                }
+                else if (strncmp(descriptor,"indexed_",8)==0){
+                        /*we are skipping these for now.*/
+                }
+                else if (strncmp(descriptor,"nodal_",8)==0){
+                        /*we are skipping these for now.*/
+                }
+
+                else if (strncmp(descriptor,"distributed_",12)==0){
+                        if (VerboseQmu())_printf0_("   updating variable " << descriptor << "\n");
+
+                        /*recover partition vector: */
+                        variable_partition=variable_partitions[variablecount];
+                        npart=variable_partitions_npart[variablecount];
+
+                        /*Variable is distributed. Determine root name of variable (ex: distributed_DragCoefficient_1 -> DragCoefficient).
+                         * Allocate distributed_values and fill the distributed_values with the next npart variables: */
+
                         memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char));
                         *strstr(root,"_")='\0';
@@ -57 +85 @@
                         }
 
-                        /*Now, pick up the parameter corresponding to root: */
-                        femmodel->parameters->FindParamInDataset(&parameter,&nrows,&ncols,QmuVariableDescriptorsEnum,StringToEnumx(root));
-
-
-                        /*We've got the parameter, we need to update it using qmu_vpart or qmu_epart (depending on whether 
-                         * it is element or vertex distributed), and the distributed_values. 
-                         In addition, the parameter can be either a static or transient (nrows+1) vector. */
-                        
-                        //_printf_("nrows: " << nrows << " numberofvertices: " << numberofvertices << " numberofelements: " << numberofelements << "\n");
-
-                        if (nrows==numberofvertices || nrows==(numberofvertices+1)){
-                                for(k=0;k<numberofvertices;k++){
-                                        for(l=0;l<ncols;l++){
-                                                *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)qmu_vpart[k]];
-                                        }
-                                }
-                        }
-                        else if (nrows==numberofelements || nrows==(numberofelements+1)){
-                                for(k=0;k<numberofelements;k++){
-                                        for(l=0;l<ncols;l++){
-                                                *(parameter+ncols*k+l)=*(parameter+ncols*k+l)*distributed_values[(int)qmu_epart[k]];
-                                        }
-                                }
-
-                        }
-                        else _error_("partitioning vector should be either elements or vertex sized!");
-
-                        #ifdef _DEBUG_
-                                PetscSynchronizedPrintf(IssmComm::GetComm(),"Parameter matrix:");
-                                PetscSynchronizedFlush(IssmComm::GetComm());
-                                for(l=0;l<ncols;l++){
-                                        PetscSynchronizedPrintf(IssmComm::GetComm()," time %i\n",l);
-                                        PetscSynchronizedFlush(IssmComm::GetComm());
-
-                                        for(k=0;k<numberofvertices;k++){
-                                                PetscSynchronizedPrintf(IssmComm::GetComm()," node %i value %g\n",k+1,*(parameter+k*ncols+l));
-                                                PetscSynchronizedFlush(IssmComm::GetComm());
-                                        }
-                                }
-                                PetscSynchronizedPrintf(IssmComm::GetComm()," descriptor: %s root %s enum: %i\n",descriptor,root,StringToEnumx(root));
-                                PetscSynchronizedFlush(IssmComm::GetComm());
-                        #endif
-
-                        /*Update inputs using the parameter matrix: */
-                        if(nrows==numberofvertices || (nrows==numberofvertices+1))
-                                InputUpdateFromMatrixDakotax(femmodel, parameter, nrows,ncols,StringToEnumx(root), VertexEnum);
-                        else
-                                InputUpdateFromMatrixDakotax(femmodel, parameter, nrows,ncols,StringToEnumx(root), ElementEnum);
+                        //for (int j=0;j<npart;j++)_printf_(j << ":" << distributed_values[j] << "\n");
+
+                        //Call specialty code:
+                        InputUpdateSpecialtyCode(femmodel,distributed_values,variable_partition,npart,root);
 
                         /*increment i to skip the distributed values just collected: */
@@ -113 +97 @@
                         xDelete<double>(distributed_values);
                 }
-                else if (strncmp(descriptor,"indexed_",8)==0){
-                        _error_("indexed variables not supported yet!");
-                }
-                else if (strncmp(descriptor,"nodal_",8)==0){
-                        _error_("nodal variables not supported yet!");
-                }
                 else{
                         /*Ok, standard variable, just update inputs using the variable: */
+                        if (VerboseQmu())_printf0_("   updating variable " << descriptor << "\n");
                         InputUpdateFromConstantx(femmodel,variables[i],StringToEnumx(descriptor));
                 }
-        }
+                variablecount++;
+        }
+
+        variablecount=0;
+        /*now deal with scaled variabes:*/
+        for(i=0;i<numdakotavariables;i++){ //these are the dakota variables, for all partitions.
+
+                descriptor=variables_descriptors[i];
+
+                /*From descriptor, figure out if the variable is scaled, indexed, distributed or just a simple variable: */
+                if (strncmp(descriptor,"scaled_",7)==0){
+
+                        if (VerboseQmu())_printf0_("   updating variable " << descriptor << "\n");
+
+                        /*recover partition vector: */
+                        variable_partition=variable_partitions[variablecount];
+                        npart=variable_partitions_npart[variablecount];
+                        nt=variable_partitions_nt[variablecount];
+
+                        /* Variable is scaled, determine its root name (ex: scaled_DragCoefficient_1 -> DragCoefficient). Allocate distributed_values and fill the
+                         * distributed_values with the next npart variables coming from Dakota: */
+                        memcpy(root,strstr(descriptor,"_")+1,(strlen(strstr(descriptor,"_")+1)+1)*sizeof(char));
+                        *strstr(root,"_")='\0';
+
+                        distributed_values=xNew<double>(npart*nt);
+                        for(j=0;j<npart*nt;j++){
+                                distributed_values[j]=variables[i+j];
+                        }
+
+                        /*Scale variable inside the inputs:*/
+                        InputScaleFromDakotax(femmodel, distributed_values, variable_partition,npart, nt, StringToEnumx(root));
+
+                        /*increment i to skip the distributed values just collected: */
+                        i+=npart*nt-1; //careful, the for loop will add 1.
+
+                        /*Free allocations: */
+                        xDelete<double>(parameter);
+                        xDelete<double>(distributed_values);
+                }
+                variablecount++;
+        }
+
+        /*Save results:*/
+        femmodel->results->AddResult(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,"uq_variables",variables,numdakotavariables,1,1,0));
 
         /*Free ressources:*/
-        xDelete<double>(qmu_vpart);
-        xDelete<double>(qmu_epart);
-}
+        for(i=0;i<variable_partitions_num;i++){
+                IssmDouble* matrix=variable_partitions[i];
+                xDelete<IssmDouble>(matrix);
+        }
+        xDelete<IssmDouble*>(variable_partitions);
+        xDelete<int>(variable_partitions_npart);
+        xDelete<int>(variable_partitions_nt);
+
+} /*}}}*/
+void  InputUpdateSpecialtyCode(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* variable_partition,int npart,char* root){ //{{{
+
+        /*Here, we put all the code that cannot be handled any other place: */
+        if (strncmp(root,"SurfaceloadModelid",18)==0){
+
+                if(VerboseQmu()){
+                        _printf0_("      SurfaceloadModelid MME, with ids: ");
+                        for (int i=0;i<npart;i++)_printf0_((int)distributed_values[i] << " ");
+                        _printf0_("\n");
+                }
+
+
+                if (femmodel->inputs->GetInputObjectEnum(SurfaceloadIceThicknessChangeEnum)==DatasetInputEnum)
+                        MmeToInput(femmodel,distributed_values,variable_partition,npart,SurfaceloadIceThicknessChangeEnum, P0Enum);
+
+                if (femmodel->inputs->GetInputObjectEnum(MaskIceLevelsetEnum)==DatasetInputEnum)
+                        MmeToInput(femmodel,distributed_values,variable_partition,npart,MaskIceLevelsetEnum, P1Enum);
+
+                if (femmodel->inputs->GetInputObjectEnum(MaskOceanLevelsetEnum)==DatasetInputEnum)
+                        MmeToInput(femmodel,distributed_values,variable_partition,npart,MaskOceanLevelsetEnum, P1Enum);
+
+        }
+        else _error_("InputUpdateSpecialtyCode error message: " << root << " not supported yet!");
+
+}       //}}}
+void  MmeToInput(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* variable_partition,int npart,int rootenum, int interpolationenum){ /*{{{*/
+
+        TransientInput* transientinput  = NULL;
+        TransientInput* transientinput2 = NULL;
+        Tria* element                    = NULL;
+        IssmDouble value;
+        IssmDouble* values               = NULL;
+        IssmDouble* times                = NULL;
+        int N;
+        int id;
+
+        /*find thickness dataset: */
+        DatasetInput* datasetinput = femmodel->inputs->GetDatasetInput(rootenum);
+
+        /*Initialize new transient input: */
+        transientinput = datasetinput->GetTransientInputByOffset(0); _assert_(transientinput);
+        transientinput->GetAllTimes(&times,&N);
+        femmodel->inputs->SetTransientInput(DummyEnum,times,N);
+        transientinput2 = femmodel->inputs->GetTransientInput(DummyEnum);
+
+        for(Object* & object : femmodel->elements->objects){
+                Tria*   element=xDynamicCast<Tria*>(object);
+
+                if((int)variable_partition[element->Sid()]==-1)id=0; //grab background field
+                else id=distributed_values[(int)variable_partition[element->Sid()]]-1; //grab partition field
+
+                /*recover the right field from the mme: */
+                transientinput = datasetinput->GetTransientInputByOffset(id); _assert_(transientinput);
+
+                /*copy values from the transientinput to the final transientinput2: */
+                for (int j=0;j<N;j++){
+                        TriaInput* tria_input=transientinput->GetTriaInput(j);
+                        element->InputServe(tria_input);
+                        if(interpolationenum==P0Enum){
+                                value=tria_input->element_values[0];
+                                transientinput2->AddTriaTimeInput( j,1,&(element->lid),&value,P0Enum);
+                        }
+                        else if(interpolationenum==P1Enum){
+
+                                /*Get values and lid list*/
+                                const int   numvertices     = element->GetNumberOfVertices();
+                                int        *vertexlids      = xNew<int>(numvertices);
+                                int        *vertexsids      = xNew<int>(numvertices);
+
+                                /*Recover vertices ids needed to initialize inputs*/
+                                element->GetVerticesLidList(&vertexlids[0]);
+                                element->GetVerticesSidList(&vertexsids[0]);
+                                values=tria_input->element_values;
+                                transientinput2->AddTriaTimeInput( j,numvertices,vertexlids,values,P1Enum);
+                        }
+                }
+        }
+
+        /*wipe out existing SurfaceloadIceThicknessChangeEnum dataset:*/
+        femmodel->inputs->ChangeEnum(DummyEnum,rootenum);
+
+        //reconfigure:
+        transientinput2->Configure(femmodel->parameters);
+}       //}}}
+void  InputScaleFromDakotax(FemModel* femmodel,IssmDouble* distributed_values,IssmDouble* partition, int npart, int nt, int name){ /*{{{*/
+
+        /*Copy input:*/
+        femmodel->inputs->DuplicateInput(name,DummyEnum);
+
+        /*Go through elements, copy input name to dummy, and scale it using the distributed_values and the partition vector:*/
+        for(Object* & object : femmodel->elements->objects){
+                Element* element=xDynamicCast<Element*>(object);
+                element->InputScaleFromDakota(distributed_values,partition,npart,nt,name);
+        }
+
+        /*We created a dummy input, which was a scaled copy of the name input. Now wipe
+         * out the name input with the new input:*/
+        femmodel->inputs->ChangeEnum(DummyEnum,name);
+
+        /*Some specialty code:*/
+        switch(name){
+                case MaterialsRheologyBEnum:
+                        femmodel->inputs->DuplicateInput(name,MaterialsRheologyBbarEnum);
+                        break;
+        }
+
+
+} /*}}}*/

@@ -23 +23 @@
 issm_ocean
 issm_dakota
+issm_post