Changeset 1907


Ignore:
Timestamp:
08/26/09 08:39:55 (15 years ago)
Author:
seroussi
Message:

removed control solution and put it in diagnostic

Location:
issm/trunk/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

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

    r1904 r1907  
    6161                        if ( *(image+i*samps+j) == 0 ){
    6262                                if ( (j > 3) && (j < samps-4) && (i > 3) && (i < lines-4)){
    63                                         for ( k = 0; k < 4; k++ ){
    64                                                 for ( l = 0; l < 4; l++ ){
     63                                        for ( k = 0; k < 30; k++ ){
     64                                                for ( l = 0; l < 30; l++ ){
    6565                                                        *(image2+samps*(i+k)+j+l)=0;
    6666                                                        *(image2+samps*(i-k)+j+l)=0;

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

    r1837 r1907  
    4242        count=parameters->Size();
    4343       
    44         /*control_type: */
     44        //control analysis?
    4545        count++;
    46         param= new Param(count,"control_type",STRING);
    47         param->SetString(iomodel->control_type);
    48         parameters->AddObject(param);
    49 
    50         /*extrude_param: */
    51         count++;
    52         param= new Param(count,"extrude_param",DOUBLE);
    53         if (strcmp(iomodel->control_type,"drag")==0)   param->SetDouble(0);
    54         else if (strcmp(iomodel->control_type,"B")==0) param->SetDouble(1);
    55         else throw ErrorException(__FUNCT__,exprintf("control_type %s not supported yet!",iomodel->control_type));
    56         parameters->AddObject(param);
    57 
    58         /*nsteps: */
    59         count++;
    60         param= new Param(count,"nsteps",INTEGER);
    61         param->SetInteger(iomodel->nsteps);
    62         parameters->AddObject(param);
    63 
    64         /*tolx: */
    65         count++;
    66         param= new Param(count,"tolx",DOUBLE);
    67         param->SetDouble(iomodel->tolx);
    68         parameters->AddObject(param);
    69 
    70         /*mincontrolconstraint: */
    71         count++;
    72         param= new Param(count,"mincontrolconstraint",DOUBLE);
    73         param->SetDouble(iomodel->mincontrolconstraint);
    74         parameters->AddObject(param);
    75 
    76         /*maxcontrolconstraint: */
    77         count++;
    78         param= new Param(count,"maxcontrolconstraint",DOUBLE);
    79         param->SetDouble(iomodel->maxcontrolconstraint);
     46        param= new Param(count,"control_analysis",INTEGER);
     47        param->SetInteger(iomodel->control_analysis);
    8048        parameters->AddObject(param);
    8149       
    82         /*epsvel: */
    83         count++;
    84         param= new Param(count,"epsvel",DOUBLE);
    85         param->SetDouble(iomodel->epsvel);
    86         parameters->AddObject(param);
    87        
    88         /*meanvel: */
    89         count++;
    90         param= new Param(count,"meanvel",DOUBLE);
    91         param->SetDouble(iomodel->meanvel);
    92         parameters->AddObject(param);
     50        if(iomodel->control_analysis){
     51                /*control_type: */
     52                count++;
     53                param= new Param(count,"control_type",STRING);
     54                param->SetString(iomodel->control_type);
     55                parameters->AddObject(param);
    9356
    94         /*Now, recover fit, optscal and maxiter as vectors: */
    95         IoModelFetchData((void**)&iomodel->fit,NULL,NULL,iomodel_handle,"fit","Matrix","Mat");
    96         IoModelFetchData((void**)&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal","Matrix","Mat");
    97         IoModelFetchData((void**)&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter","Matrix","Mat");
     57                /*extrude_param: */
     58                count++;
     59                param= new Param(count,"extrude_param",DOUBLE);
     60                if (strcmp(iomodel->control_type,"drag")==0)   param->SetDouble(0);
     61                else if (strcmp(iomodel->control_type,"B")==0) param->SetDouble(1);
     62                else throw ErrorException(__FUNCT__,exprintf("control_type %s not supported yet!",iomodel->control_type));
     63                parameters->AddObject(param);
    9864
    99         count++;
    100         param= new Param(count,"fit",DOUBLEVEC);
    101         param->SetDoubleVec(iomodel->fit,iomodel->nsteps);
    102         parameters->AddObject(param);
     65                /*nsteps: */
     66                count++;
     67                param= new Param(count,"nsteps",INTEGER);
     68                param->SetInteger(iomodel->nsteps);
     69                parameters->AddObject(param);
    10370
    104         count++;
    105         param= new Param(count,"optscal",DOUBLEVEC);
    106         param->SetDoubleVec(iomodel->optscal,iomodel->nsteps);
    107         parameters->AddObject(param);
     71                /*tolx: */
     72                count++;
     73                param= new Param(count,"tolx",DOUBLE);
     74                param->SetDouble(iomodel->tolx);
     75                parameters->AddObject(param);
    10876
    109         count++;
    110         param= new Param(count,"maxiter",DOUBLEVEC);
    111         param->SetDoubleVec(iomodel->maxiter,iomodel->nsteps);
    112         parameters->AddObject(param);
     77                /*mincontrolconstraint: */
     78                count++;
     79                param= new Param(count,"mincontrolconstraint",DOUBLE);
     80                param->SetDouble(iomodel->mincontrolconstraint);
     81                parameters->AddObject(param);
    11382
    114         xfree((void**)&iomodel->fit);
    115         xfree((void**)&iomodel->optscal);
    116         xfree((void**)&iomodel->maxiter);
     83                /*maxcontrolconstraint: */
     84                count++;
     85                param= new Param(count,"maxcontrolconstraint",DOUBLE);
     86                param->SetDouble(iomodel->maxcontrolconstraint);
     87                parameters->AddObject(param);
     88               
     89                /*epsvel: */
     90                count++;
     91                param= new Param(count,"epsvel",DOUBLE);
     92                param->SetDouble(iomodel->epsvel);
     93                parameters->AddObject(param);
     94               
     95                /*meanvel: */
     96                count++;
     97                param= new Param(count,"meanvel",DOUBLE);
     98                param->SetDouble(iomodel->meanvel);
     99                parameters->AddObject(param);
    117100
    118         /*Get vx, vx_obs, vy, vy_obs, and the parameter value: */
    119         IoModelFetchData((void**)&vx,NULL,NULL,iomodel_handle,"vx","Matrix","Mat");
    120         IoModelFetchData((void**)&vy,NULL,NULL,iomodel_handle,"vy","Matrix","Mat");
    121         IoModelFetchData((void**)&vz,NULL,NULL,iomodel_handle,"vz","Matrix","Mat");
    122         IoModelFetchData((void**)&vx_obs,NULL,NULL,iomodel_handle,"vx_obs","Matrix","Mat");
    123         IoModelFetchData((void**)&vy_obs,NULL,NULL,iomodel_handle,"vy_obs","Matrix","Mat");
    124         IoModelFetchData((void**)&control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type,"Matrix","Mat");
     101                /*Now, recover fit, optscal and maxiter as vectors: */
     102                IoModelFetchData((void**)&iomodel->fit,NULL,NULL,iomodel_handle,"fit","Matrix","Mat");
     103                IoModelFetchData((void**)&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal","Matrix","Mat");
     104                IoModelFetchData((void**)&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter","Matrix","Mat");
    125105
    126         u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
    127         if(vx)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
    128         if(vy)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
    129         if(vz)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
     106                count++;
     107                param= new Param(count,"fit",DOUBLEVEC);
     108                param->SetDoubleVec(iomodel->fit,iomodel->nsteps);
     109                parameters->AddObject(param);
    130110
    131         count++;
    132         param= new Param(count,"u_g",DOUBLEVEC);
    133         param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
    134         parameters->AddObject(param);
     111                count++;
     112                param= new Param(count,"optscal",DOUBLEVEC);
     113                param->SetDoubleVec(iomodel->optscal,iomodel->nsteps);
     114                parameters->AddObject(param);
    135115
    136         u_g_obs=(double*)xcalloc(iomodel->numberofnodes*2,sizeof(double));
    137         if(vx_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+0]=vx_obs[i]/iomodel->yts;
    138         if(vy_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+1]=vy_obs[i]/iomodel->yts;
     116                count++;
     117                param= new Param(count,"maxiter",DOUBLEVEC);
     118                param->SetDoubleVec(iomodel->maxiter,iomodel->nsteps);
     119                parameters->AddObject(param);
    139120
    140         count++;
    141         param= new Param(count,"u_g_obs",DOUBLEVEC);
    142         param->SetDoubleVec(u_g_obs,2*iomodel->numberofnodes,2);
    143         parameters->AddObject(param);
    144        
    145         param_g=(double*)xcalloc(iomodel->numberofnodes,sizeof(double));
    146         for(i=0;i<iomodel->numberofnodes;i++)param_g[i]=control_parameter[i];
     121                xfree((void**)&iomodel->fit);
     122                xfree((void**)&iomodel->optscal);
     123                xfree((void**)&iomodel->maxiter);
    147124
    148         count++;
    149         param= new Param(count,"param_g",DOUBLEVEC);
    150         param->SetDoubleVec(param_g,iomodel->numberofnodes,1);
    151         parameters->AddObject(param);
     125                /*Get vx, vx_obs, vy, vy_obs, and the parameter value: */
     126                IoModelFetchData((void**)&vx,NULL,NULL,iomodel_handle,"vx","Matrix","Mat");
     127                IoModelFetchData((void**)&vy,NULL,NULL,iomodel_handle,"vy","Matrix","Mat");
     128                IoModelFetchData((void**)&vz,NULL,NULL,iomodel_handle,"vz","Matrix","Mat");
     129                IoModelFetchData((void**)&vx_obs,NULL,NULL,iomodel_handle,"vx_obs","Matrix","Mat");
     130                IoModelFetchData((void**)&vy_obs,NULL,NULL,iomodel_handle,"vy_obs","Matrix","Mat");
     131                IoModelFetchData((void**)&control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type,"Matrix","Mat");
    152132
    153         xfree((void**)&vx);
    154         xfree((void**)&vy);
    155         xfree((void**)&vz);
    156         xfree((void**)&u_g);
    157         xfree((void**)&vx_obs);
    158         xfree((void**)&vy_obs);
    159         xfree((void**)&u_g_obs);
    160         xfree((void**)&param_g);
    161         xfree((void**)&control_parameter);
     133                u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
     134                if(vx)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
     135                if(vy)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
     136                if(vz)for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
     137
     138                count++;
     139                param= new Param(count,"u_g",DOUBLEVEC);
     140                param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
     141                parameters->AddObject(param);
     142
     143                u_g_obs=(double*)xcalloc(iomodel->numberofnodes*2,sizeof(double));
     144                if(vx_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+0]=vx_obs[i]/iomodel->yts;
     145                if(vy_obs)for(i=0;i<iomodel->numberofnodes;i++)u_g_obs[2*i+1]=vy_obs[i]/iomodel->yts;
     146
     147                count++;
     148                param= new Param(count,"u_g_obs",DOUBLEVEC);
     149                param->SetDoubleVec(u_g_obs,2*iomodel->numberofnodes,2);
     150                parameters->AddObject(param);
     151               
     152                param_g=(double*)xcalloc(iomodel->numberofnodes,sizeof(double));
     153                for(i=0;i<iomodel->numberofnodes;i++)param_g[i]=control_parameter[i];
     154
     155                count++;
     156                param= new Param(count,"param_g",DOUBLEVEC);
     157                param->SetDoubleVec(param_g,iomodel->numberofnodes,1);
     158                parameters->AddObject(param);
     159
     160                xfree((void**)&vx);
     161                xfree((void**)&vy);
     162                xfree((void**)&vz);
     163                xfree((void**)&u_g);
     164                xfree((void**)&vx_obs);
     165                xfree((void**)&vy_obs);
     166                xfree((void**)&u_g_obs);
     167                xfree((void**)&param_g);
     168                xfree((void**)&control_parameter);
     169        }
    162170
    163171        /*Assign output pointer: */

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

    r1838 r1907  
    2323        CreateParameters(pparameters,iomodel,iomodel_handle);
    2424        CreateParametersQmu(pparameters,iomodel,iomodel_handle);
     25        CreateParametersControl(pparameters,iomodel,iomodel_handle);
    2526
    2627        /*This is just a high level driver: */
    27         if (iomodel->analysis_type==ControlAnalysisEnum()){
    28 
    29                 if (iomodel->sub_analysis_type==HorizAnalysisEnum()){
    30 
    31                         CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
    32                         CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle);
    33                         CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle);
    34                         CreateParametersControl(pparameters,iomodel,iomodel_handle);
    35 
    36                 }
    37                 else if (iomodel->sub_analysis_type==VertAnalysisEnum()){
    38 
    39                         CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
    40                         CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle);
    41                         CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle);
    42 
    43                 }
    44                 else if (iomodel->sub_analysis_type==StokesAnalysisEnum()){
    45 
    46                         CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
    47                         CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle);
    48                         CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle);
    49                         CreateParametersControl(pparameters,iomodel,iomodel_handle);
    50 
    51                 }
    52                 else if (iomodel->sub_analysis_type==HutterAnalysisEnum()){
    53 
    54                         CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
    55                         CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle);
    56                         CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle);
    57                         CreateParametersControl(pparameters,iomodel,iomodel_handle);
    58 
    59                 }
    60 
    61         }
    62         else if (iomodel->analysis_type==DiagnosticAnalysisEnum()){
     28        if (iomodel->analysis_type==DiagnosticAnalysisEnum()){
    6329
    6430                if (iomodel->sub_analysis_type==HorizAnalysisEnum()){

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

    r1864 r1907  
    3838        iomodel->sub_analysis_type=0;
    3939        iomodel->qmu_analysis=0;
     40        iomodel->control_analysis=0;
    4041        iomodel->solverstring=NULL;
    4142        iomodel->numberofresponses=0;
    … …  
    286287        IoModelFetchData((void**)&iomodel->sub_analysis_type,NULL,NULL,iomodel_handle,"sub_analysis_type","Integer",NULL);
    287288        IoModelFetchData((void**)&iomodel->qmu_analysis,NULL,NULL,iomodel_handle,"qmu_analysis","Integer",NULL);
     289        IoModelFetchData((void**)&iomodel->control_analysis,NULL,NULL,iomodel_handle,"control_analysis","Integer",NULL);
    288290        IoModelFetchData((void**)&iomodel->meshtype,NULL,NULL,iomodel_handle,"type","String",NULL);
    289291        /*!Get numberofelements and numberofnodes: */

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

    r1881 r1907  
    2121        int     sub_analysis_type;
    2222        int     qmu_analysis;
     23        int     control_analysis;
    2324        char*   solverstring;
    2425

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

    r1904 r1907  
    11191119        if(inputs->Recover("temperature",&temperature_list[0],1,dofs,6,(void**)nodes)){
    11201120                if(matice && !collapse){
    1121                         B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
    1122                                                 +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
     1121                        //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
     1122                        //                      +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
     1123                        temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0;
     1124                        B_average=Paterson(temperature_average);
    11231125                        matice->SetB(B_average);
    11241126                }
    … …  
    11271129        if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){
    11281130                if(matice && collapse){
    1129                         B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
    1130                                                 +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
     1131                        temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0;
     1132                        B_average=Paterson(temperature_average);
     1133                        //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2])
     1134                        //                      +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0;
    11311135                        matice->SetB(B_average);
    11321136                }

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

    r1901 r1907  
    335335
    336336%CONTROL
    337 if md.analysis_type==ControlAnalysisEnum,
     337if md.control_analysis,
    338338
    339339        %CONTROL TYPE

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

    r1867 r1907  
    1111%   Examples:
    1212%      md=solve(md,'analysis_type','diagnostic','package','cielo');
    13 %      md=solve(md,'analysis_type','control','package','ice');
    1413%      md=solve(md,'analysis_type','thermal','sub_analysis_type','transient','package','ice');
    1514%      md=solve(md,'analysis_type','thermal','sub_analysis_type','steady','package','cielo');
    … …  
    6766        md=prognostic(md);;
    6867
    69 elseif md.analysis_type==ControlAnalysisEnum,
    70         md=control(md);
    71 
    7268elseif md.analysis_type==ThermalAnalysisEnum,
    7369        md=thermal(md);
    … …  
    9187%Check result is consistent
    9288displaystring(md.debug,'%s\n','checking result consistency');
    93 if ~isresultconsistent(md,options.analysis_type),
    94         disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
    95 end
     89%if ~isresultconsistent(md,options.analysis_type),
     90%       disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
     91%end
    9692
    9793%convert analysis type to string finally

  • issm/trunk/src/m/solutions/cielo/diagnostic.m

    r1647 r1907  
    3333        inputs=inputlist;
    3434        inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes);
     35        if md.control_analysis,
     36                inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes);
     37        end
    3538
    3639        %compute solution
    3740        if ~models.dh.parameters.qmu_analysis,
    38                 %launch core of diagnostic solution.
    39                 results=diagnostic_core(models,inputs);
    40        
    41                 %process results
    42                 if ~isstruct(md.results), md.results=struct(); end
    43                 md.results.diagnostic=processresults(models,results);
     41                if md.control_analysis,
     42                        %launch core of control solution.
     43                        results=control_core(models,inputs);
     44
     45                        %process results
     46                        if ~isstruct(md.results), md.results=struct(); end
     47                        md.results.control=processresults(models,results);
     48                else,
     49                        %launch core of diagnostic solution.
     50                        results=diagnostic_core(models,inputs);
     51
     52                        %process results
     53                        if ~isstruct(md.results), md.results=struct(); end
     54                        md.results.diagnostic=processresults(models,results);
     55                end
    4456        else
    4557                %launch dakota driver for diagnostic core solution

  • issm/trunk/src/m/solutions/cielo/gradjcompute_core.m

    r1850 r1907  
    4444%compute initial velocity from diagnostic_core (horiz+vertical)
    4545displaystring(debug,'%s','          compute initial velocity...');
    46 results_diag=diagnostic_core(models,inputs);
    47 u_g=results_diag.u_g;
     46%results_diag=diagnostic_core(models,inputs);
     47%u_g=results_diag.u_g;
    4848
    4949%Assign output
Note: See TracChangeset for help on using the changeset viewer.