Changeset 765


Ignore:
Timestamp:
06/04/09 12:17:49 (15 years ago)
Author:
Eric.Larour
Message:

New qmu in core capability. New processing of parameters

Location:
issm/trunk/src
Files:
2 added
24 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r667 r765  
    109109                                        ./shared/Elements/Paterson.cpp\
    110110                                        ./shared/Elements/GetElementNodeData.cpp\
     111                                        ./shared/String/isdistributed.cpp\
     112                                        ./shared/String/sharedstring.h\
    111113                                        ./toolkits/petsc\
    112114                                        ./toolkits/petsc/patches\
    … …  
    376378                                        ./shared/Elements/Paterson.cpp\
    377379                                        ./shared/Elements/GetElementNodeData.cpp\
     380                                        ./shared/String/isdistributed.cpp\
     381                                        ./shared/String/sharedstring.h\
    378382                                        ./toolkits/petsc\
    379383                                        ./toolkits/petsc/patches\
    … …  
    470474                                        ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\
    471475                                        ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\
     476                                        ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\
    472477                                        ./Dofx/Dofx.h\
    473478                                        ./Dofx/Dofx.cpp\
    … …  
    557562bin_PROGRAMS =
    558563else
    559 bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe
     564dnl bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe
     565bin_PROGRAMS = diagnostic.exe 
    560566endif
    561567

  • issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp

    r762 r765  
    113113
    114114        count++;
    115         param= new Param(count,"velocity_observed",DOUBLEVEC);
     115        param= new Param(count,"u_g_obs",DOUBLEVEC);
    116116        param->SetDoubleVec(u_g_obs,2*model->numberofnodes,2);
    117117        parameters->AddObject(param);
    … …  
    121121
    122122        count++;
    123         param= new Param(count,"parameter",DOUBLEVEC);
     123        param= new Param(count,"param_g",DOUBLEVEC);
    124124        param->SetDoubleVec(p_g,model->numberofnodes,1);
    125125        parameters->AddObject(param);

  • issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp

    r586 r765  
    3939                        CreateLoadsDiagnosticHoriz(ploads,model,model_handle);
    4040                        CreateParametersDiagnosticHoriz(pparameters,model,model_handle);
     41                        if(model->qmu_analysis)CreateParametersQmu(pparameters,model,model_handle);
    4142                               
    4243                }
    … …  
    7778                CreateLoadsThermal(ploads,model,model_handle);
    7879                CreateParametersThermal(pparameters,model,model_handle);
     80                if(model->qmu_analysis)CreateParametersQmu(pparameters,model,model_handle);
    7981                                       
    8082        }
    … …  
    9294                CreateLoadsPrognostic(ploads,model,model_handle);
    9395                CreateParametersPrognostic(pparameters,model,model_handle);
     96                if(model->qmu_analysis)CreateParametersQmu(pparameters,model,model_handle);
    9497                                       
    9598        }

  • issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp

    r643 r765  
    2525        int      numberofdofspernode;
    2626        int      dim;
    27         int*  epart=NULL;
    28         int*  part=NULL;
    29         double*  dpart=NULL;
    30         int      elements_width;
    3127
    3228
    … …  
    4844        parameters->AddObject(param);
    4945
    50         //qmu analysis?
    51         count++;
    52         param= new Param(count,"qmu_analysis",INTEGER);
    53         param->SetInteger(model->qmu_analysis);
    54         parameters->AddObject(param);
    55 
    56         count++;
    57         param= new Param(count,"qmu_npart",INTEGER);
    58         param->SetInteger(model->qmu_npart);
    59         parameters->AddObject(param);
    60 
    6146        //dimension 2d or 3d:
    6247        if (strcmp(model->meshtype,"2d")==0)dim=2;
    … …  
    217202        param->SetInteger(numberofdofspernode);
    218203        parameters->AddObject(param);
    219 
    220         /*Deal with responses and partition for qmu modeling: */
    221         if(model->qmu_analysis){
    222                 char** descriptors=NULL;
    223                 char*  descriptor=NULL;
    224                 char*  tag=NULL;
    225 
    226                 descriptors=(char**)xmalloc(model->numberofresponses*sizeof(char*));
    227                 tag=(char*)xmalloc((strlen("descriptori")+1)*sizeof(char));
    228 
    229                 /*Fetch descriptors: */
    230                 for(i=0;i<model->numberofresponses;i++){
    231                         sprintf(tag,"%s%i","descriptor",i);
    232                         ModelFetchData((void**)&descriptor,NULL,NULL,model_handle,tag,"String",NULL);
    233                         descriptors[i]=descriptor;
    234                 }
    235 
    236                 /*Ok, we have all the descriptors. Build a parameter with it: */
    237                 count++;
    238                 param= new Param(count,"descriptors",STRINGARRAY);
    239                 param->SetStringArray(descriptors,model->numberofresponses);
    240                 parameters->AddObject(param);
    241 
    242                 /*Free data: */
    243                 xfree((void**)&tag);
    244                 for(i=0;i<model->numberofresponses;i++){
    245                         char* descriptor=descriptors[i];
    246                         xfree((void**)&descriptor);
    247                 }
    248                 xfree((void**)&descriptors);
    249 
    250                 #ifdef _PARALLEL_
    251                 /*partition grids in model->qmu_npart parts: */
    252        
    253                 if(strcmp(model->meshtype,"2d")==0){
    254                         ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
    255                         elements_width=3; //tria elements
    256                 }
    257                 else{
    258                         ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
    259                         elements_width=6; //penta elements
    260                 }
    261 
    262                 MeshPartitionx(&epart, &part,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,model->qmu_npart);
    263 
    264                 dpart=(double*)xmalloc(model->numberofnodes*sizeof(double));
    265                 for(i=0;i<model->numberofnodes;i++)dpart[i]=part[i];
    266 
    267                 count++;
    268                 param= new Param(count,"qmu_part",DOUBLEVEC);
    269                 param->SetDoubleVec(dpart,model->numberofnodes);
    270                 parameters->AddObject(param);
    271 
    272                 /*Free elements and elements2d: */
    273                 xfree((void**)&model->elements);
    274                 xfree((void**)&model->elements2d);
    275                 xfree((void**)&epart);
    276                 xfree((void**)&part);
    277                 xfree((void**)&dpart);
    278 
    279                 #endif
    280         }
    281 
    282204       
    283205        /*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp

    r465 r765  
    2323        double* vy=NULL;
    2424        double* vz=NULL;
     25        double* ug=NULL;
    2526
    2627        /*recover parameters : */
    … …  
    3738        ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat");
    3839
    39         if(vx) for(i=0;i<model->numberofnodes;i++)vx[i]=vx[i]/model->yts;
    40         if(vy) for(i=0;i<model->numberofnodes;i++)vy[i]=vy[i]/model->yts;
    41         if(vz) for(i=0;i<model->numberofnodes;i++)vz[i]=vz[i]/model->yts;
     40        ug=(double*)xcalloc(model->numberofnodes*3,sizeof(double));
     41
     42        if(vx) for(i=0;i<model->numberofnodes;i++)ug[3*i+0]=vx[i]/model->yts;
     43        if(vy) for(i=0;i<model->numberofnodes;i++)ug[3*i+1]=vy[i]/model->yts;
     44        if(vz) for(i=0;i<model->numberofnodes;i++)ug[3*i+2]=vz[i]/model->yts;
    4245
    4346        count++;
    44         param= new Param(count,"vx",DOUBLEVEC);
    45         if(vx) param->SetDoubleVec(vx,model->numberofnodes);
    46         else param->SetDoubleVec(vx,0);
    47         parameters->AddObject(param);
    48 
    49         count++;
    50         param= new Param(count,"vy",DOUBLEVEC);
    51         if(vy) param->SetDoubleVec(vy,model->numberofnodes);
    52         else param->SetDoubleVec(vy,0);
    53         parameters->AddObject(param);
    54 
    55         count++;
    56         param= new Param(count,"vz",DOUBLEVEC);
    57         if(vz) param->SetDoubleVec(vz,model->numberofnodes);
    58         else param->SetDoubleVec(vz,0);
     47        param= new Param(count,"u_g",DOUBLEVEC);
     48        param->SetDoubleVec(ug,3*model->numberofnodes,3);
    5949        parameters->AddObject(param);
    6050
    … …  
    6252        xfree((void**)&vy);
    6353        xfree((void**)&vz);
     54        xfree((void**)&ug);
    6455       
    6556        /*Assign output pointer: */

  • issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp

    r713 r765  
    4040
    4141                count++;
    42                 param= new Param(count,"melting",DOUBLEVEC);
     42                param= new Param(count,"m_g",DOUBLEVEC);
    4343                if(melting) param->SetDoubleVec(melting,model->numberofnodes);
    4444                else param->SetDoubleVec(melting,0);

  • issm/trunk/src/c/ModelProcessorx/Model.cpp

    r653 r765  
    3939        model->solverstring=NULL;
    4040        model->numberofresponses=0;
     41        model->numberofvariables=0;
    4142        model->qmu_npart=0;
    4243        model->numberofelements=0;
    … …  
    384385        /*qmu: */
    385386        if(model->qmu_analysis){
     387                ModelFetchData((void**)&model->numberofvariables,NULL,NULL,model_handle,"numberofvariables","Integer",NULL);
    386388                ModelFetchData((void**)&model->numberofresponses,NULL,NULL,model_handle,"numberofresponses","Integer",NULL);
    387389                ModelFetchData((void**)&model->qmu_npart,NULL,NULL,model_handle,"npart","Integer",NULL);

  • issm/trunk/src/c/ModelProcessorx/Model.h

    r653 r765  
    5858        /*qmu: */
    5959        int      numberofresponses;
     60        int      numberofvariables;
    6061        int      qmu_npart;
    6162
    … …  
    237238        void    CreateParametersPrognostic(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
    238239
     240        /*qmu: */
     241        void CreateParametersQmu(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
     242
    239243
    240244#endif  /* _MODEL_H */

  • issm/trunk/src/c/ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp

    r623 r765  
    2323        double* vy=NULL;
    2424        double* vz=NULL;
     25        double* u_g=NULL;
    2526        double* pressure=NULL;
    2627        double* thickness=NULL;
    … …  
    3839        ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat");
    3940
    40         if(vx) for(i=0;i<model->numberofnodes;i++)vx[i]=vx[i]/model->yts;
    41         if(vy) for(i=0;i<model->numberofnodes;i++)vy[i]=vy[i]/model->yts;
    42         if(vz) for(i=0;i<model->numberofnodes;i++)vz[i]=vz[i]/model->yts;
     41        u_g=(double*)xcalloc(model->numberofnodes*3,sizeof(double));
     42
     43        if(vx) for(i=0;i<model->numberofnodes;i++)u_g[3*i+0]=vx[i]/model->yts;
     44        if(vy) for(i=0;i<model->numberofnodes;i++)u_g[3*i+1]=vy[i]/model->yts;
     45        if(vz) for(i=0;i<model->numberofnodes;i++)u_g[3*i+2]=vz[i]/model->yts;
    4346
    4447        count++;
    45         param= new Param(count,"vx",DOUBLEVEC);
    46         if(vx) param->SetDoubleVec(vx,model->numberofnodes);
    47         else param->SetDoubleVec(vx,0);
     48        param= new Param(count,"u_g",DOUBLEVEC);
     49        param->SetDoubleVec(u_g,3*model->numberofnodes,3);
    4850        parameters->AddObject(param);
    4951
    50         count++;
    51         param= new Param(count,"vy",DOUBLEVEC);
    52         if(vy) param->SetDoubleVec(vy,model->numberofnodes);
    53         else param->SetDoubleVec(vy,0);
    54         parameters->AddObject(param);
    55 
    56         count++;
    57         param= new Param(count,"vz",DOUBLEVEC);
    58         if(vz) param->SetDoubleVec(vz,model->numberofnodes);
    59         else param->SetDoubleVec(vz,0);
    60         parameters->AddObject(param);
    6152
    6253        xfree((void**)&vx);
    6354        xfree((void**)&vy);
    6455        xfree((void**)&vz);
     56        xfree((void**)&u_g);
    6557       
    6658        /*Get thickness: */
    … …  
    6860       
    6961        count++;
    70         param= new Param(count,"thickness",DOUBLEVEC);
    71         if(thickness) param->SetDoubleVec(thickness,model->numberofnodes);
    72         else param->SetDoubleVec(thickness,0);
     62        param= new Param(count,"h_g",DOUBLEVEC);
     63        if(thickness) param->SetDoubleVec(thickness,model->numberofnodes,1);
     64        else param->SetDoubleVec(thickness,0,0);
    7365        parameters->AddObject(param);
    7466
    … …  
    8072       
    8173        count++;
    82         param= new Param(count,"melting",DOUBLEVEC);
    83         if(melting) param->SetDoubleVec(melting,model->numberofnodes);
    84         else param->SetDoubleVec(melting,0);
     74        param= new Param(count,"m_g",DOUBLEVEC);
     75        if(melting) param->SetDoubleVec(melting,model->numberofnodes,1);
     76        else param->SetDoubleVec(melting,0,1);
    8577        parameters->AddObject(param);
    8678
    … …  
    9284       
    9385        count++;
    94         param= new Param(count,"accumulation",DOUBLEVEC);
    95         if(accumulation) param->SetDoubleVec(accumulation,model->numberofnodes);
    96         else param->SetDoubleVec(accumulation,0);
     86        param= new Param(count,"a_g",DOUBLEVEC);
     87        if(accumulation) param->SetDoubleVec(accumulation,model->numberofnodes,1);
     88        else param->SetDoubleVec(accumulation,0,0);
    9789        parameters->AddObject(param);
    9890

  • issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp

    r713 r765  
    2323        double* vy=NULL;
    2424        double* vz=NULL;
     25        double* u_g=NULL;
    2526        double* pressure=NULL;
    2627        double* temperature=NULL;
    … …  
    3637        ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat");
    3738
    38         if(vx) for(i=0;i<model->numberofnodes;i++)vx[i]=vx[i]/model->yts;
    39         if(vy) for(i=0;i<model->numberofnodes;i++)vy[i]=vy[i]/model->yts;
    40         if(vz) for(i=0;i<model->numberofnodes;i++)vz[i]=vz[i]/model->yts;
     39        u_g=(double*)xcalloc(model->numberofnodes*3,sizeof(double));
     40
     41        if(vx) for(i=0;i<model->numberofnodes;i++)u_g[3*i+0]=vx[i]/model->yts;
     42        if(vy) for(i=0;i<model->numberofnodes;i++)u_g[3*i+1]=vy[i]/model->yts;
     43        if(vz) for(i=0;i<model->numberofnodes;i++)u_g[3*i+2]=vz[i]/model->yts;
    4144
    4245        count++;
    43         param= new Param(count,"vx",DOUBLEVEC);
    44         if(vx) param->SetDoubleVec(vx,model->numberofnodes);
    45         else param->SetDoubleVec(vx,0);
    46         parameters->AddObject(param);
    47 
    48         count++;
    49         param= new Param(count,"vy",DOUBLEVEC);
    50         if(vy) param->SetDoubleVec(vy,model->numberofnodes);
    51         else param->SetDoubleVec(vy,0);
    52         parameters->AddObject(param);
    53 
    54         count++;
    55         param= new Param(count,"vz",DOUBLEVEC);
    56         if(vz) param->SetDoubleVec(vz,model->numberofnodes);
    57         else param->SetDoubleVec(vz,0);
     46        param= new Param(count,"u_g",DOUBLEVEC);
     47        param->SetDoubleVec(u_g,3*model->numberofnodes,3);
    5848        parameters->AddObject(param);
    5949
    … …  
    6151        xfree((void**)&vy);
    6252        xfree((void**)&vz);
     53        xfree((void**)&u_g);
    6354       
    6455        /*Get pressure: */
    … …  
    6657       
    6758        count++;
    68         param= new Param(count,"pressure",DOUBLEVEC);
    69         if(pressure) param->SetDoubleVec(pressure,model->numberofnodes);
    70         else param->SetDoubleVec(pressure,0);
     59        param= new Param(count,"p_g",DOUBLEVEC);
     60        if(pressure) param->SetDoubleVec(pressure,model->numberofnodes,1);
     61        else param->SetDoubleVec(pressure,0,0);
    7162        parameters->AddObject(param);
    7263
    … …  
    7465        xfree((void**)&pressure);
    7566
    76         /* get initial temperature if transient*/
     67        /* get initial temperature and melting if transient*/
    7768        if(strcmp(model->sub_analysis_type,"transient")==0){
    7869
    … …  
    8172
    8273                count++;
    83                 param= new Param(count,"temperature",DOUBLEVEC);
     74                param= new Param(count,"t_g",DOUBLEVEC);
    8475                if(temperature) param->SetDoubleVec(temperature,model->numberofnodes);
    8576                else throw ErrorException(__FUNCT__,exprintf("Missing initial temperature"));

  • issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp

    r758 r765  
    11/*!\file ProcessParamsx
    2  * \brief: process parameters using partitioning vector
     2 * \brief: process parameters using partitioning vector.
     3 * Go through all parameters in the 'parameters' dataset. For each parameter that holds a doublevec (ie, a double* vector synchronized across
     4 * the MPI ring of a cluster), partition the vector so that the new node partitioning decided in ModelProcessor is applied. Otherwise,
     5 * parameters coming directly from Matlab would be serially partitioned, which means could not be recognized by each individual node or element.
     6 * The partition needs to be the parallel partitionting.
    37 */
    48
    … …  
    1115#include "../include/macros.h"
    1216#include "../toolkits/toolkits.h"
     17#include "../DataSet/DataSet.h"
    1318#include "../EnumDefinitions/EnumDefinitions.h"
    1419
    1520void ProcessParamsx( DataSet* parameters, Vec  part){
    1621
    17         int i;
    1822       
    1923        double* partition=NULL;
    2024        Param*  param=NULL;
    2125
    22         /*diagnostic: */
    23         double* vx=NULL;
    24         double* vy=NULL;
    25         double* vz=NULL;
    26                
    27         double* u_g=NULL;
    28        
    29         /*control: */
    30         Vec     vec_vx_obs=NULL;
    31         Vec     vec_vy_obs=NULL;
    32         Vec     vec_control_parameter=NULL;
    33 
    34         double* vx_obs=NULL;
    35         double* vy_obs=NULL;
    36         double* control_parameter=NULL;
    37         double* param_g=NULL;
    38 
    39         double* u_g_obs=NULL;
    40 
    41         /*thermal*/
    42         double* pressure=NULL;
    43         double* temperature=NULL;
    44         double* melting=NULL;
    45        
    46         double* p_g=NULL;
    47         double* t_g=NULL;
    48         double* m_g=NULL;
    49 
    50         /*prognostic: */
    51         double* a_g=NULL;
    52         double* h_g=NULL;
    53         double* accumulation=NULL;
    54         double* thickness=NULL;
    55 
    5626        /*diverse: */
    5727        int     numberofnodes;
    58         int     analysis_type;
    59         int     sub_analysis_type;
    60         int     count;
     28        int i,j,k;
     29        int ndof;
    6130
    62         parameters->FindParam((void*)&analysis_type,"analysis_type");
    63         parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
    64         count=parameters->Size();
    65                
     31        /*processing of parameters: */
     32        double* parameter=NULL;
     33        double* new_parameter=NULL;
     34
     35        /*Some parameters needed: */
    6636        parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    6737
    68         /*First serialize partition vector: */
     38        /*serialize partition vector: */
    6939        if(part)VecToMPISerial(&partition,part);
    70        
    71         if (   (analysis_type==ControlAnalysisEnum()) || 
    72                    (analysis_type==DiagnosticAnalysisEnum()) ||
    73                    (analysis_type==ThermalAnalysisEnum()) ||
    74                    (analysis_type==PrognosticAnalysisEnum())
    75                    ){
    7640
    77                 parameters->FindParam((void*)&vx,"vx");
    78                 parameters->FindParam((void*)&vy,"vy");
    79                 parameters->FindParam((void*)&vz,"vz");
     41        for(i=0;i<parameters->Size();i++){
    8042
    81                 if(numberofnodes){
    82                         u_g=(double*)xcalloc(numberofnodes*3,sizeof(double));
     43                Param* param=(Param*)parameters->GetObjectByOffset(i);
    8344
    84                         if(vx){
    85                                 for(i=0;i<numberofnodes;i++){
    86                                         u_g[(int)(3*partition[i]+0)]=vx[i]; 
     45                if (param->GetType()==DOUBLEVEC){
     46
     47                        ndof=param->GetNdof();
     48
     49                        if (ndof!=0){ /*ok, we are dealing with a parameter that needs to be repartitioned, for ndof degrees of freedom: */
     50
     51                                new_parameter=(double*)xcalloc(numberofnodes*ndof,sizeof(double));
     52                                param->GetParameterValue(&parameter);
     53
     54                                for(j=0;j<ndof;j++){
     55
     56                                        for(k=0;k<numberofnodes;k++){
     57                                       
     58                                                new_parameter[(int)(ndof*partition[k]+j)]=parameter[ndof*k+j];
     59
     60                                        }
     61
    8762                                }
     63                       
     64                               
     65                                /*Now, replace old parameter with new parameter: */
     66                                param->SetDoubleVec(new_parameter,ndof*numberofnodes,ndof);
     67
     68                                /*Free ressources: */
     69                                xfree((void**)&new_parameter);
     70
    8871                        }
    89                         if(vy){
    90                                 for(i=0;i<numberofnodes;i++){
    91                                         u_g[(int)(3*partition[i]+1)]=vy[i]; 
    92                                 }
    93                         }
    94                         if(vz){
    95                                 for(i=0;i<numberofnodes;i++){
    96                                         u_g[(int)(3*partition[i]+2)]=vz[i]; 
    97                                 }
    98                         }
    99 
    100                         /*Now, create new parameters: */
    101                         count++;
    102                         param= new Param(count,"u_g",DOUBLEVEC);
    103                         param->SetDoubleVec(u_g,3*numberofnodes);
    104                         parameters->AddObject(param);
    105 
    106                         /*erase old parameters: */
    107                         param=(Param*)parameters->FindParamObject("vx");
    108                         parameters->DeleteObject((Object*)param);
    109 
    110                         param=(Param*)parameters->FindParamObject("vy");
    111                         parameters->DeleteObject((Object*)param);
    112 
    113                         param=(Param*)parameters->FindParamObject("vz");
    114                         parameters->DeleteObject((Object*)param);
    11572                }
    11673        }
    11774
    118         if(analysis_type==ControlAnalysisEnum()){
    119 
    120                 parameters->FindParam((void*)&vx_obs,"vx_obs");
    121                 parameters->FindParam((void*)&vy_obs,"vy_obs");
    122                 parameters->FindParam((void*)&control_parameter,"control_parameter");
    123 
    124                 /*Now, from vx_obs and vy_obs, build u_g_obs, correctly partitioned: */
    125 
    126                 u_g_obs=(double*)xcalloc(numberofnodes*2,sizeof(double));
    127                 param_g=(double*)xcalloc(numberofnodes*2,sizeof(double));
    128 
    129                 for(i=0;i<numberofnodes;i++){
    130                         u_g_obs[(int)(2*partition[i]+0)]=vx_obs[i]; 
    131                         u_g_obs[(int)(2*partition[i]+1)]=vy_obs[i]; 
    132                         param_g[(int)(2*partition[i]+0)]=control_parameter[i];
    133                 }
    134 
    135                 /*Now, create new parameters: */
    136                 count++;
    137                 param= new Param(count,"u_g_obs",DOUBLEVEC);
    138                 param->SetDoubleVec(u_g_obs,2*numberofnodes);
    139                 parameters->AddObject(param);
    140 
    141                 count++;
    142                 param= new Param(count,"param_g",DOUBLEVEC);
    143                 param->SetDoubleVec(param_g,2*numberofnodes);
    144                 parameters->AddObject(param);
    145 
    146                 /*erase old parameter: */
    147                 param=(Param*)parameters->FindParamObject("vx_obs");
    148                 parameters->DeleteObject((Object*)param);
    149 
    150                 param=(Param*)parameters->FindParamObject("vy_obs");
    151                 parameters->DeleteObject((Object*)param);
    152 
    153                 param=(Param*)parameters->FindParamObject("control_parameter");
    154                 parameters->DeleteObject((Object*)param);
    155         }
    156 
    157         if(analysis_type==ThermalAnalysisEnum()){
    158 
    159                 parameters->FindParam((void*)&pressure,"pressure");
    160 
    161                 /*Now, from temperature and melting, build t_g and m_g, correctly partitioned: */
    162                 p_g=(double*)xcalloc(numberofnodes,sizeof(double));
    163 
    164                 for(i=0;i<numberofnodes;i++){
    165                         p_g[(int)(partition[i])]= pressure[i]; 
    166                 }
    167 
    168                 /*Now, create new parameter: */
    169                 count++;
    170                 param= new Param(count,"p_g",DOUBLEVEC);
    171                 param->SetDoubleVec(p_g,numberofnodes);
    172                 parameters->AddObject(param);
    173 
    174                 /*erase old parameter: */
    175                 param=(Param*)parameters->FindParamObject("pressure");
    176                 parameters->DeleteObject((Object*)param);
    177 
    178                 if (sub_analysis_type==TransientAnalysisEnum()){
    179 
    180                         parameters->FindParam((void*)&temperature,"temperature");
    181 
    182                         /*Now, from temperature , build t_g , correctly partitioned: */
    183                         t_g=(double*)xcalloc(numberofnodes,sizeof(double));
    184 
    185                         for(i=0;i<numberofnodes;i++){
    186                                 t_g[(int)(partition[i])]=temperature[i];
    187                         }
    188 
    189                         /*Now, create new parameter: */
    190                         count++;
    191                         param= new Param(count,"t_g",DOUBLEVEC);
    192                         param->SetDoubleVec(t_g,numberofnodes);
    193                         parameters->AddObject(param);
    194 
    195                         /*erase old parameter: */
    196                         param=(Param*)parameters->FindParamObject("temperature");
    197                         parameters->DeleteObject((Object*)param);
    198                 }
    199         }
    200 
    201         if(analysis_type==MeltingAnalysisEnum() && sub_analysis_type==TransientAnalysisEnum()){
    202 
    203                 parameters->FindParam((void*)&melting,"melting");
    204 
    205                 /*Now, from melting , build m_g , correctly partitioned: */
    206                 m_g=(double*)xcalloc(numberofnodes,sizeof(double));
    207 
    208                 for(i=0;i<numberofnodes;i++){
    209                         m_g[(int)(partition[i])]=melting[i];
    210                 }
    211 
    212                 /*Now, create new parameter: */
    213                 count++;
    214                 param= new Param(count,"m_g",DOUBLEVEC);
    215                 param->SetDoubleVec(m_g,numberofnodes);
    216                 parameters->AddObject(param);
    217 
    218                 /*erase old parameter: */
    219                 param=(Param*)parameters->FindParamObject("melting");
    220                 parameters->DeleteObject((Object*)param);
    221         }
    222 
    223         if(analysis_type==PrognosticAnalysisEnum()){
    224 
    225                 parameters->FindParam((void*)&melting,"melting");
    226                 parameters->FindParam((void*)&thickness,"thickness");
    227                 parameters->FindParam((void*)&accumulation,"accumulation");
    228                
    229                 /*Now, from melting accumulation and thickness, build m_g, h_g, correctly partitioned: */
    230                 m_g=(double*)xcalloc(numberofnodes,sizeof(double));
    231                 h_g=(double*)xcalloc(numberofnodes,sizeof(double));
    232                 a_g=(double*)xcalloc(numberofnodes,sizeof(double));
    233 
    234                 for(i=0;i<numberofnodes;i++){
    235                         m_g[(int)(partition[i])]= melting[i]; 
    236                         h_g[(int)(partition[i])]= thickness[i]; 
    237                         a_g[(int)(partition[i])]= accumulation[i]; 
    238                 }
    239 
    240                 /*Now, create new parameter: */
    241                 count++;
    242                 param= new Param(count,"m_g",DOUBLEVEC);
    243                 param->SetDoubleVec(m_g,numberofnodes);
    244                 parameters->AddObject(param);
    245 
    246                 count++;
    247                 param= new Param(count,"h_g",DOUBLEVEC);
    248                 param->SetDoubleVec(h_g,numberofnodes);
    249                 parameters->AddObject(param);
    250 
    251                 count++;
    252                 param= new Param(count,"a_g",DOUBLEVEC);
    253                 param->SetDoubleVec(a_g,numberofnodes);
    254                 parameters->AddObject(param);
    255 
    256                 /*erase old parameter: */
    257                 param=(Param*)parameters->FindParamObject("melting");
    258                 parameters->DeleteObject((Object*)param);
    259                
    260                 param=(Param*)parameters->FindParamObject("thickness");
    261                 parameters->DeleteObject((Object*)param);
    262 
    263                 param=(Param*)parameters->FindParamObject("accumulation");
    264                 parameters->DeleteObject((Object*)param);
    265         }
    266         xfree((void**)&partition);
    267        
    268         xfree((void**)&vx);
    269         xfree((void**)&vy);
    270         xfree((void**)&vz);
    271         xfree((void**)&u_g);
    272 
    273         xfree((void**)&vx_obs);
    274         xfree((void**)&vy_obs);
    275         xfree((void**)&param_g);
    276         xfree((void**)&pressure);
    277         xfree((void**)&control_parameter);
    278         xfree((void**)&u_g_obs);
    279         xfree((void**)&p_g);
    280         xfree((void**)&a_g);
    281         xfree((void**)&h_g);
    282         xfree((void**)&m_g);
    28375}

  • issm/trunk/src/c/objects/Param.cpp

    r586 r765  
    125125               
    126126                case DOUBLEVEC:
    127                         printf("   double vector. size: %i\n",M);
     127                        printf("   double vector. size: %i ndof: %i\n",M,ndof);
    128128                        for(i=0;i<M;i++)printf("%g\n",doublevec[i]);
    129129                        break;
    … …  
    178178        memcpy(marshalled_dataset,&name,sizeof(name));marshalled_dataset+=sizeof(name);
    179179        memcpy(marshalled_dataset,&type,sizeof(type));marshalled_dataset+=sizeof(type);
     180        memcpy(marshalled_dataset,&ndof,sizeof(ndof));marshalled_dataset+=sizeof(ndof);
    180181
    181182        switch(type){
    … …  
    265266                sizeof(name)+
    266267                sizeof(type)+
     268                sizeof(ndof)+
    267269                sizeof(int); //sizeof(int) for enum type
    268270
    … …  
    333335        memcpy(&name,marshalled_dataset,sizeof(name));marshalled_dataset+=sizeof(name);
    334336        memcpy(&type,marshalled_dataset,sizeof(type));marshalled_dataset+=sizeof(type);
     337        memcpy(&ndof,marshalled_dataset,sizeof(ndof));marshalled_dataset+=sizeof(ndof);
    335338
    336339
    … …  
    575578}
    576579
     580int   Param::GetNdof(){
     581        return ndof;
     582}
     583
    577584int   Param::GetN(){
    578585        return N;
    … …  
    589596                memcpy(doublevec,value,M*sizeof(double));
    590597        }
    591 
    592 }
     598        ndof=0; //this vector will not be repartitioned.
     599
     600}
     601
     602#undef __FUNCT__
     603#define __FUNCT__ "SetDoubleVec"
     604void  Param::SetDoubleVec(double* value,int size, int numberofdofs){
     605        if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublevecfor type",type));
     606       
     607        M=size;
     608        if(M){
     609                doublevec=(double*)xmalloc(M*sizeof(double));
     610                memcpy(doublevec,value,M*sizeof(double));
     611        }
     612        ndof=numberofdofs;
     613}
     614
    593615
    594616#undef __FUNCT__

  • issm/trunk/src/c/objects/Param.h

    r586 r765  
    3030                int M; //dimensions of vectors and matrices
    3131                int N;
     32                int ndof;
    3233
    3334        public:
    … …  
    4748                void  SetDouble(double value);
    4849                void  SetDoubleVec(double* value,int size);
     50                void  SetDoubleVec(double* value,int size,int ndof);
    4951                void  SetVec(Vec value);
    5052                void  SetInteger(int value);
    … …  
    5961                int   GetM();
    6062                int   GetN();
     63                int   GetNdof();
    6164                Object* copy();
    6265

  • issm/trunk/src/c/objects/ParameterInputs.cpp

    r365 r765  
    270270        return ug;
    271271}
     272
     273#undef __FUNCT__
     274#define __FUNCT__ "ParameterInputs::UpdateFromDakota"
     275void ParameterInputs::UpdateFromDakota(double* variables,char** variables_descriptors,int numvariables,DataSet* parameters,double* qmu_part,int qmu_npart){
     276
     277        /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and
     278         * for each descriptor, take the variable value and plug it into the inputs: */
     279
     280        int     i,j,k;
     281        int     found;
     282
     283        double* distributed_values=NULL;
     284        char*   root=NULL; //root name of distributed variable, ex: thickness, drag, etc ...
     285        double* parameter=NULL;
     286        int     numberofnodes;
     287
     288        char*   descriptor=NULL;
     289
     290        /*retrive some necessary parameters first: */
     291        found=parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     292        if(!found)throw ErrorException(__FUNCT__,"coud not find numberofnodes in parameters dataset!");
     293
     294        for(i=0;i<numvariables;i++){
     295       
     296                descriptor=variables_descriptors[i];
     297
     298                /*From descriptor, figure out if the variable is distributed (distributed implies there is a numeric value at the
     299                 * end of the descriptor, for ex: thickness1, thickness10, etc .... If it is distributed, the next qmu_npart (number
     300                 * of partitions in the distributed variable) variable are the values for each partition of the distributed variable: */
     301                if (!isdistributed(&root,descriptor)){
     302                       
     303                        /*Ok, variable is not distributed, just add it to inputs: */
     304                        this->Add(descriptor,variables[i]);
     305                }
     306                else{
     307                       
     308                        /*Ok, variable is distributed. Root name of variable is also known. Now, allocate distributed_values and fill the
     309                         * distributed_values with the next qmu_npart variables: */
     310
     311                        distributed_values=(double*)xmalloc(qmu_npart*sizeof(double));
     312                        for(j=0;j<qmu_npart;j++){
     313                                distributed_values[j]=variables[i+j];
     314                        }
     315
     316                        /*Now, pick up the parameter corresponding to root: */
     317                        found=parameters->FindParam((void*)&parameter,root);
     318                        if(!found)throw ErrorException(__FUNCT__,exprintf("%s%s%s"," could not find parameter ",root," for Dakota input update"));
     319                       
     320                        /*We've got the parameter, we need to update it using qmu_part (a partitioning vector), and the distributed_values: */
     321                        for(k=0;k<numberofnodes;k++){
     322                                parameter[k]=parameter[k]*distributed_values[(int)qmu_part[k]];
     323                        }
     324
     325                        /*Add parameter to inputs: */
     326                        this->Add("root",parameter,1,numberofnodes);
     327
     328                        /*increment i to skip the distributed values just collected: */
     329                        i+=qmu_npart;
     330
     331                        /*Free allocations: */
     332                        xfree((void**)&parameter);
     333                        xfree((void**)&distributed_values);
     334                        xfree((void**)&root);
     335
     336                }
     337        }
     338
     339        /*Free ressources:*/
     340        xfree((void**)&distributed_values);
     341        xfree((void**)&root);
     342        xfree((void**)&parameter);
     343
     344}

  • issm/trunk/src/c/objects/ParameterInputs.h

    r365 r765  
    3838                void Init( void* data_handle);
    3939                void purge(char* name);
     40                void UpdateFromDakota(double* variables,char** variables_descriptors,int numvariables,DataSet* parameters, double* qmu_part,int qmu_npart);
     41
    4042};
    4143

  • issm/trunk/src/c/parallel/SpawnCore.cpp

    r659 r765  
    3636        extern int my_rank;
    3737       
     38       
    3839        /*synchronize all cpus, as CPU 0 is probably late (it is starting the entire dakota strategy!) : */
    3940        MPI_Barrier(MPI_COMM_WORLD);
    4041
    41         /*First off, recover the responses descriptors for the response functions: */
    42         param=(Param*)femmodels[0].parameters->FindParamObject("descriptors");
     42        /*First off, recover the response descriptors for the response functions: */
     43        param=(Param*)femmodels[0].parameters->FindParamObject("responsedescriptors");
     44        if(!param)throw ErrorException(__FUNCT__," could not find response descriptors!");
     45
    4346        numresponses=param->GetM();
    4447        param->GetParameterValue(&responses_descriptors);
    … …  
    4750        femmodels[0].parameters->FindParam((void*)&qmu_npart,"qmu_npart");
    4851        femmodels[0].parameters->FindParam((void*)&qmu_part,"qmu_part");
    49 
    5052        #ifdef _DEBUG_
    5153        for(i=0;i<numresponses;i++){
    52                 PetscSynchronizedPrintf(MPI_COMM_WORLD,"descriptor %i: %s\n",i,responses_descriptors[i]);
     54                PetscSynchronizedPrintf(MPI_COMM_WORLD,"response descriptor %i: %s\n",i,responses_descriptors[i]);
    5355                PetscSynchronizedFlush(MPI_COMM_WORLD);
    5456        }
    … …  
    8991        #endif
    9092
     93
    9194        _printf_("initialize results:\n");
    9295        results=new DataSet(ResultsEnum());
    9396
    9497        /*Modify core inputs to reflect the dakota variables inputs: */
    95         //inputs->UpdateFromDakota(variables,variables_descriptors,numvariables,qmu_part,qmu_npart);
     98        inputs->UpdateFromDakota(variables,variables_descriptors,numvariables,femmodels[0].parameters,qmu_part,qmu_npart); //femmodel #0 is the one holding the parameters for Dakota.
     99       
    96100
    97101        /*Run the analysis core solution sequence, with the updated inputs: */
    98102        if(analysis_type==DiagnosticAnalysisEnum()){
    99 
     103                       
    100104                _printf_("Starting diagnostic core\n");
    101105                diagnostic_core(results,femmodels,inputs);
    102106
    103107        }
    104         else{
    105                 throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type ",analysis_type," and sub_analysis_type ",sub_analysis_type," not supported yet!"));
     108        else if(analysis_type==ThermalAnalysisEnum()){
     109               
     110                _printf_("Starting thermal core\n");
     111                thermal_core(results,femmodels,inputs);
     112
    106113        }
     114        else if(analysis_type==PrognosticAnalysisEnum()){
     115
     116                _printf_("Starting prognostic core\n");
     117                prognostic_core(results,femmodels,inputs);
     118
     119        }
     120        else if(analysis_type==TransientAnalysisEnum()){
     121
     122                _printf_("Starting transient core\n");
     123                throw ErrorException(__FUNCT__,"not supported yet!");
     124
     125        }
     126        else throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type ",analysis_type," and sub_analysis_type ",sub_analysis_type," not supported yet!"));
     127       
     128               
    107129
    108130        /*Now process the outputs, before computing the dakota responses: */
    109131        _printf_("process results:\n");
    110132        ProcessResults(&results,femmodels,analysis_type);
     133       
    111134
    112135        /*compute responses on cpu 0: dummy for now! */

  • issm/trunk/src/c/parallel/prognostic.cpp

    r713 r765  
    3333        Vec     h_g=NULL;
    3434        Vec     u_g=NULL;
     35        Vec u_g=NULL;
    3536        double* u_g_serial=NULL;
    3637        double* h_g_initial=NULL;

  • issm/trunk/src/c/parallel/thermal.cpp

    r693 r765  
    9191        //erase velocities and pressure embedded in parameters
    9292        param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
     93        param=(Param*)femmodels[0].parameters->FindParamObject("velocity");
    9394        femmodels[0].parameters->DeleteObject((Object*)param);
    9495        param=(Param*)femmodels[0].parameters->FindParamObject("p_g");

  • issm/trunk/src/c/shared/String/isdistributed.cpp

    r679 r765  
    1212#endif
    1313
     14#include "stdio.h"
     15#include <string.h>
     16#include <ctype.h>
     17#include "../Alloc/alloc.h"
     18
    1419#undef __FUNCT__
    1520#define __FUNCT__ "isdistributed"
    16 
    1721int isdistributed(char** proot,char* name){
    1822

  • issm/trunk/src/m/classes/@model/model.m

    r764 r765  
    268268        md.qmu_analysis=0;
    269269        md.iresp=1;
     270        md.ivar=1;
    270271        md.npart=0;
    271272

  • issm/trunk/src/m/classes/public/displayqmu.m

    r261 r765  
    8585if isempty(md.dakotadat), disp(sprintf('         dakotadat: N/A')); else disp(sprintf('         dakotadat:    [%ix%i]    (can be accessed by typing md.dakotadat)',size(md.dakotadat)));end
    8686disp(sprintf('         npart   : %i (number of partitions for semi-descrete qmu)',md.npart));
    87 disp(sprintf('         numrifts: %i (number of rifts for semi-descrete qmu)',md.numrifts));

  • issm/trunk/src/m/classes/public/marshall.m

    r659 r765  
    159159WriteData(fid,md.qmu_analysis,'Integer','qmu_analysis');
    160160if md.qmu_analysis,
    161         responses=md.responses(md.iresp).rf;
    162         WriteData(fid,numel(responses),'Integer','numberofresponses');
    163         for i=1:numel(responses),
    164                 WriteData(fid,responses(i).descriptor,'String',['descriptor' num2str(i-1)]);
    165         end
     161        %recover variables for corresponding analysis
     162        variables=md.variables(md.ivar);
     163
     164        %count how many variables we have
     165        numvariables=0;
     166        variable_fieldnames=fieldnames(variables);
     167        for i=1:length(variable_fieldnames),
     168                field_name=variable_fieldnames{i};
     169                fieldvariables=variables.(field_name);
     170                numvariables=numvariables+numel(variables.(field_name));
     171        end
     172        %write number of variables to disk
     173        WriteData(fid,numvariables,'Integer','numberofvariables');
     174
     175        %now, for each variable, write descriptor
     176        count=0;
     177        for i=1:length(variable_fieldnames),
     178                field_name=variable_fieldnames{i};
     179                fieldvariables=variables.(field_name);
     180                for j=1:length(fieldvariables),
     181                        descriptor=fieldvariables(j).descriptor;
     182                        WriteData(fid,descriptor,'String',['variabledescriptor' num2str(count)]);
     183                        count=count+1;
     184                end
     185        end
     186       
     187        %deal with responses
     188
     189        %recover responses for corresponding analysis
     190        responses=md.responses(md.iresp);
     191
     192        %count how many responses we have
     193        numresponses=0;
     194        response_fieldnames=fieldnames(responses);
     195        for i=1:length(response_fieldnames),
     196                field_name=response_fieldnames{i};
     197                fieldresponses=responses.(field_name);
     198                numresponses=numresponses+numel(responses.(field_name));
     199        end
     200        %write number of responses to disk
     201        WriteData(fid,numresponses,'Integer','numberofresponses');
     202
     203        %now, for each response, write descriptor
     204        count=0;
     205        for i=1:length(response_fieldnames),
     206                field_name=response_fieldnames{i};
     207                fieldresponses=responses.(field_name);
     208                for j=1:length(fieldresponses),
     209                        descriptor=fieldresponses(j).descriptor;
     210                        WriteData(fid,descriptor,'String',['responsedescriptor' num2str(count)]);
     211                        count=count+1;
     212                end
     213        end
     214       
     215       
     216        %write npart to disk
    166217        WriteData(fid,md.npart,'Integer','npart');
    167218end

  • issm/trunk/src/m/classes/public/qmuparallel.m

    r586 r765  
    3737        return;
    3838end
     39
     40%reset qmu_analysis counter
     41md.qmu_analysis=0;

  • issm/trunk/src/m/solutions/dakota/qmuin.m

    r659 r765  
    8282cd ../
    8383
    84 %keep iresp for response function handling in marshall routine.
     84%keep ivar iresp for variable and response function handling in marshall routine.
    8585md.iresp=iresp;
     86md.ivar=ivar;
Note: See TracChangeset for help on using the changeset viewer.