Changeset 5997

Timestamp:

09/23/10 14:59:06 (14 years ago)

Author:

Mathieu Morlighem

Message:

merged the two solvers

Location:

issm/trunk/src/c

Files:

• Makefile.am (modified) (1 diff)
• solvers/solver_diagnostic_nonlinear.cpp (modified) (1 diff)
• solvers/solver_diagnostic_nonlinear_kff.cpp (deleted)
• solvers/solver_diagnostic_nonlinear_kgg.cpp (deleted)
• solvers/solvers.h (modified) (1 diff)

issm/trunk/src/c/Makefile.am

@@ -1146 +1146 @@
                                         ./solvers/solver_adjoint_linear.cpp\
                                         ./solvers/solver_diagnostic_nonlinear.cpp\
-                                        ./solvers/solver_diagnostic_nonlinear_kgg.cpp\
-                                        ./solvers/solver_diagnostic_nonlinear_kff.cpp\
                                         ./solvers/solver_thermal_nonlinear.cpp\
                                         ./modules/Bamgx/Bamgx.cpp\

issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp

@@ -12 +12 @@
 void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads){
 
+        /*intermediary: */
+        Mat Kgg = NULL, Kff = NULL, Kfs   = NULL;
+        Vec ug  = NULL, uf  = NULL, old_ug= NULL, old_uf = NULL;
+        Vec pg  = NULL, pf  = NULL;
+        Loads* loads=NULL;
+        int converged;
+        int constraints_converged;
+        int num_unstable_constraints;
+        int count;
+
         /*parameters:*/
+        int  verbose=0;
         bool kff=false;
+        int min_mechanical_constraints;
+        int max_nonlinear_iterations;
 
         /*Recover parameters: */
+        femmodel->parameters->FindParam(&verbose,VerboseEnum);
         femmodel->parameters->FindParam(&kff,KffEnum);
+        femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
+        femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum);
+        
+        /*Were loads requested as output? : */
+        if(conserve_loads) loads=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution
+        else               loads=(Loads*)femmodel->loads;         //modify loads  in this solution
 
-        if(!kff) solver_diagnostic_nonlinear_kgg(femmodel,conserve_loads);
-        else     solver_diagnostic_nonlinear_kff(femmodel,conserve_loads);
+        count=1;
+        converged=0;
 
+        /*Start non-linear iteration using input velocity: */
+        GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters);
+        Reducevectorgtofx(&uf, ug, femmodel->nodesets,femmodel->parameters);
+
+        //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate)
+        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
+
+        for(;;){
+
+                //save pointer to old velocity
+                VecFree(&old_ug);old_ug=ug;
+                VecFree(&old_uf);old_uf=uf;
+
+                if(kff){
+                        SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
+                        Reduceloadx(pf, Kfs, femmodel->ys, femmodel->parameters); MatFree(&Kfs);
+                }
+                else{
+                        SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
+                        Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
+                        Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters); VecFree(&pg); MatFree(&Kfs);
+                }
+                
+                Solverx(&uf, Kff, pf, old_uf, femmodel->parameters);
+
+                Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters);
+
+                InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
+
+                PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
+                if(verbose)_printf_("   number of unstable constraints: %i\n",num_unstable_constraints);
+
+                convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf);
+                
+                InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
+
+                //rift convergence
+                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>=max_nonlinear_iterations){
+                        _printf_("   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); 
+                        break;
+                }
+        }
+
+        /*clean-up*/
+        if(conserve_loads) delete loads;
+        VecFree(&uf);
+        VecFree(&ug);
+        VecFree(&old_uf);
+        VecFree(&old_ug);
 }

issm/trunk/src/c/solvers/solvers.h

@@ -14 +14 @@
 void solver_thermal_nonlinear(FemModel* femmodel);
 void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads);
-void solver_diagnostic_nonlinear_kgg(FemModel* femmodel,bool conserve_loads);
-void solver_diagnostic_nonlinear_kff(FemModel* femmodel,bool conserve_loads);
 void solver_linear(FemModel* femmodel);
 void solver_adjoint_linear(FemModel* femmodel);