Changeset 377


Ignore:
Timestamp:
05/13/09 09:15:49 (16 years ago)
Author:
Eric.Larour
Message:

New model processor for stokes. Still to be completed.
Serial and parallel control now act the same way with respect to velocity inputs.

Location:
issm/trunk/src
Files:
23 edited

Legend:

Unmodified
Added
Removed

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

    r300 r377  
    154154        param= new Param(count,"min_thermal_constraints",INTEGER);
    155155        param->SetInteger(model->min_thermal_constraints);
     156        parameters->AddObject(param);
     157
     158        /*reconditioning_number: */
     159        count++;
     160        param= new Param(count,"reconditioning_number",DOUBLE);
     161        param->SetDouble(model->reconditioning_number);
    156162        parameters->AddObject(param);
    157163

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

    r347 r377  
    9999        int penta_thermal_steadystate;
    100100        double penta_viscosity_overshoot;
     101        double penta_reconditioning_number;
    101102
    102103        /*matpar constructor input: */
    … …  
    413414                                        penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
    414415                                        penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
    415                                         penta_thermal_steadystate,penta_viscosity_overshoot);
     416                                        penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number);
    416417
    417418                        /*Add penta element to elements dataset: */

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp

    r300 r377  
    1616void    CreateConstraintsDiagnosticStokes(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
    1717
     18        int i;
     19        DataSet* constraints = NULL;
     20        Spc*    spc  = NULL;
    1821
    19         DataSet* constraints = NULL;
     22        /*spc intermediary data: */
     23        int spc_sid;
     24        int spc_node;
     25        int spc_dof;
     26        double spc_value;
     27        int count;
     28       
     29        double* gridonstokes=NULL;
     30
    2031
    2132        /*Create constraints: */
    … …  
    2435        /*Now, is the flag ishutter on? otherwise, do nothing: */
    2536        if (!model->isstokes)goto cleanup_and_return;
     37
     38        /*Fetch data: */
     39        ModelFetchData((void**)&gridonstokes,NULL,NULL,model_handle,"gridonstokes","Matrix","Mat");
     40
     41        count=0;
     42        /*Create spcs from x,y,z, as well as the spc values on those spcs: */
     43        for (i=0;i<model->numberofnodes;i++){
     44        #ifdef _PARALLEL_
     45        /*keep only this partition's nodes:*/
     46        if((model->my_grids[i]==1)){
     47        #endif
     48
     49                if ((int)!gridonstokes[i]){
    2650       
    27         cleanup_and_return:     
     51                        /*This grid will see its vx,vy and vz dofs spc'd to pattyn velocities: */
     52                        spc_sid=count;
     53                        spc_node=i+1;
     54                        spc_dof=1; //we enforce first x translation degree of freedom
     55                        spc_value=0;
     56
     57                        spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
     58                        constraints->AddObject(spc);
     59                        count++;
     60
     61                        spc_sid=count;
     62                        spc_node=i+1;
     63                        spc_dof=2; //we enforce first y translation degree of freedom
     64                        spc_value=0;
     65                       
     66                        spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
     67                        constraints->AddObject(spc);
     68                        count++;
     69                       
     70                        spc_sid=count;
     71                        spc_node=i+1;
     72                        spc_dof=3; //we enforce first y translation degree of freedom
     73                        spc_value=0;
     74                       
     75                        spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
     76                        constraints->AddObject(spc);
     77                        count++;
     78
     79                }
     80
     81        #ifdef _PARALLEL_
     82        } //if((my_grids[i]==1))
     83        #endif
     84        }
     85
     86        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     87         * datasets, it will not be redone: */
     88        constraints->Presort();
     89
     90        cleanup_and_return:
     91        /*Free data: */
     92        xfree((void**)&gridonstokes);
    2893       
    2994        /*Assign output pointer: */
    3095        *pconstraints=constraints;
    31 }
     96}       

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp

    r300 r377  
    1818
    1919
     20       
     21        /*output: int* epart, int* my_grids, double* my_bordergrids*/
     22
     23
     24        int i,j,k,n;
     25        extern int my_rank;
     26        extern int num_procs;
     27
    2028        /*DataSets: */
    2129        DataSet*    elements  = NULL;
    … …  
    2331        DataSet*    materials = NULL;
    2432       
     33        /*Objects: */
     34        Node*       node   = NULL;
     35        Penta*      penta = NULL;
     36        Matice*     matice  = NULL;
     37        Matpar*     matpar  = NULL;
     38
     39        int         analysis_type;
     40       
     41        /*output: */
     42        int* epart=NULL; //element partitioning.
     43        int* npart=NULL; //node partitioning.
     44        int* my_grids=NULL;
     45        double* my_bordergrids=NULL;
     46
     47
     48        /*intermediary: */
     49        int elements_width; //size of elements
     50        double B_avg;
     51                       
     52        /*matice constructor input: */
     53        int    matice_mid;
     54        double matice_B;
     55        double matice_n;
     56       
     57        /*penta constructor input: */
     58
     59        int penta_id;
     60        int penta_mid;
     61        int penta_mparid;
     62        int penta_g[6];
     63        double penta_h[6];
     64        double penta_s[6];
     65        double penta_b[6];
     66        double penta_k[6];
     67        int penta_friction_type;
     68        double penta_p;
     69        double penta_q;
     70        int penta_shelf;
     71        int penta_onbed;
     72        int penta_onsurface;
     73        double penta_meanvel;/*!scaling ratio for velocities*/
     74        double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
     75        int penta_collapse;
     76        double penta_melting[6];
     77        double penta_accumulation[6];
     78        double penta_geothermalflux[6];
     79        int penta_artdiff;
     80        int penta_thermal_steadystate;
     81        double penta_viscosity_overshoot;
     82        double penta_reconditioning_number;
     83
     84        /*matpar constructor input: */
     85        int     matpar_mid;
     86        double matpar_rho_ice;
     87        double matpar_rho_water;
     88        double matpar_heatcapacity;
     89        double matpar_thermalconductivity;
     90        double matpar_latentheat;
     91        double matpar_beta;
     92        double matpar_meltingpoint;
     93        double matpar_mixed_layer_capacity;
     94        double matpar_thermal_exchange_velocity;
     95        double matpar_g;
     96
     97        /* node constructor input: */
     98        int node_id;
     99        int node_partitionborder=0;
     100        double node_x[3];
     101        int node_onbed;
     102        int node_onsurface;
     103        int node_upper_node_id;
     104        int node_numdofs;
     105
     106               
     107        #ifdef _PARALLEL_
     108        /*Metis partitioning: */
     109        int  range;
     110        Vec  gridborder;
     111        int  my_numgrids;
     112        int* all_numgrids=NULL;
     113        int  gridcount;
     114        int  count;
     115        #endif
     116        int  first_grid_index;
     117
     118        /*Rifts:*/
     119        int* riftsnumpenaltypairs;
     120        double** riftspenaltypairs;
     121        int* riftsfill;
     122        double* riftsfriction;
     123        double* riftpenaltypairs=NULL;
     124        int el1,el2;
     125         
    25126        /*First create the elements, nodes and material properties: */
    26127        elements  = new DataSet(ElementsEnum());
    … …  
    31132        if (!model->isstokes)goto cleanup_and_return;
    32133
     134        /*Get analysis_type: */
     135        analysis_type=AnalysisTypeAsEnum(model->analysis_type);
     136
     137        /*Width of elements: */
     138        if(strcmp(model->meshtype,"2d")==0){
     139                elements_width=3; //tria elements
     140        }
     141        else{
     142                elements_width=6; //penta elements
     143        }
     144
     145
     146        #ifdef _PARALLEL_
     147        /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
     148        if(strcmp(model->meshtype,"2d")==0){
     149                /*load elements: */
     150                ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     151        }
     152        else{
     153                /*load elements2d: */
     154                ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
     155        }
     156
     157
     158        MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
     159
     160        /*Free elements and elements2d: */
     161        xfree((void**)&model->elements);
     162        xfree((void**)&model->elements2d);
     163
     164
     165        /*Deal with rifts, they have to be included into one partition only, not several: */
     166        FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts);
     167       
     168        for(i=0;i<model->numrifts;i++){
     169                riftpenaltypairs=model->riftspenaltypairs[i];
     170                for(j=0;j<model->riftsnumpenaltypairs[i];j++){
     171                        el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing
     172                        el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing
     173                        epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
     174                }
     175        }
     176        /*Free rifts: */
     177        xfree((void**)&riftsnumpenaltypairs);
     178        for(i=0;i<model->numrifts;i++){
     179                double* temp=riftspenaltypairs[i];
     180                xfree((void**)&temp);
     181        }
     182        xfree((void**)&riftspenaltypairs);
     183        xfree((void**)&riftsfill);
     184        xfree((void**)&riftsfriction);
     185
     186        /*Used later on: */
     187        my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int));
     188        #endif
     189
     190
     191        /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
     192
     193        /*2d mesh: */
     194        if (strcmp(model->meshtype,"2d")==0){
     195
     196                throw ErrorException(__FUNCT__," stokes elements only supported in 3d!");
     197
     198        }
     199        else{ //        if (strcmp(meshtype,"2d")==0)
     200
     201                /*Fetch data needed: */
     202                ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     203                ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
     204                ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
     205                ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
     206                ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
     207                ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat");
     208                ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat");
     209                ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat");
     210                ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat");
     211                ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat");
     212                ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
     213                ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
     214                ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
     215                ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat");
     216                ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat");
     217               
     218                for (i=0;i<model->numberofelements;i++){
     219                #ifdef _PARALLEL_
     220                /*We are using our element partition to decide which elements will be created on this node: */
     221                if(my_rank==epart[i]){
     222                #endif
     223
     224                       
     225                        /*name and id: */
     226                        penta_id=i+1; //matlab indexing.
     227                        penta_mid=i+1; //refers to the corresponding material property card
     228                        penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card
     229
     230                        /*vertices,thickness,surface,bed and drag: */
     231                        for(j=0;j<6;j++){
     232                                penta_g[j]=(int)*(model->elements+elements_width*i+j);
     233                                penta_h[j]=*(model->thickness+    ((int)*(model->elements+elements_width*i+j)-1));
     234                                penta_s[j]=*(model->surface+    ((int)*(model->elements+elements_width*i+j)-1));
     235                                penta_b[j]=*(model->bed+    ((int)*(model->elements+elements_width*i+j)-1));
     236                                penta_k[j]=*(model->drag+        ((int)*(model->elements+elements_width*i+j)-1));
     237                                penta_melting[j]=*(model->melting+        ((int)*(model->elements+elements_width*i+j)-1));
     238                                penta_accumulation[j]=*(model->accumulation+        ((int)*(model->elements+elements_width*i+j)-1));
     239                        }
     240
     241                        /*basal drag:*/
     242                        penta_friction_type=(int)model->drag_type;
     243       
     244                        penta_p=model->p[i];
     245                        penta_q=model->q[i];
     246
     247                        /*diverse: */
     248                        penta_shelf=(int)*(model->elementoniceshelf+i);
     249                        penta_onbed=(int)*(model->elementonbed+i);
     250                        penta_onsurface=(int)*(model->elementonsurface+i);
     251                        penta_meanvel=model->meanvel;
     252                        penta_epsvel=model->epsvel;
     253                       
     254                        /*viscosity_overshoot*/
     255                        penta_viscosity_overshoot=model->viscosity_overshoot;
     256
     257                        /*reconditioning_number: */
     258                        penta_reconditioning_number=model->reconditioning_number;
     259
     260                       
     261                        if (*(model->elements_type+2*i+1)==StokesEnum()){
     262                       
     263                                /*Create Penta using its constructor:*/
     264                                penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
     265                                                penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
     266                                                penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
     267                                                penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number);
     268
     269                                /*Add penta element to elements dataset: */
     270                                elements->AddObject(penta);
     271                        }
     272       
     273
     274                        /*Deal with material:*/
     275                        matice_mid=i+1; //same as the material id from the geom2 elements.
     276                        /*Average B over 6 element grids: */
     277                        B_avg=0;
     278                        for(j=0;j<6;j++){
     279                                B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1));
     280                        }
     281                        B_avg=B_avg/6;
     282                        matice_B= B_avg;
     283                        matice_n=(double)*(model->n+i);
     284                       
     285                        if (*(model->elements_type+2*i+1)==StokesEnum()){
     286       
     287                                /*Create matice using its constructor:*/
     288                                matice= new Matice(matice_mid,matice_B,matice_n);
     289               
     290                                /*Add matice element to materials dataset: */
     291                                materials->AddObject(matice);
     292                        }
     293                       
     294                        #ifdef _PARALLEL_
     295                        /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use
     296                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
     297                         into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids
     298                         will hold which grids belong to this partition*/
     299                        my_grids[(int)*(model->elements+elements_width*i+0)-1]=1;
     300                        my_grids[(int)*(model->elements+elements_width*i+1)-1]=1;
     301                        my_grids[(int)*(model->elements+elements_width*i+2)-1]=1;
     302                        my_grids[(int)*(model->elements+elements_width*i+3)-1]=1;
     303                        my_grids[(int)*(model->elements+elements_width*i+4)-1]=1;
     304                        my_grids[(int)*(model->elements+elements_width*i+5)-1]=1;
     305                        #endif
     306
     307                #ifdef _PARALLEL_
     308                }//if(my_rank==epart[i])
     309                #endif
     310
     311                }//for (i=0;i<numberofelements;i++)
     312
     313                /*Free data: */
     314                xfree((void**)&model->elements);
     315                xfree((void**)&model->thickness);
     316                xfree((void**)&model->surface);
     317                xfree((void**)&model->bed);
     318                xfree((void**)&model->drag);
     319                xfree((void**)&model->p);
     320                xfree((void**)&model->q);
     321                xfree((void**)&model->elementoniceshelf);
     322                xfree((void**)&model->elementonbed);
     323                xfree((void**)&model->elementonsurface);
     324                xfree((void**)&model->elements_type);
     325                xfree((void**)&model->n);
     326                xfree((void**)&model->B);
     327
     328        } //if (strcmp(meshtype,"2d")==0)
     329
     330        /*Add one constant material property to materials: */
     331        matpar_mid=model->numberofelements+1; //put it at the end of the materials
     332        matpar_g=model->g;
     333        matpar_rho_ice=model->rho_ice;
     334        matpar_rho_water=model->rho_water;
     335        matpar_thermalconductivity=model->thermalconductivity;
     336        matpar_heatcapacity=model->heatcapacity;
     337        matpar_latentheat=model->latentheat;
     338        matpar_beta=model->beta;
     339        matpar_meltingpoint=model->meltingpoint;
     340        matpar_mixed_layer_capacity=model->mixed_layer_capacity;
     341        matpar_thermal_exchange_velocity=model->thermal_exchange_velocity;
     342
     343        /*Create matpar object using its constructor: */
     344        matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
     345                        matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
     346                        matpar_thermal_exchange_velocity,matpar_g);
     347               
     348        /*Add to materials datset: */
     349        materials->AddObject(matpar);
     350       
     351        #ifdef _PARALLEL_
     352                /*From the element partitioning, we can determine which grids are on the inside of this cpu's
     353                 *element partition, and which are on its border with other nodes:*/
     354                gridborder=NewVec(model->numberofnodes);
     355
     356                for (i=0;i<model->numberofnodes;i++){
     357                        if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
     358                }
     359                VecAssemblyBegin(gridborder);
     360                VecAssemblyEnd(gridborder);
     361
     362                #ifdef _DEBUG_
     363                VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
     364                #endif
     365               
     366                VecToMPISerial(&my_bordergrids,gridborder);
     367
     368                #ifdef _DEBUG_
     369                if(my_rank==0){
     370                        for (i=0;i<model->numberofnodes;i++){
     371                                printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
     372                        }
     373                }
     374                #endif
     375        #endif
     376
     377        /*Partition penalties in 3d: */
     378        if(strcmp(model->meshtype,"3d")==0){
     379       
     380                /*Get penalties: */
     381                ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
     382
     383                if(model->numpenalties){
     384
     385                        model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int));
     386                        #ifdef _SERIAL_
     387                        for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=1;
     388                        #else
     389                        for(i=0;i<model->numpenalties;i++)model->penaltypartitioning[i]=-1;
     390
     391                        for(i=0;i<model->numpenalties;i++){
     392                                first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1);
     393                                if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
     394                                        /*All grids that are being penalised belong to this node's internal grid partition.:*/
     395                                        model->penaltypartitioning[i]=1;
     396                                }
     397                                if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
     398                                        model->penaltypartitioning[i]=0;
     399                                }
     400                        }
     401                        #endif
     402                }
     403
     404                /*Free penalties: */
     405                xfree((void**)&model->penalties);
     406        }
     407
     408        /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids.
     409         We can therefore determine  which grids are internal to this node's partition
     410         and which ones are shared with other nodes because they are on the border of this node's partition. Knowing
     411         that, go and create the grids*/
     412
     413        /*Create nodes from x,y,z, as well as the spc values on those grids: */
     414               
     415        /*First fetch data: */
     416        if (strcmp(model->meshtype,"3d")==0){
     417                ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
     418                ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
     419        }
     420        ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
     421        ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
     422        ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
     423        ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
     424        ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
     425        ModelFetchData((void**)&model->gridonstokes,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
     426        ModelFetchData((void**)&model->borderstokes,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
     427
     428
     429        /*Get number of dofs per node: */
     430        DistributeNumDofs(&node_numdofs,analysis_type);
     431
     432        for (i=0;i<model->numberofnodes;i++){
     433        #ifdef _PARALLEL_
     434        /*keep only this partition's nodes:*/
     435        if((my_grids[i]==1)){
     436        #endif
     437
     438                node_id=i+1; //matlab indexing
     439                       
     440               
     441               
     442                #ifdef _PARALLEL_
     443                if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
     444                        node_partitionborder=1;
     445                }
     446                else{
     447                        node_partitionborder=0;
     448                }
     449                #else
     450                        node_partitionborder=0;
     451                #endif
     452
     453
     454                node_x[0]=model->x[i];
     455                node_x[1]=model->y[i];
     456                node_x[2]=model->z[i];
     457
     458               
     459                node_onbed=(int)model->gridonbed[i];
     460                node_onsurface=(int)model->gridonsurface[i];   
     461                if (strcmp(model->meshtype,"3d")==0){
     462                        if (isnan(model->uppernodes[i])){
     463                                node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
     464                        }
     465                        else{
     466                                node_upper_node_id=(int)model->uppernodes[i];
     467                        }
     468                }
     469                else{
     470                        /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
     471                        node_upper_node_id=node_id;
     472                }
     473
     474                /*Create node using its constructor: */
     475                node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
     476
     477                /*set single point constraints.: */
     478                /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */
     479                if (model->gridonstokes[i]){
     480                        for(k=1;k<=node_numdofs;k++){
     481                                node->FreezeDof(k);
     482                        }
     483                }
     484                if (model->borderstokes[i]){
     485                        //freeze everything except pressure
     486                        node->FreezeDof(1);
     487                        node->FreezeDof(2);
     488                        node->FreezeDof(3);
     489                }
     490
     491
     492                /*Add node to nodes dataset: */
     493                nodes->AddObject(node);
     494
     495        #ifdef _PARALLEL_
     496        } //if((my_grids[i]==1))
     497        #endif
     498        }
     499
     500        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     501         * datasets, it will not be redone: */
     502        elements->Presort();
     503        nodes->Presort();
     504        materials->Presort();
     505
     506        /*Clean fetched data: */
     507        xfree((void**)&model->deadgrids);
     508        xfree((void**)&model->x);
     509        xfree((void**)&model->y);
     510        xfree((void**)&model->z);
     511        xfree((void**)&model->gridonbed);
     512        xfree((void**)&model->gridonsurface);
     513        xfree((void**)&model->uppernodes);
     514        xfree((void**)&model->gridonstokes);
     515        xfree((void**)&model->borderstokes);
     516               
    33517        cleanup_and_return:
    34518
    … …  
    37521        *pnodes=nodes;
    38522        *pmaterials=materials;
     523
     524        /*Keep partitioning information into model*/
     525        model->epart=epart;
     526        model->my_grids=my_grids;
     527        model->my_bordergrids=my_bordergrids;
     528
     529        /*Free ressources:*/
     530        #ifdef _PARALLEL_
     531        xfree((void**)&all_numgrids);
     532        xfree((void**)&npart);
     533        VecFree(&gridborder);
     534        #endif
     535
    39536}

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp

    r300 r377  
    1616void    CreateLoadsDiagnosticStokes(DataSet** ploads, Model* model,ConstDataHandle model_handle){
    1717
     18        int i,j,counter;
     19        int element;
     20
     21        extern int my_rank;
     22        extern int num_procs;
     23       
    1824        DataSet*    loads    = NULL;
    19 
     25        Icefront*   icefront = NULL;
     26        Penpair*    penpair  = NULL;
     27
     28        int segment_width;
     29        int element_type;
     30        int i1,i2,i3,i4;
     31        int i0,range;
     32        int grid1,grid2;
     33
     34        /*rift penpair: */
     35        double  normal[2];
     36        double  length;
     37        double  friction;
     38        int     fill;
     39               
     40        /*icefront intermediary data: */
     41        char icefront_type[ICEFRONTSTRING];
     42        int icefront_element_type;
     43        int     icefront_sid;
     44        int icefront_eid;
     45        int icefront_mparid;
     46        int     icefront_node_ids[MAX_ICEFRONT_GRIDS];
     47        double icefront_h[MAX_ICEFRONT_GRIDS];
     48        double  icefront_b[MAX_ICEFRONT_GRIDS];
     49
     50        /*penpair intermediary data: */
     51        int penpair_id;
     52        int penpair_node_ids[2];
     53        double penpair_penalty_offset;
     54        int penpair_numdofs;
     55        int penpair_dof;
     56        int penpair_penalty_lock;
     57        int penpair_element_ids[2];
     58        double penpair_friction;
     59        int penpair_fill;
     60        double penpair_normal[2];
     61        double penpair_length;
     62
     63        /*Rifts:*/
     64        int* riftsnumpenaltypairs=NULL;
     65        double** riftspenaltypairs=NULL;
     66        int* riftsfill=NULL;
     67        int* riftsfriction=NULL;
     68        double* riftpenaltypairs=NULL;
     69        int el1,el2;
     70       
    2071        /*Create loads: */
    2172        loads   = new DataSet(LoadsEnum());
    … …  
    2475        if (!model->isstokes)goto cleanup_and_return;
    2576
     77        /*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids
     78         * referenced by a certain load must belong to the cluster node): */
     79        ModelFetchData((void**)&model->segmentonneumann_diag,&model->numberofsegs_diag,NULL,model_handle,"segmentonneumann_diag","Matrix","Mat");
     80        ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
     81        ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
     82        ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
     83
     84        /*First load data:*/
     85        for (i=0;i<model->numberofsegs_diag;i++){
     86               
     87                if (strcmp(model->meshtype,"2d")==0){
     88                        segment_width=3;
     89                        element_type=TriaEnum();
     90                }
     91                else{
     92                        segment_width=5;
     93                        element_type=PentaEnum();
     94                }
     95
     96
     97                element=(int)(*(model->segmentonneumann_diag+segment_width*i+segment_width-1)-1); //element is in the last column
     98
     99                #ifdef _PARALLEL_
     100                if (model->epart[element]!=my_rank){
     101                        /*This load does not belong to this cluster node, as it references an element
     102                         *that does not belong to this node's partition. Drop this 'i':*/
     103                        continue;
     104                }
     105                #endif
     106       
     107                icefront_mparid=model->numberofelements+1; //matlab indexing
     108                icefront_sid=i+1; //matlab indexing
     109                icefront_eid=(int)*(model->segmentonneumann_diag+segment_width*i+segment_width-1); //matlab indexing
     110                icefront_element_type=element_type;
     111
     112                i1=(int)*(model->segmentonneumann_diag+segment_width*i+0);
     113                i2=(int)*(model->segmentonneumann_diag+segment_width*i+1);
     114                       
     115                icefront_node_ids[0]=i1;
     116                icefront_node_ids[1]=i2;
     117               
     118                icefront_h[0]=model->thickness[i1-1];
     119                icefront_h[1]=model->thickness[i2-1];
     120
     121                icefront_b[0]=model->bed[i1-1];
     122                icefront_b[1]=model->bed[i2-1];
     123               
     124                if ((int)*(model->elements_type+2*element+0)==MacAyealEnum()){ //this is a collapsed 3d element, icefront will be 2d
     125                        strcpy(icefront_type,"segment");
     126                }
     127                else if ((int)*(model->elements_type+2*element+0)==PattynEnum()){ //this is a real 3d element, icefront will be 3d.
     128                        strcpy(icefront_type,"quad");
     129                        i3=(int)*(model->segmentonneumann_diag+segment_width*i+2);
     130                        i4=(int)*(model->segmentonneumann_diag+segment_width*i+3);
     131                        icefront_node_ids[2]=i3;
     132                        icefront_node_ids[3]=i4;
     133                       
     134                        icefront_h[2]=model->thickness[i3-1];
     135                        icefront_h[3]=model->thickness[i4-1];
     136
     137                        icefront_b[2]=model->bed[i3-1];
     138                        icefront_b[3]=model->bed[i4-1];
     139                }
     140                else{
     141                        throw ErrorException(__FUNCT__," element type not supported yet");
     142                }
     143
     144                icefront = new Icefront(icefront_type,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b);
     145               
     146                loads->AddObject(icefront);
     147
     148        }
     149        /*Free data: */
     150        xfree((void**)&model->segmentonneumann_diag);
     151        xfree((void**)&model->elements_type);
     152        xfree((void**)&model->thickness);
     153        xfree((void**)&model->bed);
     154
     155               
     156        /*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */
     157
     158        /*First fetch penalties: */
     159        if (strcmp(model->meshtype,"3d")==0){
     160                ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
     161        }
     162
     163        counter=0;
     164        if(model->numpenalties){
     165
     166                /*First deal with internal grids: */
     167                for (i=0;i<model->numpenalties;i++){
     168                        if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
     169
     170                                for(j=0;j<model->numlayers-1;j++){
     171
     172                                        /*We are pairing grids along a vertical profile.*/
     173                                        grid1=(int)*(model->penalties+model->numlayers*i+j);
     174                                        grid2=(int)*(model->penalties+model->numlayers*i+j+1);
     175
     176                                        penpair_id=counter+1; //matlab indexing
     177                                        penpair_dof=1;
     178                                        penpair_node_ids[0]=grid1;
     179                                        penpair_node_ids[1]=grid2;
     180                                        penpair_penalty_offset=model->penalty_offset;
     181                                        penpair_numdofs=1;
     182
     183                                        penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
     184                                                                penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
     185                                        loads->AddObject(penpair);
     186
     187
     188                                        counter++;
     189
     190                                        penpair_id=counter+1; //matlab indexing
     191                                        penpair_dof=2;
     192                                        penpair_node_ids[0]=grid1;
     193                                        penpair_node_ids[1]=grid2;
     194                                        penpair_penalty_offset=model->penalty_offset;
     195                                        penpair_numdofs=1;
     196
     197                                        penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
     198                                                                penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
     199                                        loads->AddObject(penpair);
     200
     201
     202                                        counter++;
     203
     204                                } //for ( i=0; i<numpenalties; i++)
     205                        }
     206                }
     207        }
     208
     209        /*Free data: */
     210        xfree((void**)&model->penalties);
     211       
     212        /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
     213        /*First fetch rifts: */
     214        if(model->numrifts){
     215
     216                for(i=0;i<model->numrifts;i++){
     217                        riftpenaltypairs=model->riftspenaltypairs[i];
     218                        for(j=0;j<model->riftsnumpenaltypairs[i];j++){
     219
     220                                el1=(int)*(riftpenaltypairs+7*j+2);
     221                                #ifdef _PARALLEL_
     222                                if (model->epart[el1-1]!=my_rank){
     223                                        /*This load does not belong to this cluster node, as it references an element
     224                                         *that does not belong to this node's partition. Drop this 'j':*/
     225                                        continue;
     226                                }
     227                                #endif
     228
     229                                /*Ok, retrieve all the data needed to add a penalty between the two grids: */
     230                                el2=(int)*(riftpenaltypairs+7*j+3);
     231                                grid1=(int)*(riftpenaltypairs+7*j+0);
     232                                grid2=(int)*(riftpenaltypairs+7*j+1);
     233                                normal[0]=*(riftpenaltypairs+7*j+4);
     234                                normal[1]=*(riftpenaltypairs+7*j+5);
     235                                length=*(riftpenaltypairs+7*j+6);
     236                                friction=model->riftsfriction[i];
     237                                fill=model->riftsfill[i];
     238
     239                                penpair_id=counter+1; //matlab indexing
     240                                penpair_node_ids[0]=grid1;
     241                                penpair_node_ids[1]=grid2;
     242                                penpair_element_ids[0]=el1;
     243                                penpair_element_ids[1]=el2;
     244                                penpair_penalty_offset=model->penalty_offset;
     245                                penpair_penalty_lock=model->penalty_lock;
     246                                penpair_numdofs=2;
     247                                penpair_normal[0]=normal[0];
     248                                penpair_normal[1]=normal[1];
     249                                penpair_length=length;
     250                                penpair_friction=friction;
     251                                penpair_fill=fill;
     252
     253                                penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
     254                                                penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
     255                               
     256                                loads->AddObject(penpair);
     257
     258                                counter++;
     259                        }
     260                }
     261        }
     262
     263        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
     264         * datasets, it will not be redone: */
     265        loads->Presort();
     266
    26267        cleanup_and_return:
    27        
     268
     269        /*Free ressources:*/
     270        xfree((void**)&riftsnumpenaltypairs);
     271        for(i=0;i<model->numrifts;i++){
     272                double* temp=riftspenaltypairs[i];
     273                xfree((void**)&temp);
     274        }
     275        xfree((void**)&riftspenaltypairs);
     276        xfree((void**)&riftsfill);
     277        xfree((void**)&riftsfriction);
     278
     279
    28280        /*Assign output pointer: */
    29281        *ploads=loads;

  • issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp

    r306 r377  
    7979        int penta_thermal_steadystate;
    8080        double penta_viscosity_overshoot;
     81        double penta_reconditioning_number;
    8182
    8283        /*matpar constructor input: */
    … …  
    238239                                penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
    239240                                penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
    240                                 penta_thermal_steadystate,penta_viscosity_overshoot);
     241                                penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number);
    241242
    242243                /*Add penta element to elements dataset: */

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

    r358 r377  
    6161        model->gridonbed=NULL;
    6262        model->gridonsurface=NULL;
     63        model->gridonstokes=NULL;
     64        model->borderstokes=NULL;
    6365        model->thickness=NULL;
    6466        model->surface=NULL;
    … …  
    124126        model->yts=0;
    125127        model->viscosity_overshoot=0;
     128        model->reconditioning_number=0;
    126129        model->waitonlock=0;
    127130
    … …  
    195198        xfree((void**)&model->gridonbed);
    196199        xfree((void**)&model->gridonsurface);
     200        xfree((void**)&model->gridonstokes);
     201        xfree((void**)&model->borderstokes);
    197202        xfree((void**)&model->thickness);
    198203        xfree((void**)&model->surface);
    … …  
    328333        ModelFetchData((void**)&model->solverstring,NULL,NULL,model_handle,"solverstring","String",NULL);
    329334        ModelFetchData((void**)&model->viscosity_overshoot,NULL,NULL,model_handle,"viscosity_overshoot","Scalar",NULL);
     335        ModelFetchData((void**)&model->reconditioning_number,NULL,NULL,model_handle,"stokesreconditioning","Scalar",NULL);
    330336        ModelFetchData((void**)&model->waitonlock,NULL,NULL,model_handle,"waitonlock","Integer",NULL);
    331337

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

    r358 r377  
    5656        double* gridonbed;
    5757        double* gridonsurface;
     58        double* gridonstokes;
     59        double* borderstokes;
    5860        double* thickness;
    5961        double* surface;
    … …  
    119121        double  yts;
    120122        double  viscosity_overshoot;
     123        double  reconditioning_number;
    121124        int     waitonlock;
    122125

  • issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp

    r373 r377  
    9090        int penta_thermal_steadystate;
    9191        double penta_viscosity_overshoot;
     92        double penta_reconditioning_number;
    9293
    9394        /* node constructor input: */
    … …  
    286287                                        penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
    287288                                        penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
    288                                         penta_thermal_steadystate,penta_viscosity_overshoot);
     289                                        penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number);
    289290
    290291                        /*Add penta element to elements dataset: */
    … …  
    446447                if (strcmp(model->meshtype,"3d")==0){
    447448                        /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
    448                         if (model->gridonbed[i]==0){
     449                        if (model->deadgrids[i]){
    449450                                for(k=1;k<=node_numdofs;k++){
    450451                                        node->FreezeDof(k);

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

    r246 r377  
    1414        char* function_name;
    1515        mxArray* m;
     16        mxArray* inputs;
    1617        mxArray* p_g;
    1718        mxArray* u_g_obs;
    1819        mxArray* grad_g;
    1920        mxArray* n;
     21        mxArray* analysis_type;
    2022
    2123};

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

    r358 r377  
    2222                                double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel,
    2323                                int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6],
    24                                 int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot){
     24                                int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_reconditioning_number){
    2525       
    2626        int i;
    … …  
    5959        thermal_steadystate = penta_thermal_steadystate;
    6060        viscosity_overshoot = penta_viscosity_overshoot;
     61        reconditioning_number = penta_reconditioning_number;
    6162
    6263        return;
    … …  
    99100        printf("   thermal_steadystate: %i\n",thermal_steadystate);
    100101        printf("   viscosity_overshoot: %i\n",viscosity_overshoot);
     102        printf("   reconditioning_number: %i\n",reconditioning_number);
    101103        return;
    102104}
    … …  
    146148        memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
    147149        memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
     150        memcpy(marshalled_dataset,&reconditioning_number,sizeof(reconditioning_number));marshalled_dataset+=sizeof(reconditioning_number);
    148151       
    149152        *pmarshalled_dataset=marshalled_dataset;
    … …  
    182185                sizeof(thermal_steadystate) +
    183186                sizeof(viscosity_overshoot) +
     187                sizeof(reconditioning_number) +
    184188                sizeof(int); //sizeof(int) for enum type
    185189}
    … …  
    229233        memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
    230234        memcpy(&viscosity_overshoot,marshalled_dataset,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
     235        memcpy(&reconditioning_number,marshalled_dataset,sizeof(reconditioning_number));marshalled_dataset+=sizeof(reconditioning_number);
    231236
    232237        /*nodes and materials are not pointing to correct objects anymore:*/

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

    r358 r377  
    5353                int    thermal_steadystate;
    5454                double viscosity_overshoot;
     55                double reconditioning_number;
    5556       
    5657        public:
    … …  
    6061                                double p, double q, int    shelf, int    onbed, int    onsurface, double meanvel,double epsvel,
    6162                                int    collapse, double melting[6], double accumulation[6], double geothermalflux[6],
    62                                 int    artdiff, int    thermal_steadystate,double viscosity_overshoot);
     63                                int    artdiff, int    thermal_steadystate,double viscosity_overshoot,double reconditioning_number);
    6364                ~Penta();
    6465

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

    r8 r377  
    2323        double J;
    2424
    25         mxArray*   inputs[6];
     25        mxArray*   inputs[8];
    2626        mxArray*   psearch_scalar=NULL;
    2727        mxArray*   mxJ=NULL;
    … …  
    3232       
    3333        inputs[1]=optargs->m;
    34         inputs[2]=optargs->p_g;
    35         inputs[3]=optargs->u_g_obs;
    36         inputs[4]=optargs->grad_g;
    37         inputs[5]=optargs->n;
     34        inputs[2]=optargs->inputs;
     35        inputs[3]=optargs->p_g;
     36        inputs[4]=optargs->u_g_obs;
     37        inputs[5]=optargs->grad_g;
     38        inputs[6]=optargs->n;
     39        inputs[7]=optargs->analysis_type;
    3840
    39         mexCallMATLAB( 1, &mxJ, 6, (mxArray**)inputs, optargs->function_name);
     41        mexCallMATLAB( 1, &mxJ, 8, (mxArray**)inputs, optargs->function_name);
    4042
    4143        /*extract misfit from mxArray*/

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

    r360 r377  
    5959        md.gridonpattyn=NaN;
    6060        md.gridonstokes=NaN;
     61        md.borderstokes=NaN;
    6162
    6263        %Stokes mesh

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

    r368 r377  
    5050WriteData(fid,md.gridonbed,'Mat','gridonbed');
    5151WriteData(fid,md.gridonsurface,'Mat','gridonsurface');
     52WriteData(fid,md.gridonstokes,'Mat','gridonstokes');
     53WriteData(fid,md.borderstokes,'Mat','borderstokes');
    5254
    5355
    … …  
    126128WriteData(fid,md.solverstring,'String','solverstring');
    127129WriteData(fid,md.viscosity_overshoot,'Scalar','viscosity_overshoot');
     130WriteData(fid,md.stokesreconditioning,'Scalar','reconditioning_number');
    128131WriteData(fid,md.waitonlock,'Integer','waitonlock');
    129132

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

    r304 r377  
    120120md.deadgrids=deadgrids;
    121121
     122%figure out the border stokes grids
     123stokes_elements=find(md.elements_type(:,2)==stokesenum()); %find the elements on the stokes domain
     124borderflags=zeros(md.numberofgrids,1);
     125borderflags(md.elements(stokes_elements,:))=1; %find all the grids of the elements on stokes domain, ie stokes grids and borderstokes
     126md.borderstokes=borderflags-md.gridonstokes; %remove stokes grids from this list
     127
    122128%figure out solution types
    123129md.ishutter=double(any(md.elements_type(:,1)==hutterenum()));

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

    r370 r377  
    1 function grad_g=GradJCompute(m,inputs, u_g_obs);
     1function grad_g=GradJCompute(m,inputs, u_g_obs,analysis_type);
    22
    33%Recover solution for this stiffness and right hand side:
    … …  
    55        disp('         computing velocities...')
    66end
    7 [u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,'diagnostic_horiz');
     7[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type);
    88
    99%Buid Du, difference between observed velocity and model velocity.

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

    r276 r377  
    1919        options=ControlOptions(m.parameters);
    2020
     21        %initialize inputs, ie m.nparameters on which we invert.
     22        inputs=inputlist;
     23        inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
     24
    2125        for n=1:m.parameters.nsteps,
    2226               
    … …  
    2630                disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m.parameters.nsteps)]));
    2731
    28                 %initialize inputs, ie m.nparameters on which we invert.
    29                 inputs=inputlist;
     32                %update inputs with new parameter and fit
    3033                inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
    31                 inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
    3234                inputs=add(inputs,'fit',m.parameters.fit(n),'double');
    3335
    … …  
    3638
    3739                disp('      computing gradJ...');
    38                 c(n).grad_g=GradJCompute(m,inputs,u_g_obs);
     40                c(n).grad_g=GradJCompute(m,inputs,u_g_obs,md.analysis_type);
    3941                disp('      done.');
    4042
    … …  
    5355               
    5456                disp('      optimizing along gradient direction...');
    55                 [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,p_g,u_g_obs,c(n).grad_g,n);
    56                 %[search_scalar c(n).J]=fminbnd('objectivefunctionC',0,1,options,m,p_g,u_g_obs,c(n).grad_g,n);
     57                md.analysis_type
     58                [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m,inputs,p_g,u_g_obs,c(n).grad_g,n,md.analysis_type);
    5759                disp('      done.');
    5860
    … …  
    7274                %some temporary saving
    7375                if(mod(n,5)==0),
    74                         solution=controlfinalsol(c,m,p_g);
     76                        solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
    7577                        save temporary_control_results solution
    7678                end
    … …  
    8183        %Create final solution
    8284        disp('      preparing final velocity solution...');
    83         solution=controlfinalsol(c,m,p_g);
     85        solution=controlfinalsol(c,m,p_g,inputs,md.analysis_type);
    8486        disp('      done.');
    8587       

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

    r372 r377  
    1 function solution=controlfinalsol(c,m,p_g);
     1function solution=controlfinalsol(c,m,p_g,inputs,analysis_type);
    22
    33%From parameters, build inputs for icediagnostic_core, using the final parameters
    4 inputs=inputlist;
    5 inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
    64inputs=add(inputs,m.parameters.control_type,p_g,'doublevec',2,m.parameters.numberofnodes);
    7 u_g=diagnostic_core_nonlinear(m,inputs,'diagnostic_horiz');
     5u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
    86
    97%Build partitioning vectors to recover solution

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

    r358 r377  
    7272
    7373                %"recondition" pressure
    74                 p_g=p_g/m_ds.parameters.stokesreconditioning;
     74                p_g=p_g/m_ds.parameters.reconditioning_number;
    7575
    7676                displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
    … …  
    9595       
    9696                %"decondition" pressure
    97                 p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
     97                p_g=u_g(4:4:end)*m_dh.parameters.reconditioning_number;
    9898        end
    9999end

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

    r370 r377  
    1 function J =objectivefunctionC(search_scalar,m,p_g,u_g_obs,grad_g,n);
     1function J =objectivefunctionC(search_scalar,m,inputs,p_g,u_g_obs,grad_g,n,analysis_type);
    22       
    33%recover some parameters
    … …  
    1010
    1111%Plug parameter into inputs
    12 inputs=inputlist;
    13 inputs=add(inputs,'velocity',m.parameters.u_g,'doublevec',3,m.parameters.numberofnodes);
    1412inputs=add(inputs,m.parameters.control_type,parameter,'doublevec',2,m.parameters.numberofnodes);
    1513
    1614%Run diagnostic with updated parameters.
    17 u_g=diagnostic_core_nonlinear(m,inputs,'diagnostic_horiz');
     15u_g=diagnostic_core_nonlinear(m,inputs,analysis_type);
    1816
    1917%Compute misfit for this velocity field.
    20 inputs=add(inputs,'fit',m.parameters.fit(n),'double');
    2118J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters, u_g, u_g_obs,inputs);

  • issm/trunk/src/mex/ControlOptimization/ControlOptimization.cpp

    r8 r377  
    4545        optargs.function_name=function_name;
    4646        optargs.m=MODEL;
    47         optargs.p_g=PG;
     47        optargs.p_g=PG;
     48        optargs.inputs=INPUTS;
    4849        optargs.u_g_obs=VELOBS;
    4950        optargs.grad_g=GRADIENT;
    5051        optargs.n=STEP;
     52        optargs.analysis_type=ANALYSIS;
    5153
    5254        optpars.xmin=xmin;
    … …  
    7173{
    7274        _printf_("\n");
    73         _printf_("   usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,p_g,u_g_obs,grad_g,step)\n");
     75        _printf_("   usage: [search_scalar J] = %s(function_name,xmin,xmax,options,m,inputs,p_g,u_g_obs,grad_g,step,analysis_type)\n",__FUNCT__);
    7476        _printf_("\n");
    7577}

  • issm/trunk/src/mex/ControlOptimization/ControlOptimization.h

    r8 r377  
    2222#define OPTIONS (mxArray*)prhs[3]
    2323#define MODEL (mxArray*)prhs[4]
    24 #define PG (mxArray*)prhs[5]
    25 #define VELOBS (mxArray*)prhs[6]
    26 #define GRADIENT (mxArray*)prhs[7]
    27 #define STEP (mxArray*)prhs[8]
     24#define INPUTS (mxArray*)prhs[5]
     25#define PG (mxArray*)prhs[6]
     26#define VELOBS (mxArray*)prhs[7]
     27#define GRADIENT (mxArray*)prhs[8]
     28#define STEP (mxArray*)prhs[9]
     29#define ANALYSIS (mxArray*)prhs[10]
    2830
    2931/* serial output macros: */
    … …  
    3537#define NLHS  2
    3638#undef NRHS
    37 #define NRHS  9
     39#define NRHS  11
    3840
    3941
Note: See TracChangeset for help on using the changeset viewer.