Changeset 4954

Timestamp:

08/03/10 20:59:13 (15 years ago)

Author:

Eric.Larour

Message:

Made fully constrained model runs possibles (ie, if f-set is NILL, the solution
still goes through).
This will allow testing of features in the code that do not rely on solutions, but
on specific modules.

Location:

issm/trunk

Files:

• src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp (modified) (1 diff)
• src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp (modified) (1 diff)
• src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.cpp (modified) (1 diff)
• src/c/modules/Solverx/Solverx.cpp (modified) (1 diff)
• src/c/solutions/convergence.cpp (modified) (1 diff)
• src/m/classes/public/ismodelselfconsistent.m (modified) (1 diff)
• src/m/solutions/convergence.m (modified) (1 diff)
• src/m/solvers/solver_diagnostic_nonlinear.m (modified) (2 diffs)
• test/Verification/test01_IceShelfIceFrontM2d/Square.par (modified) (2 diffs)

issm/trunk/src/c/modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp

@@ -21 +21 @@
         /*Merge f set back into n set: */
         u_n=NewVec(nodesets->GetNSize());
-        VecMerge(u_n,uf,nodesets->GetPV_F(),nodesets->GetFSize());
+        if(uf) VecMerge(u_n,uf,nodesets->GetPV_F(),nodesets->GetFSize());
 
         /*merge single point constraints into n set: */

issm/trunk/src/c/modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp

@@ -64 +64 @@
                 if (nodesets->GetSSize()){
 
-                        VecPartition(&pf, pn, nodesets->GetPV_F(),nodesets->GetFSize());
+                        if(nodesets->GetFSize()){
+                                VecPartition(&pf, pn, nodesets->GetPV_F(),nodesets->GetFSize());
 
-                        /*pf = pf - Kfs * y_s;*/
-                        MatGetLocalSize(Kfs,&Kfsm,&Kfsn);
-                        Kfsy_s=NewVecFromLocalSize(Kfsm);
-                        if (flag_ys0){
+                                /*pf = pf - Kfs * y_s;*/
+                                MatGetLocalSize(Kfs,&Kfsm,&Kfsn);
+                                Kfsy_s=NewVecFromLocalSize(Kfsm);
+                                if (flag_ys0){
 
-                                /*Create y_s0, full of 0: */
-                                VecDuplicate(y_s,&y_s0);
-                                VecSet(y_s0,0.0);
-                                VecAssemblyBegin(y_s0);
-                                VecAssemblyEnd(y_s0);
+                                        /*Create y_s0, full of 0: */
+                                        VecDuplicate(y_s,&y_s0);
+                                        VecSet(y_s0,0.0);
+                                        VecAssemblyBegin(y_s0);
+                                        VecAssemblyEnd(y_s0);
 
-                                MatMultPatch(Kfs,y_s0,Kfsy_s);
+                                        MatMultPatch(Kfs,y_s0,Kfsy_s);
+                                }
+                                else{
+                                        MatMultPatch(Kfs,y_s,Kfsy_s);
+                                }
+
+                                a=-1;
+                                VecAXPY(pf,a,Kfsy_s);  
                         }
                         else{
-                                MatMultPatch(Kfs,y_s,Kfsy_s);
+                                pf=NULL;
                         }
-
-                        a=-1;
-                        VecAXPY(pf,a,Kfsy_s);  
                 }
                 else{

issm/trunk/src/c/modules/Reducematrixfromgtofx/Reducematrixfromgtofx.cpp

@@ -19 +19 @@
         //Reduce matrix from n-size to f-size
         if(nodesets->GetSSize()){
-                
-                MatPartition(&Kff, Knn, nodesets->GetPV_F(),nodesets->GetFSize(),nodesets->GetPV_F(),nodesets->GetFSize()); 
-                MatPartition(&Kfs, Knn, nodesets->GetPV_F(),nodesets->GetFSize(),nodesets->GetPV_S(),nodesets->GetSSize()); 
 
+                if(nodesets->GetFSize()){
+                        MatPartition(&Kff, Knn, nodesets->GetPV_F(),nodesets->GetFSize(),nodesets->GetPV_F(),nodesets->GetFSize()); 
+                        MatPartition(&Kfs, Knn, nodesets->GetPV_F(),nodesets->GetFSize(),nodesets->GetPV_S(),nodesets->GetSSize()); 
+                }
         }
         else{

issm/trunk/src/c/modules/Solverx/Solverx.cpp

@@ -29 +29 @@
         PetscTruth flag;
         int solver_type;
+
+        /*First, check that f-set is not NULL, ie model is fully constrained: */
+        if(!Kff){
+                *puf=NULL;
+                return;
+        }
 
         /*If initial guess for solution exists, use it to create uf, otherwise, 

issm/trunk/src/c/solutions/convergence.cpp

@@ -31 +31 @@
         double yts;
 
-        /*Check arguments*/
-        ISSMASSERT(Kff && pf && uf && old_uf);
+        /*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from 
+         * the get go: */
+        if(!uf){
+                *pconverged=1;
+                return;
+        }
 
         /*get convergence options*/

issm/trunk/src/m/classes/public/ismodelselfconsistent.m

@@ -188 +188 @@
 %DIAGNOSTIC{{{1
 if md.analysis_type==DiagnosticSolutionEnum,
-
-        %CHECK THAT WE ARE NOT FULLY CONSTRAINED
-        if isempty(find(~md.spcvelocity(:,1:2))),
-                error(['model not consistent: model ' md.name ' is totally constrained horizontally, no need to solve!']);
-        end
-        if (md.dim==3),
-                if isempty(find(~md.spcvelocity(:,3))),
-                        error(['model not consistent: model ' md.name ' is totally constrained vertically, no need to solve!']);
-                end
-        end
 
         %HUTTER ON ICESHELF WARNING

issm/trunk/src/m/solutions/convergence.m

@@ -1 +1 @@
 function converged=convergence(K_ff,p_f,u_f,u_f_old,parameters)
+
+%before we start, check that u_f actually holds values, otherwise, we are done converging
+if isempty(u_f), converged=true;  return; end
 
 %Get convergence options

issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m

@@ -34 +34 @@
                 displaystring(femmodel.parameters.Verbose,'%s','   Generating penalty matrices:');
                 [K_gg , p_g, kmax]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
+                
 
                 displaystring(femmodel.parameters.Verbose,'%s','   reducing matrix from g to f set');
                 [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets); 
-
+                
                 displaystring(femmodel.parameters.Verbose,'%s','   reduce load from g size to f size');
                 p_f = Reduceloadfromgtof( p_g, femmodel.Gmn, K_fs, femmodel.ys, femmodel.nodesets);
-                
+
                 displaystring(femmodel.parameters.Verbose,'%s','   solving');
                 uf=Solver(K_ff,p_f,old_uf,femmodel.parameters);
@@ -54 +55 @@
                 
                 displaystring(femmodel.parameters.Verbose,'%s%i','      number of unstable constraints: ',num_unstable_constraints);
-                
+
                 %Figure out if convergence have been reached
                 converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);

issm/trunk/test/Verification/test01_IceShelfIceFrontM2d/Square.par

@@ -2 +2 @@
 
 %dynamics
-md.verbose=0;
+md.verbose=1;
 md.dt=1; %1 year
 md.ndt=md.dt*10; 
@@ -36 +36 @@
 md.time=50;
 md.waitonlock=30;
+
+%spc all velocities
+md.spcvelocity(:,1:3)=1;
+md.spcvelocity(:,4)=0;
+md.spcvelocity(:,5)=100;
+md.spcvelocity(:,6)=0;