Changeset 643


Ignore:
Timestamp:
05/29/09 15:32:27 (15 years ago)
Author:
Eric.Larour
Message:

New framework for results output for parallel solutions . Watch out for bugs...

Location:
issm/trunk/src
Files:
5 added
5 deleted
19 edited

Legend:

Unmodified
Added
Removed

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

    r493 r643  
    1111int ParamEnum(void){                    return           5; }
    1212int RgbEnum(void){                      return           6; }
     13int ResultEnum(void){                   return           7; }
    1314
    1415/*analysis types: */
    … …  
    3839int MaterialsEnum(void){                return          44; }
    3940int ParametersEnum(void){               return          45; }
     41int ResultsEnum(void){                  return          46; }
    4042
    4143/*Elements: */

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

    r472 r643  
    2525int RgbEnum(void);
    2626int InputEnum(void);
     27int ResultEnum(void);
    2728
    2829/*formulations: */
    … …  
    6061int MaterialsEnum(void);
    6162int ParametersEnum(void);
     63int ResultsEnum(void);
    6264
    6365/*Functions on enums: */

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

    r586 r643  
    297297                                        ./objects/Node.h\
    298298                                        ./objects/Node.cpp\
     299                                        ./objects/Result.h\
     300                                        ./objects/Result.cpp\
    299301                                        ./objects/DakotaPlugin.h\
    300302                                        ./objects/DakotaPlugin.cpp\
    … …  
    534536                                        ./parallel/diagnostic_core_nonlinear.cpp\
    535537                                        ./parallel/thermal_core.cpp\
     538                                        ./parallel/thermal_core_nonlinear.cpp\
    536539                                        ./parallel/CreateFemModel.cpp\
    537                                         ./parallel/OutputDiagnostic.cpp\
    538540                                        ./parallel/WriteLockFile.cpp\
    539541                                        ./parallel/control.cpp\
    540                                         ./parallel/OutputControl.cpp\
    541                                         ./parallel/OutputThermal.cpp\
    542                                         ./parallel/OutputPrognostic.cpp\
    543542                                        ./parallel/objectivefunctionC.cpp\
    544543                                        ./parallel/GradJCompute.cpp\
    545544                                        ./parallel/qmu.cpp\
    546545                                        ./parallel/SpawnCore.cpp\
    547                                         ./parallel/DakotaResponses.cpp
     546                                        ./parallel/DakotaResponses.cpp\
     547                                        ./parallel/ProcessResults.cpp\
     548                                        ./parallel/OutputResults.cpp
    548549
    549550libpISSM_a_CXXFLAGS = -fPIC -D_PARALLEL_   -D_C_
    … …  
    552553bin_PROGRAMS =
    553554else
    554 bin_PROGRAMS = diagnostic.exe
     555bin_PROGRAMS = diagnostic.exe  thermal.exe
    555556#control.exe thermal.exe prognostic.exe
    556557endif

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

    r616 r643  
    2525        int      numberofdofspernode;
    2626        int      dim;
    27         int*     epart=NULL;
    28         int*     part=NULL;
     27        int*  epart=NULL;
     28        int*  part=NULL;
     29        double*  dpart=NULL;
     30        int      elements_width;
    2931
    3032
    … …  
    7678        param->SetInteger(model->ismacayealpattyn);
    7779        parameters->AddObject(param);
     80
    7881
    7982        count++;
    … …  
    220223                char*  descriptor=NULL;
    221224                char*  tag=NULL;
    222                 int    elements_width; //size of elements
    223225
    224226                descriptors=(char**)xmalloc(model->numberofresponses*sizeof(char*));
    … …  
    246248                xfree((void**)&descriptors);
    247249
    248                 /*Width of elements: */
    249                 if(strcmp(model->meshtype,"2d")==0){
    250                         elements_width=3; //tria elements
    251                 }
    252                 else{
    253                         elements_width=6; //penta elements
    254                 }
    255 
    256250                #ifdef _PARALLEL_
    257251                /*partition grids in model->qmu_npart parts: */
    … …  
    259253                if(strcmp(model->meshtype,"2d")==0){
    260254                        ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     255                        elements_width=3; //tria elements
    261256                }
    262257                else{
    263258                        ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
     259                        elements_width=6; //penta elements
    264260                }
    265261
    266262                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];
    267266
    268267                count++;
    269268                param= new Param(count,"qmu_part",DOUBLEVEC);
    270                 param->SetDoubleVec((double*)part,model->numberofnodes);
     269                param->SetDoubleVec(dpart,model->numberofnodes);
    271270                parameters->AddObject(param);
    272271
    … …  
    276275                xfree((void**)&epart);
    277276                xfree((void**)&part);
     277                xfree((void**)&dpart);
    278278
    279279                #endif

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

    r586 r643  
    1919#include "./Beam.h"
    2020#include "./Spc.h"
     21#include "./Result.h"
    2122#include "./Rgb.h"
    2223#include "./Icefront.h"

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

    r587 r643  
    99#endif
    1010
     11#include "../DataSet/DataSet.h"   
     12#include "../shared/shared.h"
     13
    1114#undef __FUNCT__
    1215#define __FUNCT__ "DakotaResponses"
    13    
    14 void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Vec u_g,Vec p_g,int analysis_type,int sub_analysis_type){
     16void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,DataSet* results,int analysis_type,int sub_analysis_type){
    1517
    1618        int i;

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

    r587 r643  
    3434        double* qmu_part=NULL;
    3535        int     qmu_npart;
     36        DataSet* results=NULL;
    3637
    3738        extern int my_rank;
    … …  
    9192
    9293        /*Modify core inputs to reflect the dakota variables inputs: */
    93         inputs->UpdateFromDakota(variables,variables_descriptors,numvariables,qmu_part,qmu_npart);
     94        //inputs->UpdateFromDakota(variables,variables_descriptors,numvariables,qmu_part,qmu_npart);
    9495
    9596        /*Run the analysis core solution sequence, with the updated inputs: */
    … …  
    9798
    9899                _printf_("Starting diagnostic core\n");
    99                 diagnostic_core(&u_g,&p_g,femmodels,inputs);
     100                diagnostic_core(results,femmodels,inputs);
    100101
    101102        }
    … …  
    105106
    106107        /*compute responses on cpu 0: dummy for now! */
    107         DakotaResponses(responses,responses_descriptors,numresponses,u_g,p_g,analysis_type,sub_analysis_type);
     108        DakotaResponses(responses,responses_descriptors,numresponses,results,analysis_type,sub_analysis_type);
    108109
    109110        /*Free ressources:*/

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

    r472 r643  
    158158                        inputs->Add("fit",fit[n]);
    159159                        diagnostic_core_nonlinear(&u_g,NULL,NULL,&femmodel,inputs,analysis_type,sub_analysis_type);
    160                         OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
     160                        //OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
    161161                        _printf_("%s\n","      done.");
    162162                }
    … …  
    179179       
    180180        _printf_("%s\n","      saving final results...");
    181         OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
     181        //OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
    182182        _printf_("%s\n","      done.");
    183183

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

    r615 r643  
    1515#endif
    1616
     17
    1718int main(int argc,char* *argv){
    1819       
    … …  
    2829        /*Fem models : */
    2930        FemModel femmodels[5];
     31
     32        /*Results: */
     33        DataSet* results=NULL;
     34        DataSet* newresults=NULL;
    3035       
    31         Vec u_g=NULL;
    32         Vec p_g=NULL;
    3336        ParameterInputs* inputs=NULL;
    3437        int waitonlock=0;
    … …  
    7073        CreateFemModel(&femmodels[4],fid,"slope_compute","");
    7174
     75
    7276        _printf_("initialize inputs:\n");
    7377        femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g");
    … …  
    7680        inputs=new ParameterInputs;
    7781        inputs->Add("velocity",u_g_initial,3,numberofnodes);
    78 
     82       
    7983        /*erase velocities: */
    8084        param=(Param*)femmodels[0].parameters->FindParamObject("u_g");
    8185        femmodels[0].parameters->DeleteObject((Object*)param);
    8286
     87        _printf_("initialize results:\n");
     88        results=new DataSet(ResultsEnum());
     89
    8390        /*are we running the solutoin sequence, or a qmu wrapper around it? : */
    8491        femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
    8592        if(!qmu_analysis){
     93
    8694                /*run diagnostic analysis: */
    8795                _printf_("call computational core:\n");
    88                 diagnostic_core(&u_g,&p_g,femmodels,inputs);
     96                diagnostic_core(results,femmodels,inputs);
    8997
    90                 _printf_("write results to disk:\n");
    91                 OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename);
    9298        }
    9399        else{
    94100                /*run qmu analysis: */
    95101                _printf_("calling qmu analysis on diagnostic core:\n");
     102               
    96103                qmu(qmuname,&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
    97104        }
     105
     106        _printf_("process results:\n");
     107        ProcessResults(&newresults,results,&femmodels[0],DiagnosticAnalysisEnum()); delete results;
     108       
     109        _printf_("write results to disk:\n");
     110        OutputResults(newresults,outputfilename);
    98111
    99112        _printf_("write lock file:\n");

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

    r586 r643  
    1313#include "../issm.h"
    1414
    15 void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs){
     15void diagnostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs){
    1616
    1717        extern int my_rank;
    … …  
    2323        FemModel* fem_ds=NULL;
    2424        FemModel* fem_sl=NULL;
     25
     26        /*output: */
     27        Result* result=NULL;
    2528
    2629        /*solutions: */
    … …  
    171174        }
    172175       
    173         /*Assign output pointers: */
    174         *pug=ug;
    175         *ppg=pg;
     176        /*Plug results into output dataset: */
     177        result=new Result(results->Size()+1,0,1,"u_g",ug);
     178        results->AddObject(result);
     179        result=new Result(results->Size()+1,0,1,"p_g",pg);
     180        results->AddObject(result);
    176181}

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

    r586 r643  
    1111Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel,double* u_g_obs);
    1212
    13 void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs);
     13void diagnostic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
     14
     15void thermal_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
     16void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     17
    1418void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1519void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int  analysis_type,int sub_analysis_type);
    16 void thermal_core(Vec* pt_g,double* pmelting_offset,FemModel* femmodel,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1720
    1821//int GradJOrth(WorkspaceParams* workspaceparams);
    … …  
    2831
    2932//int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
    30 void OutputDiagnostic(Vec u_g,Vec p_g, FemModel* femmodels,char* filename);
    31 void OutputThermal(Vec* t_g,Vec* m_g, double* time,FemModel* femmodels,char* filename);
    32 void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets);
    33 void OutputPrognostic(Vec h_g,FemModel* femmodel,char* filename);
     33void OutputResults(DataSet* results,char* filename);
    3434void WriteLockFile(char* filename);
    3535
    … …  
    3838void qmu(const char* dakota_input_file,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    3939void SpawnCore(double* responses,double* variables,char** variable_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    40 void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Vec u_g,Vec p_g,int analysis_type,int sub_analysis_type);
     40void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,DataSet* results,int analysis_type,int sub_analysis_type);
     41void ProcessResults(DataSet** pnewresults,DataSet* results,FemModel* fems,int analysis_type);
    4142
    4243#endif

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

    r619 r643  
    9191
    9292        _printf_("write results to disk:\n");
    93         OutputPrognostic(h_g,&fem,outputfilename);
     93        //OutputPrognostic(h_g,&fem,outputfilename);
    9494
    9595        _printf_("write lock file:\n");

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

    r614 r643  
    3232        double* u_g=NULL;
    3333        double* p_g=NULL;
    34         double* time=NULL;
    35         double  dt;
    36         double  ndt;
    37         int     nsteps;
    38         int     debug=0;
    3934
    40         /*solution vectors: */
    41         Vec* t_g=NULL;
    42         Vec* m_g=NULL;
    43 
     35        /*Results: */
     36        DataSet* results=NULL;
     37        DataSet* newresults=NULL;
     38       
    4439        ParameterInputs* inputs=NULL;
    4540        Param*           param=NULL;
     41        double  dt;
    4642
    4743        int    waitonlock=0;
    48         int    sub_analysis_type;
    49         double melting_offset;
    50        
     44               
    5145        MODULEBOOT();
    5246
    … …  
    7771        femmodels[0].parameters->FindParam((void*)&p_g,"p_g");
    7872        femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     73        femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
    7974        femmodels[0].parameters->FindParam((void*)&dt,"dt");
    80         femmodels[0].parameters->FindParam((void*)&ndt,"ndt");
    81         femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock");
    82         femmodels[0].parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
    83         femmodels[0].parameters->FindParam((void*)&debug,"debug");
    8475
    8576        inputs=new ParameterInputs;
    … …  
    9990        femmodels[1].parameters->DeleteObject((Object*)param);
    10091
    101         if(sub_analysis_type==SteadyAnalysisEnum()){
     92       
     93        _printf_("call computational core:\n");
     94        thermal_core(results,femmodels,inputs);
    10295
    103                 time=(double*)xmalloc(sizeof(double));
    104                 time[0]=0;
    105 
    106                 /*allocate t_g and m_g arrays: */
    107                 t_g=(Vec*)xmalloc(sizeof(Vec));
    108                 m_g=(Vec*)xmalloc(sizeof(Vec));
    109 
    110                 if(debug)_printf_("computing temperatures:\n");
    111                 thermal_core(&t_g[0],&melting_offset,&femmodels[0],inputs,ThermalAnalysisEnum(),SteadyAnalysisEnum());
    112                 inputs->Add("temperature",t_g[0],1,numberofnodes);
    113                 inputs->Add("melting_offset",melting_offset);
    114                
    115                 if(debug)_printf_("computing melting:\n");
    116                 diagnostic_core_linear(&m_g[0],&femmodels[1],inputs,MeltingAnalysisEnum(),SteadyAnalysisEnum());
    117         }
    118         else{
    119                
    120                 nsteps=(int)(ndt/dt);
    121                 time=(double*)xmalloc((nsteps+1)*sizeof(double));
    122                 time[0]=0;
    123 
    124                 /*allocate t_g and m_g arrays: */
    125                 t_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
    126                 m_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
    127 
    128                 //initialize temperature and melting
    129                 femmodels[0].parameters->FindParam((void*)&t_g[0],"t_g");
    130                 femmodels[1].parameters->FindParam((void*)&m_g[0],"m_g");
    131 
    132                 //erase temperature and melting embedded in parameters
    133                 param=(Param*)femmodels[0].parameters->FindParamObject("t_g");
    134                 femmodels[0].parameters->DeleteObject((Object*)param);
    135                 param=(Param*)femmodels[1].parameters->FindParamObject("m_g");
    136                 femmodels[1].parameters->DeleteObject((Object*)param);
    137 
    138                 for(i=0;i<nsteps;i++){
    139                         if(debug)_printf_("time step: %i/%i\n",i+1,nsteps);
    140                         time[i+1]=(i+1)*dt;
    141                        
    142                         if(debug)_printf_("computing temperatures:\n");
    143                         inputs->Add("temperature",t_g[i],1,numberofnodes);
    144                         thermal_core(&t_g[i+1],&melting_offset,&femmodels[0],inputs,ThermalAnalysisEnum(),TransientAnalysisEnum());
    145                        
    146                         if(debug)_printf_("computing melting:\n");
    147                         inputs->Add("temperature",t_g[i+1],1,numberofnodes);
    148                         inputs->Add("melting_offset",melting_offset);
    149                         diagnostic_core_linear(&m_g[i+1],&femmodels[1],inputs,MeltingAnalysisEnum(),TransientAnalysisEnum());
    150                 }
    151         }
     96        _printf_("process results:\n");
     97        ProcessResults(&newresults,results,&femmodels[0],ThermalAnalysisEnum()); delete results;
    15298
    15399        _printf_("write results to disk:\n");
    154         OutputThermal(t_g,m_g,time,&femmodels[0],outputfilename);
     100        OutputResults(newresults,outputfilename);
    155101
    156102        _printf_("write lock file:\n");

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

    r618 r643  
    88#include "../toolkits/toolkits.h"
    99#include "../objects/objects.h"
     10#include "../shared/shared.h"
    1011#include "../EnumDefinitions/EnumDefinitions.h"
     12#include "./parallel.h"
    1113#include "../issm.h"
    1214
    13 void thermal_core(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
     15void thermal_core(DataSet* results,FemModel* fems, ParameterInputs* inputs){
    1416
    15         /*solution : */
    16         Vec tg=NULL;
    17         Vec tf=NULL;
     17        extern int my_rank;
     18        int i;
     19
     20        /*fem models: */
     21        FemModel* fem_t=NULL;
     22        FemModel* fem_m=NULL;
     23
     24        /*output: */
     25        Result* result=NULL;
     26
     27        /*solutions vectors: */
     28        Vec* t_g=NULL;
     29        Vec* m_g=NULL;
     30        double* time=NULL;
     31
     32        /*flags: */
     33        int debug=0;
     34        int numberofdofspernode;
     35        int numberofnodes;
     36        int nsteps;
     37        double  dt;
     38        double  ndt;
     39
     40        int    sub_analysis_type;
    1841        double melting_offset;
     42       
     43        Param*           param=NULL;
    1944
    20         /*intermediary: */
    21         Mat Kgg_nopenalty=NULL;
    22         Mat Kgg=NULL;
    23         Mat Kff=NULL;
    24         Mat Kfs=NULL;
    25         Vec pg_nopenalty=NULL;
    26         Vec pg=NULL;
    27         Vec pf=NULL;
     45        /*recover fem models: */
     46        fem_t=fems+0;
     47        fem_m=fems+1;
    2848
    29         int converged;
    30         int constraints_converged;
    31         int num_unstable_constraints;
    32         int count;
    33         int numberofnodes;
    34         int min_thermal_constraints;
     49        //first recover parameters common to all solutions
     50        fem_t->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
     51        fem_t->parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
     52        fem_t->parameters->FindParam((void*)&debug,"debug");
     53        fem_t->parameters->FindParam((void*)&dt,"dt");
     54        fem_t->parameters->FindParam((void*)&ndt,"ndt");
    3555
    36         /*parameters:*/
    37         int kflag,pflag,connectivity,numberofdofspernode;
    38         char* solver_string=NULL;
    39         int debug=0;
    40         int lowmem=0;
     56        if(sub_analysis_type==SteadyAnalysisEnum()){
    4157
    42         /*Recover parameters: */
    43         kflag=1; pflag=1;
     58                time=(double*)xmalloc(sizeof(double));
     59                time[0]=0;
    4460
    45         fem->parameters->FindParam((void*)&connectivity,"connectivity");
    46         fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
    47         fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
    48         fem->parameters->FindParam((void*)&solver_string,"solverstring");
    49         fem->parameters->FindParam((void*)&debug,"debug");
    50         fem->parameters->FindParam((void*)&lowmem,"lowmem");
    51         fem->parameters->FindParam((void*)&min_thermal_constraints,"min_thermal_constraints");
     61                /*allocate t_g and m_g arrays: */
     62                t_g=(Vec*)xmalloc(sizeof(Vec));
     63                m_g=(Vec*)xmalloc(sizeof(Vec));
    5264
    53         count=1;
    54         converged=0;
    55         for(;;){
     65                if(debug)_printf_("computing temperatures:\n");
     66                thermal_core_nonlinear(&t_g[0],&melting_offset,fem_t,inputs,ThermalAnalysisEnum(),SteadyAnalysisEnum());
     67                inputs->Add("temperature",t_g[0],1,numberofnodes);
     68                inputs->Add("melting_offset",melting_offset);
     69               
     70                if(debug)_printf_("computing melting:\n");
     71                diagnostic_core_linear(&m_g[0],fem_m,inputs,MeltingAnalysisEnum(),SteadyAnalysisEnum());
     72        }
     73        else{
     74               
     75                nsteps=(int)(ndt/dt);
     76                time=(double*)xmalloc((nsteps+1)*sizeof(double));
     77                time[0]=0;
    5678
    57                 if(debug)_printf_("%s\n","starting direct shooting method");
     79                /*allocate t_g and m_g arrays: */
     80                t_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
     81                m_g=(Vec*)xmalloc((nsteps+1)*sizeof(Vec));
    5882
    59                 /*Update parameters: */
    60                 UpdateFromInputsx(fem->elements,fem->nodes,fem->loads, fem->materials,inputs);
     83                //initialize temperature and melting
     84                fem_t->parameters->FindParam((void*)&t_g[0],"t_g");
     85                fem_m->parameters->FindParam((void*)&m_g[0],"m_g");
    6186
    62                 //*Generate system matrices
    63                 if (!lowmem){
     87                //erase temperature and melting embedded in parameters
     88                param=(Param*)fem_t->parameters->FindParamObject("t_g");
     89                fem_t->parameters->DeleteObject((Object*)param);
     90                param=(Param*)fem_m->parameters->FindParamObject("m_g");
     91                fem_m->parameters->DeleteObject((Object*)param);
    6492
    65                         /*Compute Kgg_nopenalty and pg_nopenalty once for all: */
    66                         if (count==1){
    67                                 SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type);
    68                         }
     93                for(i=0;i<nsteps;i++){
     94                        if(debug)_printf_("time step: %i/%i\n",i+1,nsteps);
     95                        time[i+1]=(i+1)*dt;
     96                       
     97                        if(debug)_printf_("computing temperatures:\n");
     98                        inputs->Add("temperature",t_g[i],1,numberofnodes);
     99                        thermal_core_nonlinear(&t_g[i+1],&melting_offset,fem_t,inputs,ThermalAnalysisEnum(),TransientAnalysisEnum());
     100                       
     101                        if(debug)_printf_("computing melting:\n");
     102                        inputs->Add("temperature",t_g[i+1],1,numberofnodes);
     103                        inputs->Add("melting_offset",melting_offset);
     104                        diagnostic_core_linear(&m_g[i+1],fem_m,inputs,MeltingAnalysisEnum(),TransientAnalysisEnum());
     105                }
     106        }
     107       
     108        /*Plug results into output dataset: */
     109        if(sub_analysis_type==SteadyAnalysisEnum()){
     110                result=new Result(results->Size()+1,0,1,"t_g",t_g[0]);
     111                results->AddObject(result);
     112               
     113                result=new Result(results->Size()+1,0,1,"m_g",m_g[0]);
     114                results->AddObject(result);
     115        }
     116        else{
     117                for(i=0;i<nsteps+1;i++){
     118                        result=new Result(results->Size()+1,time[i],i+1,"t_g",t_g[i]);
     119                        results->AddObject(result);
    69120
    70                         /*Copy K_gg_nopenalty into Kgg, same for pg: */
    71                         Kgg=(Mat)xmalloc(sizeof(Mat));
    72                         MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg);
    73                         pg=(Vec)xmalloc(sizeof(Vec));
    74                         VecDuplicate(pg_nopenalty,&pg);VecCopy(pg_nopenalty,pg);
    75 
    76                         //apply penalties each time
    77                         PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type);
     121                        result=new Result(results->Size()+1,time[i],i+1,"m_g",m_g[i]);
     122                        results->AddObject(result);
    78123                }
    79                 else{
    80                         SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type,sub_analysis_type);
    81                         //apply penalties
    82                         PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,inputs,analysis_type,sub_analysis_type);
    83                 }
    84 
    85                 /*!Reduce matrix from g to f size:*/
    86                 Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
    87 
    88                 /*Free ressources: */
    89                 MatFree(&Kgg);
    90        
    91                 if (debug) _printf_("   reducing load from g to f set\n");
    92                 /*!Reduce load from g to f size: */
    93                 Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);
    94 
    95                 //no need for pg and Kfs anymore
    96                 VecFree(&pg);
    97                 MatFree(&Kfs);
    98 
    99                 /*Solve: */
    100                 if(debug)_printf_("%s\n","solving");
    101                 Solverx(&tf, Kff, pf, tf, solver_string);
    102        
    103                 //no need for Kff and pf anymore
    104                 MatFree(&Kff);VecFree(&pf);
    105 
    106                 if (debug) _printf_("   merging solution from f to g set\n");
    107                 //Merge back to g set
    108                 Mergesolutionfromftogx(&tg, tf,fem->Gmn,fem->ys,fem->nodesets);
    109 
    110                 //Deal with penalty loads
    111                 if (debug) _printf_("   penalty constraints\n");
    112                 inputs->Add("temperature",tg,numberofdofspernode,numberofnodes);
    113                
    114                 PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->loads,fem->materials,inputs,analysis_type,sub_analysis_type);
    115 
    116                 if (!converged){
    117                         if(debug)_printf_("%s%i\n","   #unstable constraints = ",num_unstable_constraints);
    118                         if (num_unstable_constraints <= min_thermal_constraints)converged=1;
    119                 }
    120                 count++;
    121                
    122                 if(converged==1)break;
    123124        }
    124125
    125         /*Assign output pointers: */
    126         *ptg=tg;
    127         *pmelting_offset=melting_offset;
    128126}

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

    r586 r643  
    276276        md.qmu_analysis=0;
    277277        md.iresp=1;
    278 
    279278        md.npart=0;
     279
     280        %results
     281        md.results=NaN;
    280282
    281283        %Ice solver string

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

    r1 r643  
    1 function field=ReadData(fid)
     1function result=ReadData(fid)
    22%READDATA - ...
    33%
    … …  
    77
    88%read field
    9 [type,count]=fread(fid,1,'double');
     9[length,count]=fread(fid,1,'int');
    1010
    1111if count==0,
    12         field=NaN;
     12        result=struct([]);
    1313else
    14         if type==0,
    15                 %string
    16                 stringsize=fread(fid,1,'double');
    17                 field=char(fread(fid,stringsize,'char'));
    18                 field=field(1:end-1)';
    19         elseif type==1,
    20                 %matrix
    21                 M=fread(fid,1,'double');
    22                 N=fread(fid,1,'double');
    23                 field=fread(fid,[M,N],'double');
    24         elseif ((type==2) || (type==3)),
    25                 field=fread(fid,1,'double');
    26         else
    27                 error('ReadData error message: data type not supported yet!');
    28         end
     14        fieldname=fread(fid,length,'char');
     15        fieldname=fieldname(1:end-1)';
     16        fieldname=char(fieldname);
     17        time=fread(fid,1,'double');
     18        step=fread(fid,1,'int');
     19        ssize=fread(fid,1,'int');
     20        field=fread(fid,ssize,'double');
     21
     22        result.fieldname=fieldname;
     23        result.time=time;
     24        result.step=step;
     25        result.field=field;
    2926end

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

    r640 r643  
    1313
    1414
    15 results=parseresultsfromdisk(filename);
    16 
    17 %First get solution type
    18 analysis_type=results{1};
    19 
    20 if ~strcmpi(analysis_type,md.analysis_type),
    21         error(['loadresultsfromdisk error message: trying  to load results from disk for analysis_type: ',md.analysis_type,' when results are from analysis_type: ',analysis_type]);
    22 end
    23 
    24 %Get part
    25 part=results{2};
    26 
    27 %now to specialized reading
    28 if strcmpi(analysis_type,'diagnostic'),
    29 
    30         %Get u_g
    31         u_g=results{3};
    32         p_g=results{4};
    33 
    34         if strcmpi(md.type,'2d'),
    35                 gsize=md.numberofgrids*2;
    36                 %Used to recover velocities
    37                 indx=1:2:gsize;
    38                 indy=2:2:gsize;
    39                 indx=indx(part);
    40                 indy=indy(part);
    41 
    42                 %Recover velocity
    43                 md.vx=u_g(indx)*md.yts;
    44                 md.vy=u_g(indy)*md.yts;
    45                 md.vel=sqrt(md.vx.^2+md.vy.^2);
    46                 md.pressure=p_g(part);
    47         else
    48                 %Used to recover velocities
    49                 gsize=length(u_g);
    50                 offset=gsize/md.numberofgrids;
    51                 indx=1:offset:gsize;
    52                 indy=2:offset:gsize;
    53                 indz=3:offset:gsize;
    54                 indx=indx(part);
    55                 indy=indy(part);
    56                 indz=indz(part);
    57 
    58                 %Recover velocity
    59                 md.vx=u_g(indx)*md.yts;
    60                 md.vy=u_g(indy)*md.yts;
    61                 md.vz=u_g(indz)*md.yts;
    62                 md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
    63                 md.pressure=p_g(part);
    64         end
    65 
    66 elseif strcmpi(analysis_type,'control'),
    67 
    68         %Get u_g
    69         u_g=results{3};
    70         gsize=length(u_g);
    71 
    72         %Used to recover velocities
    73         indx=1:2:gsize;
    74         indy=2:2:gsize;
    75         indx=indx(part);
    76         indy=indy(part);
    77 
    78         %Recover velocity
    79         md.cont_vx=u_g(indx)*md.yts;
    80         md.cont_vy=u_g(indy)*md.yts;
    81         md.cont_vel=sqrt(md.cont_vx.^2+md.cont_vy.^2);
    82        
    83         %recover parameter
    84         cont_parameter=results{4};
    85         cont_parameter=cont_parameter(indx);
    86         md.cont_parameter=cont_parameter;
    87        
    88         %read J
    89         if(md.nsteps~=results{5}),
    90                 error('output from control method incompatible with model');
    91         end
    92         md.cont_J=results{6};
    93 
    94 elseif strcmpi(analysis_type,'thermal'),
    95 
    96         %read results
    97         time=results{4};
    98         t_g=results{5};
    99         m_g=results{6};
    100 
    101         if strcmpi(md.sub_analysis_type,'steady');
    102 
    103                 %Recover temperature and melting
    104                 md.temperature=t_g(part);
    105                 md.melting=m_g(part)*md.yts;
    106 
    107         else
    108 
    109                 %Recover temperature and melting
    110                 error('not implemented yet')
    111 
    112         end
    113 
    114 elseif strcmpi(analysis_type,'prognostic'),
    115 
    116         %read results
    117         h_g=results{3};
    118         md.new_thickness=h_g(part);
    119 
    120 else
    121         error(['loadresultsfromdisk error message: unknow solution type ',analysis_type]);
    122 end
    123 
     15md.results=parseresultsfromdisk(filename);
    12416
    12517%Check result is consistent

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

    r1 r643  
    1111end
    1212
    13 results={};
     13results=struct();
    1414
    1515%Read fields until the end of the file.
    … …  
    1717
    1818        result=ReadData(fid);
    19         if isnan(result),
     19        if isempty(result),
    2020                break;
    2121        else
    22                 results{end+1}=result;
     22                eval(['results(' num2str(result.step) ').' result.fieldname '=result.field;']);
     23                eval(['results(' num2str(result.step) ').time=result.time;']);
     24                eval(['results(' num2str(result.step) ').step=result.step;']);
    2325        end
    2426end

  • issm/trunk/src/m/utils/Nightly/testsgetanalysis.m

    r517 r643  
    1515if strcmpi(string,'diagnostic'),
    1616        analysis_type='diagnostic';
    17         sub_analysis_type='';
     17        sub_analysis_type='none';
    1818
    1919elseif strcmpi(string,'prognostic'),
    2020        analysis_type='prognostic';
    21         sub_analysis_type='';
     21        sub_analysis_type='none';
    2222
    2323elseif strcmpi(string,'thermalsteady'),
    … …  
    3131elseif strcmpi(string,'transient'),
    3232        analysis_type='transient';
    33         sub_analysis_type='';
     33        sub_analysis_type='none';
    3434
    3535else
Note: See TracChangeset for help on using the changeset viewer.