Changeset 503

Timestamp:

05/19/09 15:38:40 (16 years ago)

Author:

Eric.Larour

Message:

Fixed all bugs for conduction, thermal steady solution

Location:

issm/trunk/src

Files:

• c/ConfigureObjectsx/ConfigureObjectsx.cpp (modified) (1 diff)
• c/DataSet/DataSet.cpp (modified) (1 diff)
• c/ModelProcessorx/CreateDataSets.cpp (modified) (2 diffs)
• c/ModelProcessorx/Melting/CreateParametersMelting.cpp (modified) (2 diffs)
• c/ModelProcessorx/Thermal/CreateParametersThermal.cpp (modified) (4 diffs)
• c/ProcessParamsx/ProcessParamsx.cpp (modified) (4 diffs)
• c/objects/Node.cpp (modified) (1 diff)
• c/objects/Pengrid.cpp (modified) (2 diffs)
• c/objects/Penta.cpp (modified) (4 diffs)
• c/objects/Tria.cpp (modified) (1 diff)
• m/solutions/cielo/thermal.m (modified) (3 diffs)
• m/solutions/cielo/thermal_core.m (modified) (2 diffs)
• mex/PenaltyConstraints/PenaltyConstraints.cpp (modified) (2 diffs)
• mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp (modified) (1 diff)

issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp

@@ -22 +22 @@
         
         /*intermediary: */
+
         elements->Configure(elements,loads,nodes,materials);
         loads->Configure(elements,loads,nodes,materials);

issm/trunk/src/c/DataSet/DataSet.cpp

@@ -1105 +1105 @@
 
                         pengrid=(Pengrid*)(*object);
+
                         pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type);
 

issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp

@@ -76 +76 @@
                 CreateConstraintsThermal(pconstraints,model,model_handle);
                 CreateLoadsThermal(ploads,model,model_handle);
+                CreateParametersThermal(pparameters,model,model_handle);
                                         
         }
@@ -83 +84 @@
                 CreateConstraintsMelting(pconstraints,model,model_handle);
                 CreateLoadsMelting(ploads,model,model_handle);
