Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3200) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3201) @@ -109,4 +109,10 @@ parameters->AddObject(param); + /*max_nonlinear_iterations: */ + count++; + param= new Param(count,"max_nonlinear_iterations",DOUBLE); + param->SetDouble(iomodel->max_nonlinear_iterations); + parameters->AddObject(param); + /*yts: */ count++; @@ -191,4 +197,10 @@ param= new Param(count,"min_thermal_constraints",DOUBLE); param->SetDouble(iomodel->min_thermal_constraints); + parameters->AddObject(param); + + /*min_mechanical_constraints: */ + count++; + param= new Param(count,"min_mechanical_constraints",DOUBLE); + param->SetDouble(iomodel->min_mechanical_constraints); parameters->AddObject(param); Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 3200) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 3201) @@ -61,4 +61,5 @@ int riftfront_penalty_lock; bool riftfront_active; + bool riftfront_frozen; int riftfront_counter; bool riftfront_prestable; @@ -227,8 +228,9 @@ riftfront_active=0; + riftfront_frozen=0; riftfront_counter=0; riftfront_prestable=0; - 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); + 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); loads->AddObject(riftfront); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3200) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3201) @@ -129,4 +129,5 @@ iomodel->eps_rel=0; iomodel->eps_abs=0; + iomodel->max_nonlinear_iterations=0; iomodel->dt=0; iomodel->ndt=0; @@ -151,4 +152,5 @@ iomodel->thermalconductivity=0; iomodel->min_thermal_constraints=0; + iomodel->min_mechanical_constraints=0; iomodel->stabilize_constraints=0; iomodel->mixed_layer_capacity=0; @@ -355,4 +357,5 @@ IoModelFetchData(&iomodel->eps_rel,iomodel_handle,"eps_rel"); IoModelFetchData(&iomodel->eps_abs,iomodel_handle,"eps_abs"); + IoModelFetchData(&iomodel->max_nonlinear_iterations,iomodel_handle,"max_nonlinear_iterations"); IoModelFetchData(&iomodel->dt,iomodel_handle,"dt"); IoModelFetchData(&iomodel->ndt,iomodel_handle,"ndt"); @@ -376,4 +379,5 @@ IoModelFetchData(&iomodel->thermalconductivity,iomodel_handle,"thermalconductivity"); IoModelFetchData(&iomodel->min_thermal_constraints,iomodel_handle,"min_thermal_constraints"); + IoModelFetchData(&iomodel->min_mechanical_constraints,iomodel_handle,"min_mechanical_constraints"); IoModelFetchData(&iomodel->stabilize_constraints,iomodel_handle,"stabilize_constraints"); IoModelFetchData(&iomodel->mixed_layer_capacity,iomodel_handle,"mixed_layer_capacity"); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3200) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3201) @@ -135,4 +135,5 @@ double eps_rel; double eps_abs; + double max_nonlinear_iterations; double dt,ndt; double penalty_offset; @@ -152,4 +153,5 @@ double heatcapacity,thermalconductivity; int min_thermal_constraints; + int min_mechanical_constraints; int stabilize_constraints; double mixed_layer_capacity; Index: /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 3200) +++ /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.cpp (revision 3201) @@ -25,4 +25,8 @@ int converged=0; int num_unstable_constraints=0; + int min_mechanical_constraints=0; + + /*recover parameters: */ + parameters->FindParam(&min_mechanical_constraints,"min_mechanical_constraints"); /*First, get loads configured: */ @@ -32,5 +36,5 @@ * management routine, otherwise, skip : */ if (RiftIsPresent(loads)){ - RiftConstraints(&converged,&num_unstable_constraints,loads,inputs,analysis_type,sub_analysis_type); + RiftConstraints(&converged,&num_unstable_constraints,loads,inputs,min_mechanical_constraints,analysis_type,sub_analysis_type); } else if(loads->MeltingIsPresent()){ Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 3200) +++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 3201) @@ -5,4 +5,6 @@ #include "./RiftConstraints.h" #include "../EnumDefinitions/EnumDefinitions.h" +#include "../include/macros.h" +#include "../shared/shared.h" #define _ZIGZAGCOUNTER_ @@ -10,12 +12,23 @@ #undef __FUNCT__ #define __FUNCT__ "RiftConstrain" -int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ +int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int min_mechanical_constraints,int analysis_type,int sub_analysis_type){ int num_unstable_constraints=0; int converged=0; int potential; + extern int my_rank; Constrain(&num_unstable_constraints,loads,inputs,analysis_type); if(num_unstable_constraints==0)converged=1; + + + if(IsFrozen(loads)){ + converged=1; + num_unstable_constraints=0; + } + else if(num_unstable_constraints<=min_mechanical_constraints){ + _printf_(" freezing constraints\n"); + FreezeConstraints(loads,inputs,analysis_type); + } /*Assign output pointers: */ @@ -227,4 +240,62 @@ #undef __FUNCT__ +#define __FUNCT__ "FreezeConstraints" + +void FreezeConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + + riftfront->FreezeConstraints(inputs,analysis_type); + + } + } + +} + +#undef __FUNCT__ +#define __FUNCT__ "IsFrozen" +int IsFrozen(DataSet* loads){ + + int i; + + /* generic object pointer: */ + Riftfront* riftfront=NULL; + int found=0; + int mpi_found=0; + + /*Enforce constraints: */ + for (i=0;iSize();i++){ + + if (RiftfrontEnum()==loads->GetEnum(i)){ + + riftfront=(Riftfront*)loads->GetObjectByOffset(i); + if (riftfront->IsFrozen()){ + found=1; + break; + } + } + } + + /*Is there just one found? that would mean we have frozen! : */ + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_DOUBLE,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; +} + +#undef __FUNCT__ #define __FUNCT__ "MaxPenetrationInInputs" Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h (revision 3200) +++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.h (revision 3201) @@ -10,5 +10,5 @@ #include "../DataSet/DataSet.h" -int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +int RiftConstraints(int* pconverged, int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int min_mechanical_constraints,int analysis_type,int sub_analysis_type); int RiftIsPresent(DataSet* loads); @@ -22,4 +22,6 @@ int Constrain(int* pnum_unstable_constraints,DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); +void FreezeConstraints(DataSet* loads,ParameterInputs* inputs,int analysis_type); + int MaxPenetrationInInputs(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); @@ -27,3 +29,5 @@ int IsMaterialStable(DataSet* loads,ParameterInputs* inputs,int analysis_type_enum); + +int IsFrozen(DataSet* loads); #endif Index: /issm/trunk/src/c/objects/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Riftfront.cpp (revision 3200) +++ /issm/trunk/src/c/objects/Riftfront.cpp (revision 3201) @@ -28,5 +28,5 @@ /*}}}1*/ /*FUNCTION Riftfront::creation {{{1*/ -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){ +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,bool riftfront_frozen, int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){ int i; @@ -61,4 +61,5 @@ penalty_lock=riftfront_penalty_lock; active=riftfront_active; + frozen=riftfront_frozen; counter=riftfront_counter; prestable=riftfront_prestable; @@ -114,4 +115,5 @@ memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active); + memcpy(marshalled_dataset,&frozen,sizeof(frozen));marshalled_dataset+=sizeof(frozen); memcpy(marshalled_dataset,&counter,sizeof(counter));marshalled_dataset+=sizeof(counter); memcpy(marshalled_dataset,&prestable,sizeof(prestable));marshalled_dataset+=sizeof(prestable); @@ -144,4 +146,5 @@ sizeof(penalty_lock)+ sizeof(active)+ + sizeof(frozen)+ sizeof(counter)+ sizeof(prestable)+ @@ -182,4 +185,5 @@ memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active); + memcpy(&frozen,marshalled_dataset,sizeof(frozen));marshalled_dataset+=sizeof(frozen); memcpy(&counter,marshalled_dataset,sizeof(counter));marshalled_dataset+=sizeof(counter); memcpy(&prestable,marshalled_dataset,sizeof(prestable));marshalled_dataset+=sizeof(prestable); @@ -237,4 +241,9 @@ inputs=(ParameterInputs*)vinputs; + /*Is this constraint frozen? In which case we don't touch: */ + if (this->frozen){ + *punstable=0; + return 1; + } /*First recover parameter inputs: */ @@ -322,4 +331,5 @@ printf("penalty_lock %i\n",penalty_lock); printf("active %i\n",active); + printf("frozen %i\n",frozen); printf("counter %i\n",counter); printf("prestable %i\n",prestable); @@ -348,4 +358,5 @@ printf("penalty_lock %i\n",penalty_lock); printf("active %i\n",active); + printf("frozen %i\n",frozen); printf("counter %i\n",counter); printf("prestable %i\n",prestable); @@ -709,5 +720,6 @@ } /*}}}1*/ -/*FUNCTION Riftfront::PotentialUnstableConstraint {{{1*/ +/*FUNCTION Riftfront::PotentialUnstableC + * nstraint {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Riftfront::PotentialUnstableConstraint" @@ -831,2 +843,22 @@ } /*}}}1*/ +/*FUNCTION Riftfront::FreezeConstraints{{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::FreezeConstraints" +void Riftfront::FreezeConstraints(void* vinputs, int analysis_type){ + + /*Just set frozen flag to 1: */ + this->frozen=1; + +} +/*}}}1*/ +/*FUNCTION Riftfront::IsFrozen{{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::IsFrozen" +bool Riftfront::IsFrozen(void){ + + /*Just set frozen flag to 1: */ + if(this->frozen)return 1; + else return 0; +} +/*}}}1*/ Index: /issm/trunk/src/c/objects/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Riftfront.h (revision 3200) +++ /issm/trunk/src/c/objects/Riftfront.h (revision 3201) @@ -50,4 +50,5 @@ /*computational: */ bool active; + bool frozen; int counter; bool prestable; @@ -58,5 +59,5 @@ Riftfront(); - 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); + 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); ~Riftfront(); @@ -81,4 +82,6 @@ int PreConstrain(int* punstable, void* inputs, int analysis_type); int Constrain(int* punstable, void* inputs, int analysis_type); + void FreezeConstraints(void* inputs, int analysis_type); + bool IsFrozen(void); int Penetration(double* ppenetration, void* inputs, int analysis_type); int MaxPenetration(double* ppenetration, void* inputs, int analysis_type); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3200) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3201) @@ -4772,5 +4772,5 @@ /*Update internal data if inputs holds new values: */ - if (id==1) printf("WARNING if QMU: no hydrostatic equilibrium is applied here (conflict with prognostic, which does not have matpar)\n"); + //if (id==1) printf("WARNING if QMU: no hydrostatic equilibrium is applied here (conflict with prognostic, which does not have matpar)\n"); //For now inputs->Recover("thickness",&h[0],1,dofs,3,(void**)nodes); Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 3200) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 3201) @@ -33,4 +33,6 @@ int count; int numberofnodes; + int min_mechanical_constraints; + int max_nonlinear_iterations; /*parameters:*/ @@ -47,5 +49,7 @@ fem->FindParam(&solver_string,"solverstring"); fem->FindParam(&verbose,"verbose"); - + fem->FindParam(&min_mechanical_constraints,"min_mechanical_constraints"); + fem->FindParam(&max_nonlinear_iterations,"max_nonlinear_iterations"); + /*Were loads requested as output? : */ if(!input_loads){ @@ -129,11 +133,16 @@ //rift convergence - if (!constraints_converged) converged=0; + if (!constraints_converged) { + if (converged){ + if (num_unstable_constraints <= min_mechanical_constraints) converged=1; + else converged=0; + } + } /*Increase count: */ count++; if(converged==1)break; - if(count>=1000){ - _printf_(" maximum number of iterations exceeded\n"); + if(count>=max_nonlinear_iterations){ + _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); break; } @@ -155,5 +164,5 @@ } - + /*Delete loads only if no ouput was requested: */ if(!input_loads)delete loads; Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 3200) +++ /issm/trunk/src/m/classes/@model/model.m (revision 3201) @@ -116,4 +116,5 @@ md.thermal_exchange_velocity=0; md.min_thermal_constraints=0; + md.min_mechanical_constraints=0; md.stabilize_constraints=0; @@ -165,4 +166,5 @@ md.eps_rel=0; md.eps_abs=0; + md.max_nonlinear_iterations=0; md.sparsity=0; md.connectivity=0; Index: /issm/trunk/src/m/classes/@model/setdefaultparameters.m =================================================================== --- /issm/trunk/src/m/classes/@model/setdefaultparameters.m (revision 3200) +++ /issm/trunk/src/m/classes/@model/setdefaultparameters.m (revision 3201) @@ -65,4 +65,7 @@ md.eps_abs=10; +%maximum of non-linear iterations. +md.max_nonlinear_iterations=100; + %sparsity md.sparsity=0.001; @@ -114,4 +117,5 @@ %parameter is often used. md.min_thermal_constraints=0; +md.min_mechanical_constraints=0; %Transient parameters Index: /issm/trunk/src/m/classes/public/display/displaydiagnostic.m =================================================================== --- /issm/trunk/src/m/classes/public/display/displaydiagnostic.m (revision 3200) +++ /issm/trunk/src/m/classes/public/display/displaydiagnostic.m (revision 3201) @@ -15,4 +15,5 @@ fielddisplay(md,'eps_rel','velocity relative convergence criterion, NaN -> not applied'); fielddisplay(md,'eps_abs','velocity absolute convergence criterion, NaN -> not applied'); +fielddisplay(md,'max_nonlinear_iterations','maximum number of nonlinear iterations'); fielddisplay(md,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'); @@ -23,4 +24,5 @@ disp(sprintf('\n %s','Penalties:')); fielddisplay(md,'penalty_offset','offset used by penalties: penalty = Kmax*10^offset'); +fielddisplay(md,'min_mechanical_constraints','threshold for instability of mechanical constraints'); disp(sprintf('\n %s','Memory management:')); Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 3200) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 3201) @@ -122,4 +122,5 @@ WriteData(fid,md.eps_rel,'Scalar','eps_rel'); WriteData(fid,md.eps_abs,'Scalar','eps_abs'); +WriteData(fid,md.max_nonlinear_iterations,'Scalar','max_nonlinear_iterations'); WriteData(fid,md.dt,'Scalar','dt'); WriteData(fid,md.ndt,'Scalar','ndt'); @@ -143,4 +144,5 @@ WriteData(fid,md.thermalconductivity,'Scalar','thermalconductivity'); WriteData(fid,md.min_thermal_constraints,'Integer','min_thermal_constraints'); +WriteData(fid,md.min_mechanical_constraints,'Integer','min_mechanical_constraints'); WriteData(fid,md.mixed_layer_capacity,'Scalar','mixed_layer_capacity'); WriteData(fid,md.thermal_exchange_velocity,'Scalar','thermal_exchange_velocity');