Changeset 1037


Ignore:
Timestamp:
06/22/09 08:04:45 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added core routine in control parallel

Location:
issm/trunk/src/c
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

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

    r962 r1037  
    559559                                        ./parallel/WriteLockFile.cpp\
    560560                                        ./parallel/control.cpp\
     561                                        ./parallel/control_core.cpp\
    561562                                        ./parallel/objectivefunctionC.cpp\
    562563                                        ./parallel/GradJCompute.cpp\

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

    r928 r1037  
    146146                fem_dh=fems+0; //load u_g
    147147                fem_c=fems+0;  //load param_g
     148
     149                /*some flags needed: */
     150                fem_dh->parameters->FindParam((void*)&dim,"dim");
     151                fem_dh->parameters->FindParam((void*)&ishutter,"ishutter");
     152                fem_dh->parameters->FindParam((void*)&isstokes,"isstokes");
     153                fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
    148154        }
    149155

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

    r962 r1037  
    1717int main(int argc,char* *argv){
    1818
    19         int n;
    20         int i;
    21 
    2219        /*I/O: */
    2320        FILE* fid=NULL;
    … …  
    2522        char* outputfilename=NULL;
    2623        char* lockname=NULL;
    27         int   analysis_type;
    28         int   sub_analysis_type;
     24        int   numberofnodes;
     25        int   qmu_analysis=0;
     26
     27        /*Fem models : */
     28        FemModel femmodels[1];
     29
     30        /*Results: */
     31        DataSet* results=NULL;
     32
     33        ParameterInputs* inputs=NULL;
     34        int waitonlock=0;
    2935
    3036        /*Intermediary: */
    31         FemModel femmodel;
    32         Vec u_g=NULL;
    33         ParameterInputs* inputs=NULL;
    34         int waitonlock=0;
    35         double search_scalar;
    36         char* control_type=NULL;
    37         double* fit=NULL;
    38         double* optscal=NULL;
    39         double* u_g_obs=NULL;
    4037        double* u_g_initial=NULL;
    41         int  gsize;
    42         int  numberofnodes;
    43         double*    maxiter=NULL;
    44         double  tolx;
    45         double*   param_g=NULL;
    46         Vec       grad_g=NULL;
    47         Vec       new_grad_g=NULL;
    48         Vec       grad_g_old=NULL;
    49         double*   grad_g_double=NULL;
    50         double    mincontrolconstraint;
    51         double    maxcontrolconstraint;
    52         int       nsteps;
    53         double*   J=NULL;
    54         OptArgs   optargs;
    55         OptPars   optpars;
    56         Param*    param=NULL;
     38        Param*  param=NULL;
    5739
    5840
    … …  
    7658        /*Open handle to data on disk: */
    7759        fid=pfopen(inputfilename,"rb");
    78        
     60
    7961        _printf_("read and create finite element model:\n");
    80         CreateFemModel(&femmodel,fid,"control","");
     62        _printf_("\n   reading diagnostic horiz model data:\n");
     63        CreateFemModel(&femmodels[0],fid,"control","");
    8164
    82         /*Recover parameters used throughout the solution:*/
    83         femmodel.parameters->FindParam((void*)&nsteps,"nsteps");
    84         femmodel.parameters->FindParam((void*)&control_type,"control_type");
    85         femmodel.parameters->FindParam((void*)&fit,"fit");
    86         femmodel.parameters->FindParam((void*)&optscal,"optscal");
    87         femmodel.parameters->FindParam((void*)&maxiter,"maxiter");
    88         femmodel.parameters->FindParam((void*)&tolx,"tolx");
    89         femmodel.parameters->FindParam((void*)&waitonlock,"waitonlock");
    90         femmodel.parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint");
    91         femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint");
    92         femmodel.parameters->FindParam((void*)&param_g,"param_g");
    93         femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs");
    94         femmodel.parameters->FindParam((void*)&analysis_type,"analysis_type");
    95         femmodel.parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
    96         gsize=femmodel.nodes->NumberOfDofs();
    97 
    98         /*Initialize misfit: */
    99         J=(double*)xmalloc(nsteps*sizeof(double));
    100 
    101         /*Initialize inputs:*/
    102         femmodel.parameters->FindParam((void*)&u_g_initial,"u_g");
    103         femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     65        _printf_("initialize inputs:\n");
     66        femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g");
     67        femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    10468
    10569        inputs=new ParameterInputs;
    10670        inputs->Add("velocity",u_g_initial,3,numberofnodes);
    107        
    108         /*erase useless parameters: */
    109         param=(Param*)femmodel.parameters->FindParamObject("param_g");
    110         femmodel.parameters->DeleteObject((Object*)param);
    11171
    112         param=(Param*)femmodel.parameters->FindParamObject("u_g_obs");
    113         femmodel.parameters->DeleteObject((Object*)param);
     72        /*erase velocities: */
     73        param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
     74        femmodels[0].parameters->DeleteObject((Object*)param);
    11475
    115         param=(Param*)femmodel.parameters->FindParamObject("u_g");
    116         femmodel.parameters->DeleteObject((Object*)param);
     76        _printf_("initialize results:\n");
     77        results=new DataSet(ResultsEnum());
    11778
     79        /*are we running the solution sequence, or a qmu wrapper around it? : */
     80        femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
     81        if(!qmu_analysis){
    11882
    119         /*Start looping: */
    120         for(n=0;n<nsteps;n++){
    121                        
    122                 _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
    123                 inputs->Add(control_type,param_g,2,numberofnodes);
    124                 inputs->Add("fit",fit[n]);
     83                /*run control analysis: */
     84                _printf_("call computational core:\n");
     85                control_core(results,femmodels,inputs);
    12586
    126                 /*Update parameters: */
    127                 UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
     87        }
     88        else{
     89                /*run qmu analysis: */
     90                _printf_("calling qmu analysis on control core:\n");
    12891
    129                 _printf_("%s\n","      computing gradJ...");
    130                 grad_g=GradJCompute(inputs,&femmodel,u_g_obs);
    131                 _printf_("%s\n","      done.");
    132                        
    133                 _printf_("%s\n","      normalizing directions...");
    134                 Orthx(&new_grad_g,grad_g,grad_g_old);
    135                 VecFree(&grad_g); grad_g=new_grad_g;
    136                 VecFree(&grad_g_old); grad_g_old=grad_g;
    137                 VecToMPISerial(&grad_g_double,grad_g);
    138                 _printf_("%s\n","      done.");
    139        
    140                 _printf_("%s\n","      optimizing along gradient direction...");
    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;
    142                 optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];
    143                 BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
    144                 _printf_("%s\n","      done.");
    145        
    146                 _printf_("%s\n","      updating parameter using optimized search scalar...");
    147                 for(i=0;i<gsize;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
    148                 _printf_("%s\n","      done.");
    149 
    150                 _printf_("%s\n","      constraining the new distribution...");   
    151                 ControlConstrainx( param_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
    152                 _printf_("%s\n","      done.");
    153 
    154                 //some temporary saving
    155                 /*if (((n+1)%5)==0){
    156                         _printf_("%s\n","      saving temporary results...");
    157                         inputs->Add(control_type,param_g,2,numberofnodes);
    158                         inputs->Add("fit",fit[n]);
    159                         UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
    160                         diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
    161                         OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
    162                         _printf_("%s\n","      done.");
    163                 }*/
    164 
    165                 _printf_("%s%i%s%g\n","      value of misfit J after optimization #",n,": ",J[n]);
    166                
    167                 /*some freeing:*/
    168                 xfree((void**)&grad_g_double);
     92#ifdef _HAVE_DAKOTA_
     93                Qmux(&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
     94#else
     95                throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!");
     96#endif
    16997        }
    17098
    171         /*Write results to disk: */
    172         _printf_("%s\n","      preparing final velocity solution...");
    173         inputs->Add(control_type,param_g,2,numberofnodes);
    174         inputs->Add("fit",fit[n]);
     99        _printf_("process results:\n");
     100        ProcessResults(&results,&femmodels[0],ControlAnalysisEnum());
    175101
    176         /*Update parameters: */
    177         UpdateFromInputsx(femmodel.elements,femmodel.nodes,femmodel.loads, femmodel.materials,inputs);
    178        
    179         diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
    180        
    181         _printf_("%s\n","      saving final results...");
    182         OutputControl(u_g,param_g,J,nsteps,&femmodel,outputfilename);
    183         _printf_("%s\n","      done.");
     102        _printf_("write results to disk:\n");
     103        OutputResults(results,outputfilename);
    184104
    185         /*Write lock file if requested: */
     105        _printf_("write lock file:\n");
     106        femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
    186107        if (waitonlock){
    187108                WriteLockFile(lockname);
    188109        }
    189                        
     110
    190111        _printf_("closing MPI and Petsc\n");
    191112        PetscFinalize();

  • issm/trunk/src/c/parallel/parallel.h

    r962 r1037  
    1313void diagnostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
    1414void prognostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
     15void control_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
    1516
    1617void thermal_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
    … …  
    3738//int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
    3839void OutputResults(DataSet* results,char* filename);
    39 void OutputControl(Vec u_g,double* p_g,double* J,int nsteps,FemModel* fem,char* outputfilename);
    4040void WriteLockFile(char* filename);
    4141
Note: See TracChangeset for help on using the changeset viewer.