Changeset 851

Timestamp:

06/08/09 16:45:56 (16 years ago)

Author:

Eric.Larour

Message:

New transient 3d parallel solution

Location:

issm/trunk/src

Files:

• c/Makefile.am (modified) (2 diffs)
• c/parallel/transient_core.cpp (modified) (1 diff)
• c/parallel/transient_core_2d.cpp (modified) (1 diff)
• c/parallel/transient_core_3d.cpp (modified) (2 diffs)
• m/solutions/cielo/transient3d.m (modified) (1 diff)

issm/trunk/src/c/Makefile.am

@@ -224 +224 @@
                                         ./UpdateGeometryx/UpdateGeometryx.h\
                                         ./UpdateGeometryx/UpdateGeometryx.cpp\
+                                        ./UpdateNodePositionsx/UpdateNodePositionsx.h\
+                                        ./UpdateNodePositionsx/UpdateNodePositionsx.cpp\
                                         ./ConfigureObjectsx/ConfigureObjectsx.h\
                                         ./ConfigureObjectsx/ConfigureObjectsx.cpp\
@@ -492 +494 @@
                                         ./UpdateGeometryx/UpdateGeometryx.h\
                                         ./UpdateGeometryx/UpdateGeometryx.cpp\
+                                        ./UpdateNodePositionsx/UpdateNodePositionsx.h\
+                                        ./UpdateNodePositionsx/UpdateNodePositionsx.cpp\
                                         ./ConfigureObjectsx/ConfigureObjectsx.h\
                                         ./ConfigureObjectsx/ConfigureObjectsx.cpp\

issm/trunk/src/c/parallel/transient_core.cpp

@@ -21 +21 @@
 
         int dim=-1;
-        
-        
+
+
         fem_p=fems+5;
 

issm/trunk/src/c/parallel/transient_core_2d.cpp

@@ -150 +150 @@
                 result=new Result(results->Size()+1,time,step,"b_g",b_g); results->AddObject(result);
         }
+
+        /*Free ressources:*/
+        xfree((void**)&u_g_serial);
+        xfree((void**)&p_g_serial);
+        xfree((void**)&h_g_serial);
+        xfree((void**)&s_g_serial);
+        xfree((void**)&b_g_serial);
 }

issm/trunk/src/c/parallel/transient_core_3d.cpp

@@ -4 +4 @@
 
 #undef __FUNCT__ 
-#define __FUNCT__ "transient_core_3d"
+#define __FUNCT__ "transient_3d_core"
 
 #include "../toolkits/toolkits.h"
