Changeset 3201


Ignore:
Timestamp:
03/05/10 13:47:09 (15 years ago)
Author:
Eric.Larour
Message:

Added freezing of constraints for rifts that zigzag.
Also added max_nonlinear_iterations, to control how many iterations can be carried out in Newton iterations.
Add min_mechanical_constraints, acts like min_thermal_constraints, but for rift constraints.

Location:
issm/trunk/src
Files:
15 edited

Legend:

Unmodified
Added
Removed

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

    r2383 r3201  
    109109        parameters->AddObject(param);
    110110
     111        /*max_nonlinear_iterations: */
     112        count++;
     113        param= new Param(count,"max_nonlinear_iterations",DOUBLE);
     114        param->SetDouble(iomodel->max_nonlinear_iterations);
     115        parameters->AddObject(param);
     116
    111117        /*yts: */
    112118        count++;
    … …  
    191197        param= new Param(count,"min_thermal_constraints",DOUBLE);
    192198        param->SetDouble(iomodel->min_thermal_constraints);
     199        parameters->AddObject(param);
     200
     201        /*min_mechanical_constraints: */
     202        count++;
     203        param= new Param(count,"min_mechanical_constraints",DOUBLE);
     204        param->SetDouble(iomodel->min_mechanical_constraints);
    193205        parameters->AddObject(param);
    194206

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

    r3105 r3201  
    6161        int riftfront_penalty_lock;
    6262        bool riftfront_active;
     63        bool riftfront_frozen;
    6364        int  riftfront_counter;
    6465        bool riftfront_prestable;
    … …  
    227228
    228229                        riftfront_active=0;
     230                        riftfront_frozen=0;
    229231                        riftfront_counter=0;
    230232                        riftfront_prestable=0;
    231233                       
    232                         riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
     234                        riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_frozen,riftfront_counter,riftfront_prestable,riftfront_shelf);
    233235
    234236                        loads->AddObject(riftfront);

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

    r2959 r3201  
    129129        iomodel->eps_rel=0;
    130130        iomodel->eps_abs=0;
     131        iomodel->max_nonlinear_iterations=0;
    131132        iomodel->dt=0;
    132133        iomodel->ndt=0;
    … …  
    151152        iomodel->thermalconductivity=0;
    152153        iomodel->min_thermal_constraints=0;
     154        iomodel->min_mechanical_constraints=0;
    153155        iomodel->stabilize_constraints=0;
    154156        iomodel->mixed_layer_capacity=0;
    … …  
    355357        IoModelFetchData(&iomodel->eps_rel,iomodel_handle,"eps_rel");
    356358        IoModelFetchData(&iomodel->eps_abs,iomodel_handle,"eps_abs");
     359        IoModelFetchData(&iomodel->max_nonlinear_iterations,iomodel_handle,"max_nonlinear_iterations");
    357360        IoModelFetchData(&iomodel->dt,iomodel_handle,"dt");
    358361        IoModelFetchData(&iomodel->ndt,iomodel_handle,"ndt");
    … …  
    376379        IoModelFetchData(&iomodel->thermalconductivity,iomodel_handle,"thermalconductivity");
    377380        IoModelFetchData(&iomodel->min_thermal_constraints,iomodel_handle,"min_thermal_constraints");
     381        IoModelFetchData(&iomodel->min_mechanical_constraints,iomodel_handle,"min_mechanical_constraints");
    378382        IoModelFetchData(&iomodel->stabilize_constraints,iomodel_handle,"stabilize_constraints");
    379383        IoModelFetchData(&iomodel->mixed_layer_capacity,iomodel_handle,"mixed_layer_capacity");

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

    r2959 r3201  
    135135        double  eps_rel;
    136136        double  eps_abs;
     137        double  max_nonlinear_iterations;
    137138        double  dt,ndt;
    138139        double  penalty_offset;
    … …  
    152153        double  heatcapacity,thermalconductivity;
    153154        int    min_thermal_constraints;
     155        int    min_mechanical_constraints;
    154156        int    stabilize_constraints;
    155157        double mixed_layer_capacity;

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

    r2333 r3201  
    2525        int converged=0;
    2626        int num_unstable_constraints=0;
     27        int min_mechanical_constraints=0;
     28
     29        /*recover parameters: */
     30        parameters->FindParam(&min_mechanical_constraints,"min_mechanical_constraints");
    2731
    2832        /*First, get loads configured: */
    … …  
    3236         * management routine, otherwise, skip : */
    3337        if (RiftIsPresent(loads)){
    34                 RiftConstraints(&converged,&num_unstable_constraints,loads,inputs,analysis_type,sub_analysis_type);
     38                RiftConstraints(&converged,&num_unstable_constraints,loads,inputs,min_mechanical_constraints,analysis_type,sub_analysis_type);
    3539        }
    3640        else if(loads->MeltingIsPresent()){

  • issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp

    r1904 r3201  
    55#include "./RiftConstraints.h"
    66#include "../EnumDefinitions/EnumDefinitions.h"
     7#include "../include/macros.h"
     8#include "../shared/shared.h"
    79
    810#define _ZIGZAGCOUNTER_
    … …  
    1012#undef __FUNCT__
    1113#define __FUNCT__ "RiftConstrain"
    12 int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
     14int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int min_mechanical_constraints,int analysis_type,int sub_analysis_type){
    1315
    1416        int num_unstable_constraints=0;
    1517        int converged=0;
    1618        int potential;
     19        extern int my_rank;
    1720
    1821        Constrain(&num_unstable_constraints,loads,inputs,analysis_type);
    1922        if(num_unstable_constraints==0)converged=1;
     23       
     24       
     25        if(IsFrozen(loads)){
     26                converged=1;
     27                num_unstable_constraints=0;
     28        }
     29        else if(num_unstable_constraints<=min_mechanical_constraints){
     30                _printf_("   freezing constraints\n");
     31                FreezeConstraints(loads,inputs,analysis_type);
     32        }
    2033
    2134        /*Assign output pointers: */
    … …  
    227240
    228241#undef __FUNCT__
     242#define __FUNCT__ "FreezeConstraints"
     243
     244void FreezeConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type){
     245
     246        int                     i;
     247       
     248        /* generic object pointer: */
     249        Riftfront* riftfront=NULL;
     250
     251        /*Enforce constraints: */
     252        for (i=0;i<loads->Size();i++){
     253
     254                if (RiftfrontEnum()==loads->GetEnum(i)){
     255
     256                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     257
     258                        riftfront->FreezeConstraints(inputs,analysis_type);
     259
     260                }
     261        }
     262
     263}
     264
     265#undef __FUNCT__
     266#define __FUNCT__ "IsFrozen"
     267int IsFrozen(DataSet* loads){
     268
     269        int                     i;
     270       
     271        /* generic object pointer: */
     272        Riftfront* riftfront=NULL;
     273        int found=0;
     274        int mpi_found=0;
     275
     276        /*Enforce constraints: */
     277        for (i=0;i<loads->Size();i++){
     278
     279                if (RiftfrontEnum()==loads->GetEnum(i)){
     280
     281                        riftfront=(Riftfront*)loads->GetObjectByOffset(i);
     282                        if (riftfront->IsFrozen()){
     283                                found=1;
     284                                break;
     285                        }
     286                }
     287        }
     288       
     289        /*Is there just one found? that would mean we have frozen! : */
     290        #ifdef _PARALLEL_
     291        MPI_Reduce (&found,&mpi_found,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD );
     292        MPI_Bcast(&mpi_found,1,MPI_DOUBLE,0,MPI_COMM_WORLD);               
     293        found=mpi_found;
     294        #endif
     295
     296        return found;
     297}
     298
     299#undef __FUNCT__
    229300#define __FUNCT__ "MaxPenetrationInInputs"
    230301

  • issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h

    r1784 r3201  
    1010#include "../DataSet/DataSet.h"
    1111
    12 int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     12int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int min_mechanical_constraints,int analysis_type,int sub_analysis_type);
    1313
    1414int RiftIsPresent(DataSet* loads);
    … …  
    2222int Constrain(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type_enum);
    2323
     24void FreezeConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type);
     25
    2426int MaxPenetrationInInputs(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum);
    2527
    … …  
    2729
    2830int IsMaterialStable(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum);
     31
     32int IsFrozen(DataSet* loads);
    2933#endif

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

    r2908 r3201  
    2828/*}}}1*/
    2929/*FUNCTION Riftfront::creation {{{1*/
    30 Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
     30Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,bool riftfront_frozen, int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
    3131
    3232        int i;
    … …  
    6161        penalty_lock=riftfront_penalty_lock;
    6262        active=riftfront_active;
     63        frozen=riftfront_frozen;
    6364        counter=riftfront_counter;
    6465        prestable=riftfront_prestable;
    … …  
    114115        memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
    115116        memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
     117        memcpy(marshalled_dataset,&frozen,sizeof(frozen));marshalled_dataset+=sizeof(frozen);
    116118        memcpy(marshalled_dataset,&counter,sizeof(counter));marshalled_dataset+=sizeof(counter);
    117119        memcpy(marshalled_dataset,&prestable,sizeof(prestable));marshalled_dataset+=sizeof(prestable);
    … …  
    144146                sizeof(penalty_lock)+
    145147                sizeof(active)+
     148                sizeof(frozen)+
    146149                sizeof(counter)+
    147150                sizeof(prestable)+
    … …  
    182185        memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
    183186        memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
     187        memcpy(&frozen,marshalled_dataset,sizeof(frozen));marshalled_dataset+=sizeof(frozen);
    184188        memcpy(&counter,marshalled_dataset,sizeof(counter));marshalled_dataset+=sizeof(counter);
    185189        memcpy(&prestable,marshalled_dataset,sizeof(prestable));marshalled_dataset+=sizeof(prestable);
    … …  
    237241        inputs=(ParameterInputs*)vinputs;
    238242
     243        /*Is this constraint frozen? In which case we don't touch: */
     244        if (this->frozen){
     245                *punstable=0;
     246                return 1;
     247        }
    239248
    240249        /*First recover parameter inputs: */
    … …  
    322331        printf("penalty_lock %i\n",penalty_lock);
    323332        printf("active %i\n",active);
     333        printf("frozen %i\n",frozen);
    324334        printf("counter %i\n",counter);
    325335        printf("prestable %i\n",prestable);
    … …  
    348358        printf("penalty_lock %i\n",penalty_lock);
    349359        printf("active %i\n",active);
     360        printf("frozen %i\n",frozen);
    350361        printf("counter %i\n",counter);
    351362        printf("prestable %i\n",prestable);
    … …  
    709720}
    710721/*}}}1*/
    711 /*FUNCTION Riftfront::PotentialUnstableConstraint {{{1*/
     722/*FUNCTION Riftfront::PotentialUnstableC
     723 * nstraint {{{1*/
    712724#undef __FUNCT__
    713725#define __FUNCT__ "Riftfront::PotentialUnstableConstraint"
    … …  
    831843}
    832844/*}}}1*/
     845/*FUNCTION Riftfront::FreezeConstraints{{{1*/
     846#undef __FUNCT__
     847#define __FUNCT__ "Riftfront::FreezeConstraints"
     848void   Riftfront::FreezeConstraints(void* vinputs, int analysis_type){
     849
     850        /*Just set frozen flag to 1: */
     851        this->frozen=1;
     852
     853}
     854/*}}}1*/
     855/*FUNCTION Riftfront::IsFrozen{{{1*/
     856#undef __FUNCT__
     857#define __FUNCT__ "Riftfront::IsFrozen"
     858bool   Riftfront::IsFrozen(void){
     859
     860        /*Just set frozen flag to 1: */
     861        if(this->frozen)return 1;
     862        else return 0;
     863}
     864/*}}}1*/

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

    r1805 r3201  
    5050                /*computational: */
    5151                bool        active;
     52                bool        frozen;
    5253                int         counter;
    5354                bool        prestable;
    … …  
    5859
    5960                Riftfront();
    60                 Riftfront(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double fraction, double fractionincrement, double penalty_offset, int penalty_lock,bool active,int counter,bool prestable,bool shelf);
     61                Riftfront(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double fraction, double fractionincrement, double penalty_offset, int penalty_lock,bool active,bool frozen, int counter,bool prestable,bool shelf);
    6162                ~Riftfront();
    6263
    … …  
    8182                int   PreConstrain(int* punstable, void* inputs, int analysis_type);
    8283                int   Constrain(int* punstable, void* inputs, int analysis_type);
     84                void  FreezeConstraints(void* inputs, int analysis_type);
     85                bool  IsFrozen(void);
    8386                int   Penetration(double* ppenetration, void* inputs, int analysis_type);
    8487                int   MaxPenetration(double* ppenetration, void* inputs, int analysis_type);

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

    r3189 r3201  
    47724772
    47734773        /*Update internal data if inputs holds new values: */
    4774         if (id==1) printf("WARNING if QMU: no hydrostatic equilibrium is applied here (conflict with prognostic, which does not have matpar)\n");
     4774        //if (id==1) printf("WARNING if QMU: no hydrostatic equilibrium is applied here (conflict with prognostic, which does not have matpar)\n");
    47754775        //For now
    47764776        inputs->Recover("thickness",&h[0],1,dofs,3,(void**)nodes);

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

    r3186 r3201  
    3333        int count;
    3434        int numberofnodes;
     35        int min_mechanical_constraints;
     36        int max_nonlinear_iterations;
    3537
    3638        /*parameters:*/
    … …  
    4749        fem->FindParam(&solver_string,"solverstring");
    4850        fem->FindParam(&verbose,"verbose");
    49 
     51        fem->FindParam(&min_mechanical_constraints,"min_mechanical_constraints");
     52        fem->FindParam(&max_nonlinear_iterations,"max_nonlinear_iterations");
     53       
    5054        /*Were loads requested as output? : */
    5155        if(!input_loads){
    … …  
    129133
    130134                //rift convergence
    131                 if (!constraints_converged) converged=0;
     135                if (!constraints_converged) {
     136                        if (converged){
     137                                if (num_unstable_constraints <= min_mechanical_constraints) converged=1;
     138                                else converged=0;
     139                        }
     140                }
    132141
    133142                /*Increase count: */
    134143                count++;
    135144                if(converged==1)break;
    136                 if(count>=1000){
    137                         _printf_("   maximum number of iterations exceeded\n");
     145                if(count>=max_nonlinear_iterations){
     146                        _printf_("   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
    138147                        break;
    139148                }
    … …  
    155164
    156165        }
    157        
     166
    158167        /*Delete loads only if no ouput was requested: */
    159168        if(!input_loads)delete loads;

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

    r3136 r3201  
    116116        md.thermal_exchange_velocity=0;
    117117        md.min_thermal_constraints=0;
     118        md.min_mechanical_constraints=0;
    118119        md.stabilize_constraints=0;
    119120       
    … …  
    165166        md.eps_rel=0;
    166167        md.eps_abs=0;
     168        md.max_nonlinear_iterations=0;
    167169        md.sparsity=0;
    168170        md.connectivity=0;

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

    r2892 r3201  
    6565md.eps_abs=10;
    6666
     67%maximum of non-linear iterations.
     68md.max_nonlinear_iterations=100;
     69
    6770%sparsity
    6871md.sparsity=0.001;
    … …  
    114117%parameter is often used.
    115118md.min_thermal_constraints=0;
     119md.min_mechanical_constraints=0;
    116120
    117121%Transient parameters

  • issm/trunk/src/m/classes/public/display/displaydiagnostic.m

    r2330 r3201  
    1515fielddisplay(md,'eps_rel','velocity relative convergence criterion, NaN -> not applied');
    1616fielddisplay(md,'eps_abs','velocity absolute convergence criterion, NaN -> not applied');
     17fielddisplay(md,'max_nonlinear_iterations','maximum number of nonlinear iterations');
    1718fielddisplay(md,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)');
    1819
    … …  
    2324disp(sprintf('\n      %s','Penalties:'));
    2425fielddisplay(md,'penalty_offset','offset used by penalties: penalty = Kmax*10^offset');
     26fielddisplay(md,'min_mechanical_constraints','threshold for instability of mechanical constraints');
    2527
    2628disp(sprintf('\n      %s','Memory management:'));

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

    r2959 r3201  
    122122WriteData(fid,md.eps_rel,'Scalar','eps_rel');
    123123WriteData(fid,md.eps_abs,'Scalar','eps_abs');
     124WriteData(fid,md.max_nonlinear_iterations,'Scalar','max_nonlinear_iterations');
    124125WriteData(fid,md.dt,'Scalar','dt');
    125126WriteData(fid,md.ndt,'Scalar','ndt');
    … …  
    143144WriteData(fid,md.thermalconductivity,'Scalar','thermalconductivity');
    144145WriteData(fid,md.min_thermal_constraints,'Integer','min_thermal_constraints');
     146WriteData(fid,md.min_mechanical_constraints,'Integer','min_mechanical_constraints');
    145147WriteData(fid,md.mixed_layer_capacity,'Scalar','mixed_layer_capacity');
    146148WriteData(fid,md.thermal_exchange_velocity,'Scalar','thermal_exchange_velocity');
Note: See TracChangeset for help on using the changeset viewer.