Changeset 5698

Timestamp:

09/08/10 08:05:59 (15 years ago)

Author:

Mathieu Morlighem

Message:

Fixed new thermal

Location:

issm/trunk/src

Files:

• c/io/FetchData.cpp (modified) (4 diffs)
• c/modules/Solverx/Solverx.cpp (modified) (1 diff)
• c/solvers/solver_thermal_nonlinear.cpp (modified) (1 diff)
• m/classes/public/solversettomumps.m (modified) (1 diff)
• m/solvers/solver_thermal_nonlinear.m (modified) (1 diff)
• mex/Solver/Solver.cpp (modified) (5 diffs)
• mex/SystemMatrices/SystemMatrices.cpp (modified) (4 diffs)
• mex/SystemMatrices/SystemMatrices.h (modified) (2 diffs)

issm/trunk/src/c/io/FetchData.cpp

@@ -14 +14 @@
 
 #ifdef _SERIAL_
-/***************** **************** **************** **************** **************** **************** **************** **************** **************** ****************
-                                                                                                          Serial Fetch Data Routines, all overloaded.
-**************** **************** **************** **************** **************** **************** **************** **************** **************** *****************/
 #include <mex.h>
 /*FUNCTION FetchData(DataSet** pdataset,const mxArray* dataref){{{1*/
@@ -298 +295 @@
 void FetchData(bool* pboolean,const mxArray* dataref){
 
-        bool boolean;
+        bool* mxbool_ptr=NULL;
 
         if (mxIsClass(dataref,"logical")){
-                /*Recover the double: */
-                boolean=(bool)mxGetScalar(dataref);
+                if(mxGetM(dataref)!=1) ISSMERROR("input data is not of size 1x1");
+                if(mxGetN(dataref)!=1) ISSMERROR("input data is not of size 1x1");
+                mxbool_ptr=mxGetLogicals(dataref);
         }
         else{
@@ -308 +306 @@
         }
 
-        /*Assign output pointers:*/
-        *pboolean=boolean;
+        *pboolean=*mxbool_ptr;
 }
 /*}}}*/
@@ -315 +312 @@
 
 #if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
-/***************** **************** **************** **************** **************** **************** **************** **************** **************** ****************
-                                                                                                          Parallel Fetch Data Routines, all overloaded.
-**************** **************** **************** **************** **************** **************** **************** **************** **************** *****************/
 /*FUNCTION FetchData(double** pmatrix, int* pM,int* pN,FILE* fid){{{1*/
 void FetchData(double** pmatrix, int* pM,int* pN,FILE* fid){

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

@@ -98 +98 @@
                 }
         }
-        KSPSolve(ksp,pf,uf);
+        if(verbose){
+                double max,min;
+                KSPSetComputeSingularValues(ksp,PETSC_TRUE);
+                KSPSetUp(ksp);
+                KSPSolve(ksp,pf,uf);
+                KSPComputeExtremeSingularValues(ksp,&max,&min);
+                printf("Condition number = %g %g\n",max,min);
+        }
+        else{
+                KSPSolve(ksp,pf,uf);
+        }
         
         /*Check convergence*/

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

@@ -55 +55 @@
 
                 /*===================== DEBUGING ====================*/