@@ -15 +15 @@
 void transient_core_3d(DataSet* results,FemModel* fems, ParameterInputs* inputs){
 
-        throw ErrorException(__FUNCT__," 3d transient not supported yet!");
-
-
+        extern int my_rank;
+
+        /*fem models: */
+        FemModel* fem_dh=NULL;
+        FemModel* fem_dv=NULL;
+        FemModel* fem_dhu=NULL;
+        FemModel* fem_ds=NULL;
+        FemModel* fem_sl=NULL;
+        FemModel* fem_p=NULL;
+        FemModel* fem_t=NULL;
+        FemModel* fem_m=NULL;
+
+        /*output: */
+        Result* result=NULL;
+
+        /*solutions: */
+        int step;
+        double time;
+        
+        Vec u_g=NULL;
+        double* u_g_serial=NULL;
+        Vec p_g=NULL;
+        double* p_g_serial=NULL;
+        Vec h_g=NULL;
+        Vec h_g_intermediary=NULL;
+        double* h_g_serial=NULL;
+        Vec s_g=NULL;
+        double* s_g_serial=NULL;
+        Vec b_g=NULL;
+        double* b_g_serial=NULL;
+        Vec t_g=NULL;
+        Vec t_g_average=NULL;
+        double* t_g_serial=NULL;
+        Vec m_g=NULL;
+        double* m_g_serial=NULL;
+
+        /*intermediary: */
+        double finaltime;
+        double dt;
+        double yts;
+        DataSet* temperature_results=NULL;
+        DataSet* diagnostic_results=NULL;
+        DataSet* prognostic_results=NULL;
+
+
+        /*flags: */
+        int debug=0;
+        int numberofnodes;
+
+        /*dof recovery: */
+        int dof01[2]={0,1};
+        int dof2[1]={2};
+        int dof012[3]={0,1,2};
+        int dof3[1]={3};
+
+        /*recover fem models: */
+        fem_dh=fems+0;
+        fem_dv=fems+1;
+        fem_ds=fems+2;
+        fem_dhu=fems+3;
+        fem_sl=fems+4;
+        fem_p=fems+5;
+
+        fem_t=fems+6;
+        fem_m=fems+7;
+
+        //first recover parameters common to all solutions
+        fem_dh->parameters->FindParam((void*)&debug,"debug");
+        fem_dh->parameters->FindParam((void*)&finaltime,"ndt");
+        fem_dh->parameters->FindParam((void*)&dt,"dt");
+        fem_dh->parameters->FindParam((void*)&yts,"yts");
+        fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes");
+
+        /*initialize: */
+        step=1;
+        time=0;
+        
+        fem_p->parameters->FindParam((void*)&u_g_serial,"u_g");
+        u_g=SerialToVec(u_g_serial,2*numberofnodes); xfree((void**)&u_g_serial);
+
+        fem_p->parameters->FindParam((void*)&p_g_serial,"p_g");
+        p_g=SerialToVec(p_g_serial,1*numberofnodes); xfree((void**)&p_g_serial);
+
+        fem_p->parameters->FindParam((void*)&h_g_serial,"h_g");
+        h_g=SerialToVec(h_g_serial,1*numberofnodes); xfree((void**)&h_g_serial);
+
+        fem_p->parameters->FindParam((void*)&s_g_serial,"s_g");
+        s_g=SerialToVec(s_g_serial,1*numberofnodes); xfree((void**)&s_g_serial);
+
+        fem_p->parameters->FindParam((void*)&b_g_serial,"b_g");
+        b_g=SerialToVec(b_g_serial,1*numberofnodes); xfree((void**)&b_g_serial);
+
+        fem_p->parameters->FindParam((void*)&t_g_serial,"t_g");
+        t_g=SerialToVec(t_g_serial,1*numberofnodes); xfree((void**)&t_g_serial);
+
+        fem_p->parameters->FindParam((void*)&m_g_serial,"m_g");
+        m_g=SerialToVec(m_g_serial,1*numberofnodes); xfree((void**)&m_g_serial);
+
+        /*Add all these initial vectors in the results!: */
+        result=new Result(results->Size()+1,time,step,"u_g",u_g); results->AddObject(result);
+        result=new Result(results->Size()+1,time,step,"p_g",p_g); results->AddObject(result);
+        result=new Result(results->Size()+1,time,step,"h_g",h_g); results->AddObject(result);
+        result=new Result(results->Size()+1,time,step,"s_g",s_g); results->AddObject(result);
+        result=new Result(results->Size()+1,time,step,"b_g",b_g); results->AddObject(result);
+        result=new Result(results->Size()+1,time,step,"t_g",t_g); results->AddObject(result);
+        result=new Result(results->Size()+1,time,step,"m_g",m_g); results->AddObject(result);
+
+        while(time<finaltime+dt){ //make sure we run up to finaltime.
+        
+                if(debug)_printf_("%s%g%s%i%s%g\n","time [yr]: ",time/yts,"    iteration number: ",step,"/",floor(finaltime/dt));
+
+                step+=1;
+                time+=dt;
+
+                //update inputs
+                inputs->Add("thickness",h_g,1,numberofnodes);
+                inputs->Add("surface",s_g,1,numberofnodes);
+                inputs->Add("bed",b_g,1,numberofnodes);
+                inputs->Add("velocity",u_g,2,numberofnodes);
+                inputs->Add("pressure",p_g,1,numberofnodes);
+                inputs->Add("temperature",t_g,1,numberofnodes);
+
+                //Deal with temperature first 
+                if(debug)_printf_("%s\n","computing temperature");
+                temperature_results=new DataSet(ResultsEnum()); 
+                thermal_core(results,fems+6, inputs);
+                
+                //extract t_g and m_g from thermal, and erase thermal results;
+                temperature_results->FindResult(&t_g,"t_g");
+                temperature_results->FindResult(&m_g,"m_g");
+
+                //Compute depth averaged temperature and add to inputs
+                if(debug)_printf_("%s\n","computing depth average temperature");
+                VecDuplicatePatch(&t_g_average,t_g); 
+                FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->loads, fem_t->materials,"temperature");
+                inputs->Add("temperature_average",t_g_average,1,numberofnodes);
+                VecFree(&t_g_average); //not needed anymore
+
+                //Deal with velocities.
+                diagnostic_results=new DataSet(ResultsEnum()); 
+                diagnostic_core(diagnostic_results,fems, inputs);
+
+                //extract u_g and p_g from diagnostic_results, and erase diagnostic_results;
+                diagnostic_results->FindResult(&u_g,"u_g");
+                diagnostic_results->FindResult(&p_g,"p_g");
+                delete diagnostic_results;
+
+                //compute new thickness
+                if(debug)_printf_("%s\n","computing new thickness");
+                
+                inputs->Add("velocity",u_g,2,numberofnodes);
+                prognostic_results=new DataSet(ResultsEnum()); prognostic_core(prognostic_results,fem_p, inputs);
+        
+                //extract h_g prognostic_results, and erase prognostic_results;
+                prognostic_results->FindResult(&h_g_intermediary,"h_g");
+                delete prognostic_results;
+
+                //update surface and bed using the new thickness
+                if(debug)_printf_("   updating geometry\n");
+                UpdateGeometryx(&h_g,&b_g,&s_g, 
+                                fem_p->elements, fem_p->nodes,fem_p->loads, fem_p->materials, 
+                                h_g_intermediary,b_g,s_g); 
+                VecFree(&h_g_intermediary);
+                
+                if(debug)_printf_("%s\n","updating node positions");
+                UpdateNodePositionsx( fem_dh->elements,fem_dh->nodes,fem_dh->loads,fem_dh->materials,h_g,b_g);
+                UpdateNodePositionsx( fem_dv->elements,fem_dv->nodes,fem_dv->loads,fem_dv->materials,h_g,b_g);
+                UpdateNodePositionsx( fem_ds->elements,fem_ds->nodes,fem_ds->loads,fem_ds->materials,h_g,b_g);
+                UpdateNodePositionsx( fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,h_g,b_g);
+                UpdateNodePositionsx( fem_p->elements,fem_p->nodes,fem_p->loads,fem_p->materials,h_g,b_g);
+                UpdateNodePositionsx( fem_p->elements,fem_p->nodes,fem_p->loads,fem_p->materials,h_g,b_g);
+                UpdateNodePositionsx( fem_t->elements,fem_t->nodes,fem_t->loads,fem_t->materials,h_g,b_g);
+                UpdateNodePositionsx( fem_m->elements,fem_m->nodes,fem_m->loads,fem_m->materials,h_g,b_g);
+
+                //plug into results.
+                result=new Result(results->Size()+1,time,step,"u_g",u_g); results->AddObject(result);
+                result=new Result(results->Size()+1,time,step,"p_g",p_g); results->AddObject(result);
+                result=new Result(results->Size()+1,time,step,"h_g",h_g); results->AddObject(result);
+                result=new Result(results->Size()+1,time,step,"s_g",s_g); results->AddObject(result);
+                result=new Result(results->Size()+1,time,step,"b_g",b_g); results->AddObject(result);
+        }
+
+        /*Free ressources:*/
+        xfree((void**)&u_g_serial);
+        xfree((void**)&p_g_serial);
+        xfree((void**)&h_g_serial);
+        xfree((void**)&s_g_serial);
+        xfree((void**)&b_g_serial);
+        xfree((void**)&t_g_serial);
+        xfree((void**)&m_g_serial);
 }

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

@@ -115 +115 @@
         %end
 
-        %update nodes
+        %update node positions
         displaystring(md.debug,'   updating node positions...'));
         nodes_dh=UpdateNodePositions(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,new_bed,new_thickness);