+                CreateParametersMelting(pparameters,model,model_handle);
         }
         else if (strcmp(model->analysis_type,"prognostic")==0){

issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp

@@ -1 @@
-/*!\file: CreateParametersDiagnosticHoriz.cpp
+/*!\file: CreateParametersMelting.cpp
  * \brief driver for creating parameters dataset, for diagnostic horiz analysis.
  */ 
 
 #undef __FUNCT__ 
-#define __FUNCT__ "CreateParametersDiagnosticHoriz"
+#define __FUNCT__ "CreateParametersMelting"
 
 #include "../../DataSet/DataSet.h"
@@ -13 +13 @@
 #include "../Model.h"
 
-void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
+void CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
         
         Param*   param = NULL;

issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp

@@ -1 @@
-/*!\file: CreateParametersDiagnosticHoriz.cpp
+/*!\file: CreateParametersThermal.cpp
  * \brief driver for creating parameters dataset, for diagnostic horiz analysis.
  */ 
 
 #undef __FUNCT__ 
-#define __FUNCT__ "CreateParametersDiagnosticHoriz"
+#define __FUNCT__ "CreateParametersThermal"
 
 #include "../../DataSet/DataSet.h"
@@ -13 +13 @@
 #include "../Model.h"
 
-void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
+void CreateParametersThermal(DataSet** pparameters,Model* model,ConstDataHandle model_handle){
         
         Param*   param = NULL;
@@ -23 +23 @@
         double* vy=NULL;
         double* vz=NULL;
+        double* pressure=NULL;
 
         /*recover parameters : */
@@ -63 +64 @@
         xfree((void**)&vz);
         
+        /*Get pressure: */
+        ModelFetchData((void**)&pressure,NULL,NULL,model_handle,"pressure","Matrix","Mat");
+        
+        count++;
+        param= new Param(count,"pressure",DOUBLEVEC);
+        if(pressure) param->SetDoubleVec(pressure,model->numberofnodes);
+        else param->SetDoubleVec(pressure,0);
+        parameters->AddObject(param);
+
+        /*Free pressure: */
+        xfree((void**)&pressure);
+
+        
         /*Assign output pointer: */
         *pparameters=parameters;

issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp

@@ -24 +24 @@
         double* vy=NULL;
         double* vz=NULL;
+        double* pressure=NULL;
         
         double* u_g=NULL;
@@ -54 +55 @@
 
         
-        if (   (analysis_type==ControlAnalysisEnum()) ||  (analysis_type==DiagnosticAnalysisEnum())){
+        if (   (analysis_type==ControlAnalysisEnum()) ||  (analysis_type==DiagnosticAnalysisEnum()) || (analysis_type==ThermalAnalysisEnum())){
 
                 parameters->FindParam((void*)&vx,"vx");
@@ -141 +142 @@
         }
 
+        if(analysis_type==ThermalAnalysisEnum()){
+
+                parameters->FindParam((void*)&pressure,"pressure");
+                
+                /*Now, from pressure, build p_g, correctly partitioned: */
+                p_g=(double*)xcalloc(numberofnodes,sizeof(double));
+
+                for(i=0;i<numberofnodes;i++){
+                        p_g[(int)(partition[i])]= pressure[i];  
+                }
+
+                /*Now, create new parameter: */
+                count++;
+                param= new Param(count,"p_g",DOUBLEVEC);
+                param->SetDoubleVec(p_g,numberofnodes);
+                parameters->AddObject(param);
+
+                /*erase old parameter: */
+                param=(Param*)parameters->FindParamObject("pressure");
+                parameters->DeleteObject((Object*)param);
+        }
+
         xfree((void**)&partition);
         
@@ -150 +173 @@
         xfree((void**)&vx_obs);
         xfree((void**)&vy_obs);
+        xfree((void**)&pressure);
         xfree((void**)&control_parameter);
         xfree((void**)&u_g_obs);

issm/trunk/src/c/objects/Node.cpp

@@ -444 +444 @@
         ParameterInputs* inputs=NULL;
         Node* node=this;
-        int* dof={0};
+        int dof[1]={0};
 
         /*Recover parameter inputs: */

issm/trunk/src/c/objects/Pengrid.cpp

@@ -529 +529 @@
         double meltingpoint;
         int    new_active;
-        int* dofs1={0};
+        int  dofs1[1]={0};
         int  unstable=0;
         
@@ -569 +569 @@
                 unstable=1;
         }
-        
+
         //Set penalty flag
         active=new_active;

issm/trunk/src/c/objects/Penta.cpp

@@ -3539 +3539 @@
         double         vz_list[numgrids];
         double         vxvyvz_list[numgrids][3];
-        double         pressure_list[3];
-        double         pressure;
         double         temperature_list[numgrids];
         double         temperature;
@@ -3599 +3597 @@
         if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!");
 
-        found=inputs->Recover("dt",&dt);
-        if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!");
-        
-        found=inputs->Recover("pressure",&pressure_list[0],1,dofs1,numgrids,(void**)nodes);
-        if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!");
-
-        found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes);
-        if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
-
+        if(sub_analysis_type==TransientAnalysisEnum()){
+                found=inputs->Recover("dt",&dt);
+                if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
+
+                found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes);
+                if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
+        }
+        
         for(i=0;i<numgrids;i++){
                 vx_list[i]=vxvyvz_list[i][0];
@@ -3662 +3659 @@
                         /* Build transient now */
                         if(sub_analysis_type==TransientAnalysisEnum()){
+                                GetParameterValue(&temperature, &temperature_list[0],gauss_coord);
                                 scalar_transient=temperature*Jdet*gauss_weight;
                                 for(i=0;i<numgrids;i++){
@@ -3676 +3674 @@
         if(onbed && shelf){
 
+
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
                 tria->CreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type);
                 delete tria;
-                return;
         }
 
         /* Geothermal flux on ice sheet base and basal friction */
         if(onbed && !shelf){
-                
+
                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
                 tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type);
                 delete tria;
-                return;
         }
 
         cleanup_and_return:
-        xfree((void**)first_gauss_area_coord);
-        xfree((void**)second_gauss_area_coord);
-        xfree((void**)third_gauss_area_coord);
-        xfree((void**)vert_gauss_coord);
-        xfree((void**)area_gauss_weights);
-        xfree((void**)vert_gauss_weights);
+        xfree((void**)&first_gauss_area_coord);
+        xfree((void**)&second_gauss_area_coord);
+        xfree((void**)&third_gauss_area_coord);
+        xfree((void**)&vert_gauss_coord);
+        xfree((void**)&area_gauss_weights);
+        xfree((void**)&vert_gauss_weights);
+
 }
 

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

@@ -2598 +2598 @@
 
         /*recover extra inputs from users, dt: */
-        inputs->Recover("dt",&dt);
+        found=inputs->Recover("dt",&dt);
         if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
 

issm/trunk/src/m/solutions/cielo/thermal.m

@@ -19 +19 @@
 
         %initialize inputs
-        displaystring(debug,'\n%s',['setup inputs...']);
+        displaystring(md.debug,'\n%s',['setup inputs...']);
         inputs=inputlist;
         inputs=add(inputs,'velocity',m_t.parameters.u_g,'doublevec',3,m_t.parameters.numberofnodes);
