Changeset 46


Ignore:
Timestamp:
04/24/09 15:31:51 (16 years ago)
Author:
Eric.Larour
Message:

Parallel control method

Location:
issm/trunk/src/c
Files:
3 added
3 deleted
14 edited

Legend:

Unmodified
Added
Removed

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

    r9 r46  
    7171                                        ./shared/Numerics/numerics.h\
    7272                                        ./shared/Numerics/GaussPoints.h\
     73                                        ./shared/Numerics/GaussPoints.cpp\
    7374                                        ./shared/Numerics/BrentSearch.cpp\
    7475                                        ./shared/Numerics/OptFunc.cpp\
    7576                                        ./shared/Numerics/extrema.cpp\
    76                                         ./shared/Numerics/GaussPoints.cpp\
    7777                                        ./shared/Exceptions/exceptions.h\
    7878                                        ./shared/Exceptions/Exceptions.cpp\
    … …  
    419419                                        ./parallel/CreateFemModel.cpp\
    420420                                        ./parallel/OutputDiagnostic.cpp\
    421                                         ./parallel/WriteLockFile.cpp
     421                                        ./parallel/WriteLockFile.cpp\
     422                                        ./parallel/control.cpp\
     423                                        ./parallel/OutputControl.cpp\
     424                                        ./parallel/objectivefunctionC.cpp\
     425                                        ./parallel/GradJCompute.cpp
    422426
    423427libpISSM_a_CXXFLAGS = -fPIC -D_PARALLEL_   -D_C_
    … …  
    426430bin_PROGRAMS =
    427431else
    428 bin_PROGRAMS = diagnostic.exe
    429 #control.exe thermalsteady.exe
     432bin_PROGRAMS = diagnostic.exe control.exe
     433#thermalsteady.exe
    430434endif
    431435

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

    r1 r46  
    6464        double tria_q;
    6565        int    tria_shelf;
    66         double tria_fit;
    6766        double tria_meanvel;/*!scaling ratio for velocities*/
    6867        double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
    69         int    tria_acceleration;
    7068
    7169        /*matice constructor input: */
    … …  
    9088        int penta_onbed;
    9189        int penta_onsurface;
    92         double penta_fit;
    9390        double penta_meanvel;/*!scaling ratio for velocities*/
    9491        double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
    95         int penta_acceleration;
    9692        int penta_collapse;
    9793        double penta_melting[6];
    … …  
    151147        int     grid_id;
    152148
     149        /*Get analysis_type: */
     150        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
    153151
    154152        /*First create the elements, nodes and material properties: */
    … …  
    156154        nodes     = new DataSet(NodesEnum());
    157155        materials = new DataSet(MaterialsEnum());
    158 
    159156
    160157        /*Width of elements: */
    … …  
    226223                ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
    227224                ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
    228                 ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat");
    229225                ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
    230226                ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
    231        
    232227               
    233228                for (i=0;i<model->numberofelements;i++){
    … …  
    275270                        tria_shelf=(int)*(model->elementoniceshelf+i);
    276271
    277                         /*type of fitting? : */
    278                         tria_fit=model->fit[0];
    279272                        tria_meanvel=model->meanvel;
    280273                        tria_epsvel=model->epsvel;
    281274
    282                         /*diverse:*/
    283                         tria_acceleration=0;
    284 
    285275                        /*Create tria element using its constructor:*/
    286                         tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_fit, tria_meanvel, tria_epsvel, tria_acceleration);
     276                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel);
    287277
    288278                        /*Add tria element to elements dataset: */
    … …  
    333323                xfree((void**)&model->q);
    334324                xfree((void**)&model->elementoniceshelf);
    335                 xfree((void**)&model->fit);
    336325                xfree((void**)&model->B);
    337326                xfree((void**)&model->n);*/
    … …  
    348337                ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
    349338                ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
    350                 ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat");
    351339                ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
    352340                ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
    … …  
    384372                        /*diverse: */
    385373                        penta_shelf=(int)*(model->elementoniceshelf+i);
    386                         penta_fit=(int)model->fit[0];
    387374                        penta_onbed=(int)*(model->elementonbed+i);
    388375                        penta_onsurface=(int)*(model->elementonsurface+i);
    389376                        penta_meanvel=model->meanvel;
    390377                        penta_epsvel=model->epsvel;
    391                         penta_acceleration=0;
    392378
    393379                        if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base.
    … …  
    401387                        /*Create Penta using its constructor:*/
    402388                        penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
    403                                         penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_fit,penta_meanvel,penta_epsvel,
    404                                         penta_acceleration,penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
     389                                        penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
     390                                        penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
    405391                                        penta_thermal_steadystate);
    406392
    … …  
    446432
    447433                /*Free data: */
    448                 /*xfree((void**)&model->elements);
     434                xfree((void**)&model->elements);
    449435                xfree((void**)&model->thickness);
    450436                xfree((void**)&model->surface);
    … …  
    454440                xfree((void**)&model->q);
    455441                xfree((void**)&model->elementoniceshelf);
    456                 xfree((void**)&model->fit);
    457442                xfree((void**)&model->elementonbed);
    458443                xfree((void**)&model->elementonsurface);
    459444                xfree((void**)&model->elements_type);
    460445                xfree((void**)&model->n);
    461                 xfree((void**)&model->B);*/
     446                xfree((void**)&model->B);
    462447
    463448        } //if (strcmp(meshtype,"2d")==0)
    … …  
    555540        ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
    556541
    557         /*Get analysis_type: */
    558         analysis_type=AnalysisTypeAsEnum(model->analysis_type);
    559 
     542       
    560543        /*Get number of dofs per node: */
    561544        DistributeNumDofs(&node_numdofs,analysis_type);

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

    r8 r46  
    2121};
    2222#else
     23
     24#include "./FemModel.h"
     25#include "./ParameterInputs.h"
     26
    2327struct OptArgs{
    24 
     28        FemModel* femmodel;
     29        double* p_g;
     30        double* u_g_obs;
     31        double* grad_g;
     32        ParameterInputs* inputs;
     33        int n;
    2534};
    2635#endif

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

    r8 r46  
    1616
    1717#endif
    18 
    19 
    20                                                    

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

    r1 r46  
    1919}
    2020Penta::Penta( int penta_id, int penta_mid, int penta_mparid, int penta_node_ids[6], double penta_h[6], double penta_s[6], double penta_b[6], double penta_k[6], int penta_friction_type,
    21                                 double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_fit, double penta_meanvel,double penta_epsvel,
    22                                 int penta_acceleration, int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6],
     21                                double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel,
     22                                int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6],
    2323                                int penta_artdiff, int penta_thermal_steadystate){
    2424       
    … …  
    4646        onbed = penta_onbed;
    4747        onsurface = penta_onsurface;
    48         fit = penta_fit;
    4948        meanvel = penta_meanvel;
    5049        epsvel = penta_epsvel;
    51         acceleration = penta_acceleration;
    5250        collapse = penta_collapse;
    5351        artdiff = penta_artdiff;
    … …  
    7977        printf("   onbed: %i\n",onbed);
    8078        printf("   onsurface: %i\n",onsurface);
    81         printf("   fit: %g\n",fit);
    8279        printf("   meanvel: %g\n",meanvel);
    8380        printf("   epsvel: %g\n",epsvel);
    84         printf("   acceleration: %i\n",acceleration);
    8581        printf("   collapse: %i\n",collapse);
    8682       
    … …  
    122118        memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
    123119        memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
    124         memcpy(marshalled_dataset,&fit,sizeof(fit));marshalled_dataset+=sizeof(fit);
    125120        memcpy(marshalled_dataset,&meanvel,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
    126121        memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
    127         memcpy(marshalled_dataset,&acceleration,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration);
    128122        memcpy(marshalled_dataset,&collapse,sizeof(collapse));marshalled_dataset+=sizeof(collapse);
    129123        memcpy(marshalled_dataset,&melting,sizeof(melting));marshalled_dataset+=sizeof(melting);
    … …  
    139133int   Penta::MarshallSize(){
    140134
    141         return sizeof(id)+sizeof(mid)+sizeof(mparid)+sizeof(node_ids)+sizeof(h)+sizeof(s)+sizeof(b)+sizeof(k)+sizeof(friction_type)+sizeof(p)+sizeof(q)+sizeof(shelf)+sizeof(onbed)+sizeof(onsurface)+sizeof(fit)+sizeof(meanvel)+sizeof(epsvel)+sizeof(acceleration)+sizeof(collapse)+sizeof(melting)+sizeof(accumulation)+sizeof(geothermalflux)+sizeof(artdiff)+sizeof(thermal_steadystate) +sizeof(int); //sizeof(int) for enum type
     135        return sizeof(id)+sizeof(mid)+sizeof(mparid)+sizeof(node_ids)+sizeof(h)+sizeof(s)+sizeof(b)+sizeof(k)+sizeof(friction_type)+sizeof(p)+sizeof(q)+sizeof(shelf)+sizeof(onbed)+sizeof(onsurface)+sizeof(meanvel)+sizeof(epsvel)+sizeof(collapse)+sizeof(melting)+sizeof(accumulation)+sizeof(geothermalflux)+sizeof(artdiff)+sizeof(thermal_steadystate) +sizeof(int); //sizeof(int) for enum type
    142136}
    143137
    … …  
    170164        memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
    171165        memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface);
    172         memcpy(&fit,marshalled_dataset,sizeof(fit));marshalled_dataset+=sizeof(fit);
    173166        memcpy(&meanvel,marshalled_dataset,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
    174167        memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
    175         memcpy(&acceleration,marshalled_dataset,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration);
    176168        memcpy(&collapse,marshalled_dataset,sizeof(collapse));marshalled_dataset+=sizeof(collapse);
    177169        memcpy(&melting,marshalled_dataset,sizeof(melting));marshalled_dataset+=sizeof(melting);

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

    r1 r46  
    4242                int    onbed;
    4343                int    onsurface;
    44                 double fit;
    4544                double meanvel;/*!scaling ratio for velocities*/
    4645                double epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
    47                 int    acceleration;
    4846                int    collapse;
    4947                double melting[6];
    … …  
    5856                Penta();
    5957                Penta( int id, int mid, int mparid, int g[6], double h[6], double s[6], double b[6], double k[6], int    friction_type,
    60                                 double p, double q, int    shelf, int    onbed, int    onsurface, double fit, double meanvel,double epsvel,
    61                                 int    acceleration, int    collapse, double melting[6], double accumulation[6], double geothermalflux[6],
     58                                double p, double q, int    shelf, int    onbed, int    onsurface, double meanvel,double epsvel,
     59                                int    collapse, double melting[6], double accumulation[6], double geothermalflux[6],
    6260                                int    artdiff, int    thermal_steadystate);
    6361                ~Penta();

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

    r1 r46  
    3333
    3434Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],
    35                                 int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_fit,double tria_meanvel,double tria_epsvel,int tria_acceleration){
     35                                int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel){
    3636       
    3737        int i;
    … …  
    5757        q=tria_q;
    5858        shelf=tria_shelf;
    59         fit=tria_fit;
    6059        meanvel=tria_meanvel;
    6160        epsvel=tria_epsvel;
    62         acceleration=tria_acceleration;
    6361        onbed=1;
    6462       
    … …  
    8785        printf("   q: %g\n",q);
    8886        printf("   shelf: %i\n",shelf);
    89         printf("   fit: %g\n",fit);
    9087        printf("   meanvel: %g\n",meanvel);
    9188        printf("   epsvel: %g\n",epsvel);
    92         printf("   acceleration: %i\n",acceleration);
    9389        printf("   onbed: %i\n",onbed);
    9490        printf("   nodes: \n");
    … …  
    136132        memcpy(marshalled_dataset,&q,sizeof(q));marshalled_dataset+=sizeof(q);
    137133        memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
    138         memcpy(marshalled_dataset,&fit,sizeof(fit));marshalled_dataset+=sizeof(fit);
    139134        memcpy(marshalled_dataset,&meanvel,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
    140135        memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
    141         memcpy(marshalled_dataset,&acceleration,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration);
    142136       
    143137        *pmarshalled_dataset=marshalled_dataset;
    … …  
    165159                +sizeof(q)
    166160                +sizeof(shelf)
    167                 +sizeof(fit)
    168161                +sizeof(meanvel)
    169162                +sizeof(epsvel)
    170                 +sizeof(acceleration)
    171163                +sizeof(int); //sizeof(int) for enum type
    172164}
    … …  
    206198        memcpy(&q,marshalled_dataset,sizeof(q));marshalled_dataset+=sizeof(q);
    207199        memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
    208         memcpy(&fit,marshalled_dataset,sizeof(fit));marshalled_dataset+=sizeof(fit);
    209200        memcpy(&meanvel,marshalled_dataset,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
    210201        memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
    211         memcpy(&acceleration,marshalled_dataset,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration);
    212202
    213203        /*nodes and materials are not pointing to correct objects anymore:*/
    … …  
    338328        double MOUNTAINKEXPONENT=10;
    339329
    340 
    341         /*First check that acceleration is not on: */
    342         if (acceleration){
    343                 /*do nothing: */
    344                 return;
    345         }
    346330
    347331        /*recover extra inputs from users, at current convergence iteration: */
    … …  
    650634        double  thickness;
    651635
    652         /*First check that acceleration is not on: */
    653         if (acceleration){
    654                 return;
    655         }
    656 
    657636        /* Get node coordinates and dof list: */
    658637        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    13131292                fit=*fit_param;
    13141293        }
     1294        else throw ErrorException(__FUNCT__," missing fit input parameter");
    13151295
    13161296        /*Initialize velocities: */
    … …  
    18441824                fit=*fit_param;
    18451825        }
     1826        else ErrorException(__FUNCT__," missing fit input parameter");
    18461827
    18471828        /*Initialize velocities: */

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

    r1 r46  
    4040                double q;
    4141                int    shelf;
    42                 double fit;
    4342                double meanvel;/*!scaling ratio for velocities*/
    4443                double epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
    45                 int    acceleration;
    4644                int    onbed;
    4745       
    … …  
    5048                Tria();
    5149                Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],
    52                                 int friction_type,double p,double q,int shelf,double fit,double meanvel,double epsvel,int acceleration);
     50                                int friction_type,double p,double q,int shelf,double meanvel,double epsvel);
    5351                ~Tria();
    5452

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

    r1 r46  
    6868        ConfigureObjectsx(elements, loads, nodes, materials);
    6969       
     70        _printf_("   process parameters:\n");
     71        ProcessParamsx( parameters, partition);
     72       
    7073        _printf_("   free ressources:\n");
    7174        DeleteModel(&model);
    … …  
    8689        femmodel->ys=ys;
    8790        femmodel->ys0=ys0;
     91
    8892       
    8993}

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

    r1 r46  
    1313void OutputDiagnostic(Vec u_g,Vec partition,char* filename,NodeSets* nodesets,char* analysis_type){
    1414
    15         /* error handling: */
    1615        int i;
    17 
    1816        extern int my_rank;
    1917       
    … …  
    2321        /* standard output: */
    2422        double* serial_partition=NULL;
    25         int     serial_partition_rows;
    2623        int    analysis_size;
    2724
    2825        double* serial_u_g=NULL;
    29         int     serial_u_g_rows;
    30         int     dummy=1;
     26        int     one=1;
    3127        int     gsize;
    3228
    3329        /*serialize outputs: */
    3430        VecShift(partition,1.0); //matlab indexing
    35         VecGetSize(partition,&serial_partition_rows);
    3631        VecToMPISerial(&serial_partition,partition);
    3732
    38         VecGetSize(u_g,&serial_u_g_rows);
    3933        VecToMPISerial(&serial_u_g,u_g);
    4034
    … …  
    5246
    5347                /*Write partition: */
    54                 WriteDataToDisk(serial_partition,&serial_partition_rows,&dummy,"Mat",fid);
     48                WriteDataToDisk(serial_partition,&one,&one,"Mat",fid);
    5549               
    5650                /*Write solution to disk: */
    57                 WriteDataToDisk(serial_u_g,&serial_u_g_rows,&dummy,"Mat",fid);
     51                WriteDataToDisk(serial_u_g,&gsize,&one,"Mat",fid);
    5852       
    5953                /*Close file: */

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

    r1 r46  
    1 /*
    2  * cielocontrol.c:
    3  */
     1/*!\file:  control.cpp
     2 * \brief: control solution
     3 */ 
    44
    5 #include "../include/cielo.h"
    6 #include "../modules.h"
     5#include "../issm.h"
    76#include "./parallel.h"
    87
    98#undef __FUNCT__
    10 #define __FUNCT__ "cielocontrol"
     9#define __FUNCT__ "control"
     10
     11#ifdef HAVE_CONFIG_H
     12        #include "config.h"
     13#else
     14#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     15#endif
    1116
    1217int main(int argc,char* *argv){
    13        
    14         int i,n;
    15        
     18
     19        int n;
     20        int i;
     21
    1622        /*I/O: */
    1723        FILE* fid=NULL;
    … …  
    2127        char* analysis_type="control";
    2228
    23         /*Finite element model: */
    24         #if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
     29        /*Intermediary: */
    2530        FemModel femmodel;
     31        Vec u_g=NULL;
     32        ParameterInputs* inputs=NULL;
     33        int waitonlock=0;
     34        double search_scalar;
     35        char* control_type=NULL;
     36        int   gsize;
     37        double* fit=NULL;
     38        double* optscal=NULL;
     39        double* u_g_obs=NULL;
     40        int*    maxiter=NULL;
     41        double  tolx;
     42        double*   p_g=NULL;
     43        Vec       grad_g=NULL;
     44        Vec       new_grad_g=NULL;
     45        Vec       grad_g_old=NULL;
     46        double*   grad_g_double=NULL;
     47        double    mincontrolconstraint;
     48        double    maxcontrolconstraint;
     49        int       nsteps;
     50        double*   J=NULL;
     51        OptArgs   optargs;
     52        OptPars   optpars;
     53
     54
     55        MODULEBOOT();
     56
     57        #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
     58        throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!");
    2659        #endif
    2760
    28         /*control : */
    29         #if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
    30         Vec*             u_g=NULL;
    31         #endif
    32         double*          search_vector=NULL;
    33         int              status;
    34         ParameterInputs* inputs=NULL;
    35         int              reloop;
     61        PetscInitialize(&argc,&argv,(char *)0,""); 
     62
     63        /*Size and rank: */
     64        MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); 
     65        MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
     66
     67        _printf_("recover , input file name and output file name:\n");
     68        inputfilename=argv[2];
     69        outputfilename=argv[3];
     70        lockname=argv[4];
     71
     72        /*Open handle to data on disk: */
     73        fid=fopen(inputfilename,"rb");
     74        if(fid==NULL) throw ErrorException(__FUNCT__,exprintf("%s%s%s","could not open file ",inputfilename," for binary reading"));
    3675       
    37         #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
    38                 _printf_("%s%s\n",__FUNCT__," error message: parallel executable was compiled without support of parallel libraries!");
    39                 return 1;
    40         #else
     76        _printf_("read and create finite element model:\n");
     77        CreateFemModel(&femmodel,fid,analysis_type);
     78
     79        /*Recover parameters used throughout the solution:*/
     80        femmodel.parameters->FindParam((void*)&nsteps,"nsteps");
     81        femmodel.parameters->FindParam((void*)&control_type,"control_type");
     82        femmodel.parameters->FindParam((void*)&fit,"fit");
     83        femmodel.parameters->FindParam((void*)&optscal,"optscal");
     84        femmodel.parameters->FindParam((void*)&maxiter,"maxiter");
     85        femmodel.parameters->FindParam((void*)&tolx,"tolx");
     86        femmodel.parameters->FindParam((void*)&waitonlock,"waitonlock");
     87        femmodel.parameters->FindParam((void*)&mincontrolconstraint,"mincontrolconstraint");
     88        femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint");
     89        gsize=femmodel.nodes->NumberOfDofs();
    4190       
    42                 /*Initialize MPI environment: */
    43                 PetscInitialize(&argc,&argv,(char *)0,""); 
     91        femmodel.parameters->FindParam((void*)&p_g,"p_g");
     92        femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs");
    4493
    45                 /*Size and rank: */
    46                 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); 
    47                 MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
    4894
    49                 /*Some checks on size of cluster*/
    50                 if (num_procs<=1){
    51                         _printf_("\nSize of MPI COMM WORLD is 1, needs to be at least 2. Include more nodes\n");
    52                         PetscFinalize();
    53                         return 0;
     95        /*Start looping: */
     96        for(n=0;n<nsteps;n++){
     97                       
     98                _printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
     99
     100                //initialize inputs, ie parameters on which we invert.
     101                DeleteParameterInputs(&inputs); inputs=NewParameterInputs();
     102                ParameterInputsAddFromMat(inputs,p_g,gsize,control_type);
     103                ParameterInputsAddFromDouble(inputs,fit[n],"fit");
     104
     105                _printf_("%s\n","      computing gradJ...");
     106                grad_g=GradJCompute(inputs,&femmodel,u_g_obs);
     107                _printf_("%s\n","      done.");
     108                       
     109               
     110                _printf_("%s\n","      normalizing directions...");
     111                Orthx(&new_grad_g,grad_g,grad_g_old);
     112                VecFree(&grad_g); grad_g=new_grad_g;
     113                VecFree(&grad_g_old); grad_g_old=grad_g;
     114                VecToMPISerial(&grad_g_double,grad_g);
     115                _printf_("%s\n","      done.");
     116
     117       
     118                _printf_("%s\n","      optimizing along gradient direction...");
     119                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;
     120                optpars.xmin=-1; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=maxiter[n];
     121                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
     122                _printf_("%s\n","      done.");
     123       
     124       
     125               
     126                _printf_("%s\n","      updating parameter using optimized search scalar...");
     127                for(i=0;i<gsize;i++)p_g[i]=p_g[i]+search_scalar*optscal[n]*grad_g_double[i];
     128                _printf_("%s\n","      done.");
     129
     130                _printf_("%s\n","      constraining the new distribution...");   
     131                ControlConstrainx( p_g,gsize,mincontrolconstraint,maxcontrolconstraint,control_type);
     132                _printf_("%s\n","      done.");
     133       
     134                //some temporary saving
     135                if ((n%5)==0){
     136                        _printf_("%s\n","      saving temporary results...");
     137                        OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
     138                        _printf_("%s\n","      done.");
    54139                }
    55140
    56                 /*Recover dbdir, input file name and output file name: */
    57                 dbdir=argv[1];
    58                 inputfilename=argv[2];
    59                 outputfilename=argv[3];
    60                 lockname=argv[4];
     141                _printf_("%s\n","      value of misfit J after optimization #",n,": ",J[n]);
     142               
     143                /*some freeing:*/
     144                xfree((void**)&grad_g_double);
     145        }
    61146
    62                 /*Open handle to data on disk: */
    63                 fid=fopen(inputfilename,"rb");
    64                 if(fid==NULL){
    65                         _printf_("%s%s\n",__FUNCT__," error message: could not open file ",inputfilename," for binary reading");
    66                         return 0;
    67                 }
     147        /*Write results to disk: */
     148        _printf_("%s\n","       preparing final velocity solution...");
     149        DeleteParameterInputs(&inputs); inputs=NewParameterInputs();
     150        ParameterInputsAddFromMat(inputs,p_g,gsize,control_type);
     151        ParameterInputsAddFromDouble(inputs,fit[n],"fit");
     152       
     153        diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel);
     154       
     155        _printf_("%s\n","      saving final results...");
     156        OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets);
     157        _printf_("%s\n","      done.");
    68158
    69                 /*Read and create finite element model: */
    70                 if(!CreateFemModel(&femmodel,fid,analysis_type)){
    71                         _printf_("%s\n",__FUNCT__," error message: could not read finite element model!\n");
    72                         return 0;
    73                 }
     159        /*Write lock file if requested: */
     160        if (waitonlock){
     161                WriteLockFile(lockname);
     162        }
     163                       
     164        _printf_("closing MPI and Petsc\n");
     165        MPI_Barrier(MPI_COMM_WORLD);
    74166
    75                 /*Initialize inputs: */
    76                 inputs=NewParameterInputs();
    77 
    78                 /*For now, only one parameter is allowed: */
    79                 if (femmodel.workspaceparams->num_control_parameters>1){
    80                         _printf_("%s%s\n",__FUNCT__," error message: multiple control parameters not implemented yet!");
    81                         PetscFinalize();
    82                         return 1;
    83                 }
    84 
    85                 /*Initialize search vector: */
    86                 search_vector=xmalloc(femmodel.workspaceparams->num_control_parameters*sizeof(double));
    87 
    88                 /*Start looping: */
    89                 for(n=0;n<femmodel.workspaceparams->nsteps;n++){
    90                        
    91                         _printf_("\n%s%i%s%i\n","   control method step ",(n+1),"/",femmodel.workspaceparams->nsteps);
    92 
    93                         //initialize inputs, ie parameters on which we invert.
    94                         DeleteParameterInputs(&inputs); inputs=NewParameterInputs();
    95                        
    96                         for(i=0;i<femmodel.workspaceparams->num_control_parameters;i++){
    97                                 char* control_type=femmodel.workspaceparams->control_types[i];
    98                                 ParameterInputsAddFromMat(inputs,WorkspaceParamsGetParameter(femmodel.workspaceparams,control_type),femmodel.workspaceparams->gsize,control_type);
    99                         }
    100                         ParameterInputsAddFromDouble(inputs,femmodel.workspaceparams->fit[n],"fit");
    101 
    102                         //update velocity:
    103                         _printf_("%s\n","      updating velocity...");
    104                         VecFree(&u_g); cielodiagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel);
    105                         _printf_("%s\n","      done.");
    106 
    107                         //call gradJ module
    108                         _printf_("%s\n","      computing gradJ...");
    109                         GradJCompute(inputs,&femmodel);
    110                         _printf_("%s\n","      done.");
    111                        
    112                         //normalize directions
    113                         _printf_("%s\n","      normalizing directions...");
    114                         GradJOrth(femmodel.workspaceparams);
    115                         _printf_("%s\n","      done.");
     167        /*Close MPI libraries: */
     168        PetscFinalize();
    116169
    117170
    118                         //search along the direction for the parameter vector that minimizes the misfit J
    119                         _printf_("%s\n","      optimizing...");
    120                         status=GradJSearch(search_vector,&femmodel,n);
    121                         _printf_("%s\n","      done.");
    122                        
    123                         _printf_("%s%g\n","      misfit: ",femmodel.workspaceparams->J[n]);
    124 
    125                         //Check the GradJSearch actually resulting in improving the misfit. Otherwise, reloop.
    126                         reloop=GradJCheck(femmodel.workspaceparams,n,status);
    127                         if (reloop){
    128                                 printf("reloop\n",reloop);
    129                                 n=n-1;
    130                                 continue;
    131                         }
    132 
    133                         //update parameters with new optimized values
    134                         _printf_("%s\n","      updating parameters...");
    135                         ParameterUpdate(search_vector,n,femmodel.workspaceparams,femmodel.batchparams);
    136                         _printf_("%s\n","      done.");
    137 
    138                         //temporary saving
    139                         if ((n%1)==0){
    140                                 _printf_("%s\n","      saving temporary results...");
    141                                 OutputControl(femmodel.workspaceparams,femmodel.batchparams,u_g,femmodel.tpartition,outputfilename,"control");
    142                                 _printf_("%s\n","      done.");
    143                         }
    144 
    145                 }
    146                
    147                 /*Write results to disk: */
    148                 _printf_("%s\n","      saving final results...");
    149                 OutputControl(femmodel.workspaceparams,femmodel.batchparams,u_g,femmodel.tpartition,outputfilename,"control");
    150                 _printf_("%s\n","      done.");
    151 
    152                 /*Write lock file if requested: */
    153                 if (femmodel.batchparams->waitonlock){
    154                         WriteLockFile(lockname);
    155                 }
    156                        
    157                 cleanup_and_return:
    158                 _printf_("exiting.\n");
    159                
    160                 /*Synchronize everyone before exiting: */
    161                 MPI_Barrier(MPI_COMM_WORLD);
    162 
    163                 /*Close MPI libraries: */
    164                 PetscFinalize();
    165                
    166                 return 0; //unix success return;
    167         #endif
     171        /*end module: */
     172        MODULEEND();
     173       
     174        return 0; //unix success return;
    168175}

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

    r1 r46  
    5959
    6060        _printf_("call computational core:\n");
    61         diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,femmodel);
     61        diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel);
    6262
    6363        _printf_("write results to disk:\n");

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

    r1 r46  
    1111#include "../issm.h"
    1212
    13 void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, ParameterInputs* inputs,FemModel fem){
     13void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, ParameterInputs* inputs,FemModel* fem){
    1414
    1515
    … …  
    4545        /*Recover parameters: */
    4646        kflag=1; pflag=1;
    47         fem.parameters->FindParam((void*)&connectivity,"connectivity");
    48         fem.parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
    49         fem.parameters->FindParam((void*)&solver_string,"solverstring");
    50         fem.parameters->FindParam((void*)&eps_rel,"eps_rel");
    51         fem.parameters->FindParam((void*)&debug,"debug");
     47        fem->parameters->FindParam((void*)&connectivity,"connectivity");
     48        fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
     49        fem->parameters->FindParam((void*)&solver_string,"solverstring");
     50        fem->parameters->FindParam((void*)&eps_rel,"eps_rel");
     51        fem->parameters->FindParam((void*)&debug,"debug");
    5252
    53         fem.parameters->FindParam((void*)&analysis_type_string,"analysis_type");
     53        fem->parameters->FindParam((void*)&analysis_type_string,"analysis_type");
    5454    analysis_type=AnalysisTypeAsEnum(analysis_type_string);
    5555
    5656        /*Copy loads for backup: */
    57         loads=fem.loads->Copy();
     57        loads=fem->loads->Copy();
    5858
    5959        count=1;
    … …  
    7171
    7272                /*Update parameters: */
    73                 UpdateFromInputsx(fem.elements,fem.nodes,loads, fem.materials,inputs);
     73                UpdateFromInputsx(fem->elements,fem->nodes,loads, fem->materials,inputs);
    7474
    7575                if (debug) _printf_("   Generating matrices\n");
    7676                //*Generate system matrices
    77                 SystemMatricesx(&Kgg, &pg,fem.elements,fem.nodes,loads,fem.materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
     77                SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
    7878
    7979                if (debug) _printf_("   Generating penalty matrices\n");
    8080                //*Generate penalty system matrices
    81                 PenaltySystemMatricesx(Kgg, pg,fem.elements,fem.nodes,loads,fem.materials,kflag,pflag,inputs,analysis_type);
     81                PenaltySystemMatricesx(Kgg, pg,fem->elements,fem->nodes,loads,fem->materials,kflag,pflag,inputs,analysis_type);
    8282
    8383                if (debug) _printf_("   reducing matrix from g to f set\n");
    8484                /*!Reduce matrix from g to f size:*/
    85                 Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem.Gmn,fem.nodesets);
     85                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
    8686
    8787                /*Free ressources: */
    … …  
    9090                if (debug) _printf_("   reducing load from g to f set\n");
    9191                /*!Reduce load from g to f size: */
    92                 Reduceloadfromgtofx(&pf, pg, fem.Gmn, Kfs, fem.ys, fem.nodesets);
     92                Reduceloadfromgtofx(&pf, pg, fem->Gmn, Kfs, fem->ys, fem->nodesets);
    9393
    9494                //no need for pg and Kfs anymore
    … …  
    106106                if (debug) _printf_("   merging solution from f to g set\n");
    107107                //Merge back to g set
    108                 Mergesolutionfromftogx(&ug, uf,fem.Gmn,fem.ys,fem.nodesets);
     108                Mergesolutionfromftogx(&ug, uf,fem->Gmn,fem->ys,fem->nodesets);
    109109
    110110                //Deal with penalty loads
    … …  
    112112                ParameterInputsAddFromVec(inputs,ug,"velocity");
    113113               
    114                 PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem.elements,fem.nodes,loads,fem.materials,inputs,analysis_type);
     114                PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type);
    115115
    116116                //Figure out if convergence is reached.
    … …  
    143143                kflag=1; pflag=0; //stiffness generation only
    144144       
    145                 SystemMatricesx(&Kgg, &pg,fem.elements,fem.nodes,fem.loads,fem.materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
     145                SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->loads,fem->materials,kflag,pflag,connectivity,numberofdofspernode,inputs,analysis_type);
    146146       
    147                 Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem.Gmn,fem.nodesets);
     147                Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
    148148               
    149149                MatFree(&Kgg);VecFree(&pg);

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

    r1 r46  
    88#include "../objects/objects.h"
    99
    10 int GradJCompute(ParameterInputs* inputs,FemModel* femmodel);
     10Vec GradJCompute(ParameterInputs* inputs,FemModel* femmodel,double* u_g_obs);
    1111
    12 void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, ParameterInputs* inputs,FemModel fem);
     12void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, ParameterInputs* inputs,FemModel* fem);
    1313int cielothermal_core(Vec** pt_g,ParameterInputs* inputs,FemModel* femmodel);
    1414
    … …  
    1919int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*,ParameterInputs*),FemModel* femmodel);
    2020       
    21 double objectivefunctionC(double* search_vector,double fit,double optscal,FemModel*  femmodel,ParameterInputs* inputs);
     21double objectivefunctionC(double search_scalar,OptArgs* optargs);
    2222
    2323int GradJSearch(double* search_vector,FemModel* femmodel,int step);
    … …  
    2727void OutputDiagnostic(Vec u_g,Vec tpartition,char* filename,NodeSets* nodesets,char* analysis_type);
    2828int OutputThermal(Vec* t_g,Vec* tpartition,char* filename,char* analysis_type);
    29 //int OutputControl(WorkspaceParams* workspaceparams,BatchParams* batchparams,Vec* u_g,Vec* tpartition, char* filename,char* analysis_type);
     29void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets);
    3030void WriteLockFile(char* filename);
    3131
Note: See TracChangeset for help on using the changeset viewer.