Changeset 758


Ignore:
Timestamp:
06/04/09 10:18:31 (16 years ago)
Author:
Mathieu Morlighem
Message:

fixed control method output

Location:
issm/trunk/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r713 r758  
    3535        double* vy_obs=NULL;
    3636        double* control_parameter=NULL;
     37        double* param_g=NULL;
    3738
    3839        double* u_g_obs=NULL;
    … …  
    124125
    125126                u_g_obs=(double*)xcalloc(numberofnodes*2,sizeof(double));
    126                 p_g=(double*)xcalloc(numberofnodes*2,sizeof(double));
     127                param_g=(double*)xcalloc(numberofnodes*2,sizeof(double));
    127128
    128129                for(i=0;i<numberofnodes;i++){
    129130                        u_g_obs[(int)(2*partition[i]+0)]=vx_obs[i]; 
    130131                        u_g_obs[(int)(2*partition[i]+1)]=vy_obs[i]; 
    131                         p_g[(int)(2*partition[i]+0)]=control_parameter[i];
     132                        param_g[(int)(2*partition[i]+0)]=control_parameter[i];
    132133                }
    133134
    … …  
    139140
    140141                count++;
    141                 param= new Param(count,"p_g",DOUBLEVEC);
    142                 param->SetDoubleVec(p_g,2*numberofnodes);
     142                param= new Param(count,"param_g",DOUBLEVEC);
     143                param->SetDoubleVec(param_g,2*numberofnodes);
    143144                parameters->AddObject(param);
    144145
    … …  
    272273        xfree((void**)&vx_obs);
    273274        xfree((void**)&vy_obs);
     275        xfree((void**)&param_g);
    274276        xfree((void**)&pressure);
    275277        xfree((void**)&control_parameter);

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

    r465 r758  
    1515        mxArray* m;
    1616        mxArray* inputs;
    17         mxArray* p_g;
     17        mxArray* param_g;
    1818        mxArray* u_g_obs;
    1919        mxArray* grad_g;
    … …  
    3131struct OptArgs{
    3232        FemModel* femmodel;
    33         double* p_g;
     33        double* param_g;
    3434        double* u_g_obs;
    3535        double* grad_g;

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

    r643 r758  
    2020Result::Result(){
    2121        return;
     22}
     23
     24Result::Result(const Result& result){
     25
     26        id=result.id;
     27        time=result.time;
     28        step=result.step;
     29        size=result.size;
     30       
     31        /*copy constructor: copy dynamically allocated fields: */
     32        if(result.fieldname){
     33                fieldname=(char*)xmalloc((strlen(result.fieldname)+1)*sizeof(char));
     34                strcpy(fieldname,result.fieldname);
     35        }
     36        if(result.field){
     37                VecDuplicatePatch(&field,result.field);
     38                dfield=NULL;
     39        }
     40        if(result.dfield){
     41                dfield=(double*)xmalloc(result.size*sizeof(double));
     42                memcpy(dfield,result.dfield,result.size*sizeof(double));
     43                field=NULL;
     44        }
    2245}
    2346
    … …  
    144167        /*First write field name :*/
    145168        length=(strlen(fieldname)+1)*sizeof(char);
    146 
    147169        fwrite(&length,sizeof(int),1,fid);
    148170        fwrite(fieldname,length,1,fid);
    … …  
    169191        return new Result(*this);
    170192}
     193               
     194

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

    r643 r758  
    2424
    2525                Result();
     26                Result(const Result& result);
    2627                Result(int result_id,double result_time,int result_step,char* result_fieldname,Vec result_field);
    2728                Result(int result_id,double result_time,int result_step,char* result_fieldname,double* result_field,int result_size);

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

    r669 r758  
    1616#define __FUNCT__ "OutputControl"
    1717
    18 void OutputControl(Vec u_g,double* p_g,double* J,int nsteps,FemModel* fem,char* outputfilename){
     18void OutputControl(Vec u_g,double* param_g,double* J,int nsteps,FemModel* fem,char* outputfilename){
    1919
    2020        /*intermediary: */
    … …  
    2626        fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    2727
     28        /*Initialize results*/
    2829        results=new DataSet(ResultsEnum());
    2930
    … …  
    3233        results->AddObject(result);
    3334
    34         result=new Result(results->Size()+1,0,1,"param_g",p_g,numberofnodes);
     35        result=new Result(results->Size()+1,0,1,"param_g",param_g,numberofnodes);
    3536        results->AddObject(result);
    3637
    … …  
    4041        //process results
    4142        ProcessResults(&results,fem,ControlAnalysisEnum());
    42        
     43
    4344        //write results to disk:
    4445        OutputResults(results,outputfilename);
    … …  
    4647        /*Free ressources:*/
    4748        delete results;
    48 
    4949}
    50 

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

    r695 r758  
    4545        /*fem prognostic models: */
    4646        FemModel* fem_p=NULL;
     47
     48        /*fem control models: */
     49        FemModel* fem_c=NULL;
    4750
    4851        /*some parameters*/
    … …  
    6265        double* p_g_serial=NULL;
    6366        double* pressure=NULL;
    64         double* parameter=NULL;
    6567        double* partition=NULL;
    6668        double  yts;
     69
     70        double* param_g=NULL;
     71        double* parameter=NULL;
    6772
    6873        Vec     h_g=NULL;
    … …  
    108113        if(analysis_type==ThermalAnalysisEnum()){
    109114                fem_t=fems+0;
     115        }
     116
     117        if(analysis_type==ControlAnalysisEnum()){
     118                fem_dh=fems+0; //load u_g
     119                fem_c=fems+0;  //load param_g
    110120        }
    111121
    … …  
    209219                        xfree((void**)&partition);
    210220                }
    211 
    212221                else if(strcmp(result->GetFieldName(),"p_g")==0){
    213222                        /*easy, p_g is of size numberofnodes, on 1 dof, just repartition: */
    … …  
    244253                        xfree((void**)&partition);
    245254                }
    246 
    247255                else if(strcmp(result->GetFieldName(),"t_g")==0){
    248256                        /*easy, t_g is of size numberofnodes, on 1 dof, just repartition: */
    … …  
    266274                        xfree((void**)&partition);
    267275                }
    268 
    269276                else if(strcmp(result->GetFieldName(),"m_g")==0){
    270277                        /*easy, m_g is of size numberofnodes, on 1 dof, just repartition: */
    … …  
    289296                        xfree((void**)&partition);
    290297                }
    291 
    292298                else if(strcmp(result->GetFieldName(),"h_g")==0){
    293299                        /*easy, h_g is of size numberofnodes, on 1 dof, just repartition: */
    … …  
    311317                        xfree((void**)&partition);
    312318                }
    313 
    314319                else if(strcmp(result->GetFieldName(),"param_g")==0){
    315320                        /*easy, param_g is of size numberofnodes, on 1 dof, just repartition: */
    316                         result->GetField(&p_g);
    317                         VecToMPISerial(&p_g_serial,p_g);
    318                         fem_p->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    319                         VecToMPISerial(&partition,fem_p->partition);
     321                        result->GetField(&param_g);
     322                        fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     323                        VecToMPISerial(&partition,fem_dh->partition);
    320324
    321325                        parameter=(double*)xmalloc(numberofnodes*sizeof(double));
    322326
    323327                        for(i=0;i<numberofnodes;i++){
    324                                 parameter[i]=p_g_serial[(int)partition[i]];
     328                                parameter[i]=param_g[2*(int)partition[i]];
    325329                        }
    326330                       
    … …  
    330334
    331335                        /*do some cleanup: */
    332                         xfree((void**)&p_g_serial);
    333                         xfree((void**)&partition);
    334                 }
    335 
     336                        xfree((void**)&partition);
     337                }
    336338                else{
    337339                        /*Just copy the result into the new results dataset: */
    338                         newresults->AddObject(result);
     340                        newresults->AddObject(result->copy());
    339341                }
    340342        }

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

    r667 r758  
    4343        double*    maxiter=NULL;
    4444        double  tolx;
    45         double*   p_g=NULL;
     45        double*   param_g=NULL;
    4646        Vec       grad_g=NULL;
    4747        Vec       new_grad_g=NULL;
    … …  
    9090        femmodel.parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint");
    9191        femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint");
    92         femmodel.parameters->FindParam((void*)&p_g,"p_g");
     92        femmodel.parameters->FindParam((void*)&param_g,"param_g");
    9393        femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs");
    9494        femmodel.parameters->FindParam((void*)&analysis_type,"analysis_type");
    … …  
    107107       
    108108        /*erase useless parameters: */
    109         param=(Param*)femmodel.parameters->FindParamObject("p_g");
     109        param=(Param*)femmodel.parameters->FindParamObject("param_g");
    110110        femmodel.parameters->DeleteObject((Object*)param);
    111111
    … …  
    121121                       
    122122                _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
    123                 inputs->Add(control_type,p_g,2,numberofnodes);
     123                inputs->Add(control_type,param_g,2,numberofnodes);
    124124                inputs->Add("fit",fit[n]);
    125125
    … …  
    139139       
    140140                _printf_("%s\n","      optimizing along gradient direction...");
    141                 optargs.femmodel=&femmodel; optargs.p_g=p_g; optargs.u_g_obs=u_g_obs; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n;
     141                optargs.femmodel=&femmodel; optargs.param_g=param_g; optargs.u_g_obs=u_g_obs; optargs.grad_g=grad_g_double; optargs.inputs=inputs;optargs.n=n;
    142142                optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];
    143143                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
    … …  
    145145       
    146146                _printf_("%s\n","      updating parameter using optimized search scalar...");
    147                 for(i=0;i<gsize;i++)p_g[i]=p_g[i]+search_scalar*optscal[n]*grad_g_double[i];
     147                PetscSynchronizedPrintf(MPI_COMM_WORLD,"gsize=%i\n",gsize);
     148                PetscSynchronizedFlush(MPI_COMM_WORLD);
     149
     150                for(i=0;i<gsize;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
    148151                _printf_("%s\n","      done.");
    149152
    150153                _printf_("%s\n","      constraining the new distribution...");   
    151                 ControlConstrainx( p_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
     154                ControlConstrainx( param_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
    152155                _printf_("%s\n","      done.");
    153156       
    … …  
    155158                if (((n+1)%5)==0){
    156159                        _printf_("%s\n","      saving temporary results...");
    157                         inputs->Add(control_type,p_g,2,numberofnodes);
     160                        inputs->Add(control_type,param_g,2,numberofnodes);
    158161                        inputs->Add("fit",fit[n]);
    159162                        diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
    160                         OutputControl(u_g,p_g,J,nsteps,&femmodel,outputfilename);
     163                        OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
    161164                        _printf_("%s\n","      done.");
    162165                }
    … …  
    170173        /*Write results to disk: */
    171174        _printf_("%s\n","      preparing final velocity solution...");
    172         inputs->Add(control_type,p_g,2,numberofnodes);
     175        inputs->Add(control_type,param_g,2,numberofnodes);
    173176        inputs->Add("fit",fit[n]);
    174177
    … …  
    179182       
    180183        _printf_("%s\n","      saving final results...");
    181         OutputControl(u_g,p_g,J,nsteps,&femmodel,outputfilename);
     184        OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
    182185        _printf_("%s\n","      done.");
    183186

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

    r465 r758  
    2424        /*parameters: */
    2525        FemModel* femmodel=NULL;
    26         double* p_g=NULL;
     26        double* param_g=NULL;
    2727        double* u_g_obs=NULL;
    2828        double* grad_g=NULL;
    … …  
    3737        double  maxcontrolconstraint;
    3838        char*   control_type=NULL;
    39         double* p_g_copy=NULL;
     39        double* param_g_copy=NULL;
    4040        int     analysis_type;
    4141        int     sub_analysis_type;
    … …  
    4747        /*Recover parameters: */
    4848        femmodel=optargs->femmodel;
    49         p_g=optargs->p_g;
     49        param_g=optargs->param_g;
    5050        u_g_obs=optargs->u_g_obs;
    5151        grad_g=optargs->grad_g;
    … …  
    6363        femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    6464
    65         /*First copy p_g so we don't modify it: */
    66         p_g_copy=(double*)xmalloc(gsize*sizeof(double));
    67         memcpy(p_g_copy,p_g,gsize*sizeof(double));
     65        /*First copy param_g so we don't modify it: */
     66        param_g_copy=(double*)xmalloc(gsize*sizeof(double));
     67        memcpy(param_g_copy,param_g,gsize*sizeof(double));
    6868
    69         /*First, update p_g using search_scalar: */
    70         for(i=0;i<gsize;i++)p_g_copy[i]=p_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
     69        /*First, update param_g using search_scalar: */
     70        for(i=0;i<gsize;i++)param_g_copy[i]=param_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
    7171
    7272        /*Constrain:*/
    73         ControlConstrainx( p_g_copy,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
     73        ControlConstrainx( param_g_copy,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
    7474
    7575        /*Add new parameter to inputs: */
    76         inputs->Add(control_type,p_g_copy,2,numberofnodes);
     76        inputs->Add(control_type,param_g_copy,2,numberofnodes);
    7777
    7878        //Run diagnostic with updated parameters.
    … …  
    8989        xfree((void**)&optscal);
    9090        xfree((void**)&control_type);
    91         xfree((void**)&p_g_copy);
     91        xfree((void**)&param_g_copy);
    9292        xfree((void**)&u_g_double);
    9393

  • issm/trunk/src/c/shared/Numerics/OptFunc.cpp

    r465 r758  
    3232        inputs[1]=optargs->m;
    3333        inputs[2]=optargs->inputs;
    34         inputs[3]=optargs->p_g;
     34        inputs[3]=optargs->param_g;
    3535        inputs[4]=optargs->u_g_obs;
    3636        inputs[5]=optargs->grad_g;
Note: See TracChangeset for help on using the changeset viewer.