-                if (count==1) SystemMatricesx(&Kgg_nopenalty,&pg_nopenalty,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+                if (count==1) SystemMatricesx(&Kgg_nopenalty,&pg_nopenalty,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,true,false,false);
                 MatDuplicate(Kgg_nopenalty,MAT_COPY_VALUES,&Kgg);
                 VecDuplicatePatch(&pg,pg_nopenalty);

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

@@ -5 +5 @@
 %      md=solversettomumps(md)
 
-md.solverstring=' -mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 120 -pc_factor_shift_positive_definite true';
+md.solverstring='-mat_type aijmumps -ksp_type preonly -pc_type lu -mat_mumps_icntl_14 120 -pc_factor_shift_positive_definite true';
 
 %optional

issm/trunk/src/m/solvers/solver_thermal_nonlinear.m

@@ -21 +21 @@
                 % ================= DEBUG FOR NOW ====================
                 if count==1 
-                        [K_gg_nopenalty,p_g_nopenalty,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+                        [K_gg_nopenalty,p_g_nopenalty,melting_offset]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,true,true,false,false);
                 end
                 [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);

issm/trunk/src/mex/Solver/Solver.cpp

@@ -8 +8 @@
 
         /*input datasets: */
-        Mat Kff=NULL;
-        Vec pf=NULL;
-        Vec uf0=NULL;
-        char* solver_string=NULL;
-        Parameters* parameters=NULL;
+        Mat         Kff           = NULL;
+        Vec         pf            = NULL;
+        Vec         uf0           = NULL;
+        Vec         uf            = NULL;
+        char       *solver_string = NULL;
+        Parameters *parameters    = NULL;
+        int         verbose;
                 
         /*Matlab solver: */
@@ -18 +20 @@
         char* matlabstring="-ksp_type matlab";
 
-        /* output datasets: */
-        Vec uf=NULL;
-
-        /*Boot module: */
         MODULEBOOT();
 
@@ -30 +28 @@
         FetchParams(&parameters,PARAMETERS);
         parameters->FindParam(&solver_string,SolverStringEnum);
+        int verbose; parameters->FindParam(&verbose,VerboseEnum);
         
         /*Fetch rest of data only if not running the matlab solver: */
@@ -44 +43 @@
         else{
                 /*Matlab solver: */
+                if (verbose) printf("   Solving\n");
                 array[0]=KFF;
                 array[1]=PF;
@@ -57 +57 @@
         xfree((void**)&solver_string);
 
-        /*end module: */
         MODULEEND();
 }

issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp

@@ -14 +14 @@
         Materials  *materials  = NULL;
         Parameters *parameters = NULL;
+        bool        kflag,pflag,penalty_kflag,penalty_pflag;
         
         /* output datasets: */
@@ -24 +25 @@
 
         /*checks on arguments on the matlab side: */
-        CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&SystemMatricesUsage);
+        if((nlhs!=NLHS) || (nrhs!=6 && nrhs!=10)){
+                SystemMatricesUsage();
+                ISSMERROR(" usage. See above");
+        }
 
         /*Input datasets: */
@@ -41 +45 @@
 
         /*!Generate internal degree of freedom numbers: */
-        SystemMatricesx(&Kgg,&pg,&kmax,elements,nodes,vertices,loads,materials,parameters);
+        if(nrhs==10){
+                FetchData(&kflag,KFLAG);
+                FetchData(&pflag,PFLAG);
+                FetchData(&penalty_kflag,PENALTYKFLAG);
+                FetchData(&penalty_pflag,PENALTYPFLAG);
+                SystemMatricesx(&Kgg,&pg,&kmax,elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);
+        }
+        else
+         SystemMatricesx(&Kgg,&pg,&kmax,elements,nodes,vertices,loads,materials,parameters);
 
         /*write output datasets: */
@@ -66 +78 @@
         _printf_("\n");
         _printf_("   usage: [Kgg,pg,kmax] = %s(elements,nodes,vertices,loads,materials,parameters);\n",__FUNCT__);
+        _printf_("   usage: [Kgg,pg,kmax] = %s(elements,nodes,vertices,loads,materials,parameters,kflag,pflag,penalty_kflag,penalty_pflag);\n",__FUNCT__);
         _printf_("\n");
 }

issm/trunk/src/mex/SystemMatrices/SystemMatrices.h

@@ -18 +18 @@
 
 /* serial input macros: */
-#define ELEMENTS   (mxArray *)prhs[0]
-#define NODES      (mxArray *)prhs[1]
-#define VERTICES   (mxArray *)prhs[2]
-#define LOADS      (mxArray *)prhs[3]
-#define MATERIALS  (mxArray *)prhs[4]
-#define PARAMETERS (mxArray *)prhs[5]
+#define ELEMENTS     (mxArray *)prhs[0]
+#define NODES        (mxArray *)prhs[1]
+#define VERTICES     (mxArray *)prhs[2]
+#define LOADS        (mxArray *)prhs[3]
+#define MATERIALS    (mxArray *)prhs[4]
+#define PARAMETERS   (mxArray *)prhs[5]
+#define KFLAG        (mxArray *)prhs[6]
+#define PFLAG        (mxArray *)prhs[7]
+#define PENALTYKFLAG (mxArray *)prhs[8]
+#define PENALTYPFLAG (mxArray *)prhs[9]
 
 /* serial output macros: */
@@ -34 +38 @@
 #define NLHS  3
 #undef NRHS
-#define NRHS  6
+#define NRHS  10
 
 #endif  /* _SYSTEMMATRICES_H */