@@ -27 +27 @@
         if strcmpi(md.sub_analysis_type,'steady'),
         
-                displaystring(debug,'\n%s',['computing temperatures...']);
+                displaystring(md.debug,'\n%s',['computing temperatures...']);
                 [t_g m_t.loads melting_offset]=thermal_core(m_t,inputs,'thermal','steady');
                 
-                displaystring(debug,'\n%s',['computing melting...']);
+                displaystring(md.debug,'\n%s',['computing melting...']);
                 inputs=add(inputs,'melting_offset',melting_offset,'double');
                 inputs=add(inputs,'temperature',t_g,'doublevec',1,m_t.parameters.numberofnodes);
                 melting_g=diagnostic_core_linear(m_m,inputs,'melting','steady');
                 
-                displaystring(debug,'\n%s',['load results...']);
+                displaystring(md.debug,'\n%s',['load results...']);
                 md.temperature=t_g;
                 md.melting=melting_g*md.yts; %from m/s to m/a
@@ -53 +53 @@
                 for n=1:nsteps, 
 
-                        displaystring(debug,'\n%s%i/%i\n','time step: ',n,nsteps);
+                        displaystring(md.debug,'\n%s%i/%i\n','time step: ',n,nsteps);
                         
-                        displaystring(debug,'\n%s',['    computing temperatures...']);
+                        displaystring(md.debug,'\n%s',['    computing temperatures...']);
                         inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes);
                         [soln(n+1).t_g m_t.loads melting_offset]=thermal_core(m_t,inputs);

issm/trunk/src/m/solutions/cielo/thermal_core.m

@@ -26 +26 @@
 
                 %system matrices 
-                [K_gg_ , p_g_]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs);
-                [K_gg_ , p_g_]=SystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs);
+                displaystring(m.parameters.debug,'%s',['   system matrices']);
+                [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
+                displaystring(m.parameters.debug,'%s',['   penalty system matrices']);
+                [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
 
                 %Reduce tangent matrix from g size to f size
@@ -36 +38 @@
 
                 %Solve  
-                displaystring(m.parameters.debug,'\n%s',['   condition number of stiffness matrix: ',condest(K_ff)]);
+                displaystring(m.parameters.debug,'%s',['   condition number of stiffness matrix: ',condest(K_ff)]);
                 [t_f]=Solver(K_ff,p_f,[],m.parameters);
 
                 %Merge back to g set
+                displaystring(m.parameters.debug,'%s',['   merging solution back to g set']);
                 [t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); 
 
                 %Update inputs in datasets
                 inputs=add(inputs,'temperature',t_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes);
+                displaystring(m.parameters.debug,'%s',['   update inputs']);
                 [m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs);
         
                 %penalty constraints
-                [m.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, m.loads, m.materials,m.parameters,inputs);
+                displaystring(m.parameters.debug,'%s',['   penalty constraints']);
+                [m.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, m.loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type);
         
                 if ~converged,
                         displaystring(m.parameters.debug,'\n%s%i','#unstable constraints ',num_unstable_constraints);
+                        m.parameters.min_thermal_constraints
                         
-                        if num_unstable_constraints<m.parameters.min_thermal_constraints,
+                        if num_unstable_constraints<=m.parameters.min_thermal_constraints,
                                 converged=1;
                         end

issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp

@@ -36 +36 @@
         FetchData((void**)&loads,NULL,NULL,LOADSIN,"DataSet",NULL);
         FetchData((void**)&materials,NULL,NULL,MATERIALS,"DataSet",NULL);
-        
+
         /*parameters: */
         FetchData((void**)&analysis_type_string,NULL,NULL,ANALYSIS,"String",NULL);
@@ -73 +73 @@
 {
         _printf_("\n");
-        _printf_("   usage: [loads, constraints_converged, num_unstable_constraints] = %s(elements,nodes,loads,materials,params,inputs);\n",__FUNCT__);
+        _printf_("   usage: [loads, constraints_converged, num_unstable_constraints] = %s(elements,nodes,loads,materials,params,inputs,analysis_type,sub_analysis_type);\n",__FUNCT__);
         _printf_("\n");
 }

issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp

@@ -81 +81 @@
 {
         _printf_("\n");
-        _printf_("   usage: [Kgg,pg] = %s(Kggin,pgin,eleemnts,nodes,loads,materials,params,inputs);\n",__FUNCT__);
+        _printf_("   usage: [Kgg,pg] = %s(Kggin,pgin,eleemnts,nodes,loads,materials,params,inputs,analysis_type,sub_analysis_type);\n",__FUNCT__);
         _printf_("\n");
 }