yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03589
[Branch ~yade-dev/yade/trunk] Rev 2076: - Introduced timing information from engines and functors
------------------------------------------------------------
revno: 2076
committer: Emanuele Catalano <ecatalano@r2balme>
branch nick: trunk
timestamp: Thu 2010-03-11 15:48:26 +0100
message:
- Introduced timing information from engines and functors
- Reorganized the code and cleaned from useless stuff
modified:
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.h
pkg/dem/Engine/PartialEngine/FlowEngine.cpp*
pkg/dem/Engine/PartialEngine/FlowEngine.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2010-03-01 17:55:59 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2010-03-11 14:48:26 +0000
@@ -28,7 +28,7 @@
const double RELAX = 1.9;
-const double TOLERANCE = 1e-06;
+
const double ONE_THIRD = 1.0/3.0;
//! Use this factor, or minLength, to reduce max permeability values (see usage below))
const double MAXK_DIV_KMEAN = 5000;
@@ -70,6 +70,7 @@
for (int i=0;i<6;i++) boundsIds[i] = 0;
minPermLength=-1;
SLIP_ON_LATERALS = true;//no-slip/symmetry conditions on lateral boundaries
+ TOLERANCE = 1e-06;
}
void FlowBoundingSphere::Compute_Action()
@@ -161,7 +162,8 @@
clock.top("DisplayStatistics");
/** START GAUSS SEIDEL */
// Boundary_Conditions ( Tri );
- Initialize_pressures();
+ double P_zero = abs((boundary(y_min_id).value-boundary(y_max_id).value)/2);
+ Initialize_pressures( P_zero );
GaussSeidel();
clock.top("GaussSeidel");
@@ -300,7 +302,7 @@
if (V->info().isFictious) {
++pass;
//FIXME : remove the isFictious flag and use cell->info().fictious() instead
- Boundary& bi = boundary(V->info().id());
+// Boundary& bi = boundary(V->info().id());
// Boundary& bi = boundaries [V->info().id()];
// if (bi.flowCondition) {
@@ -1507,15 +1509,13 @@
return rayon2 * fast_solid_angle(STA1,STA2,STA3,PTA1);
}
-void FlowBoundingSphere::Initialize_pressures()
+void FlowBoundingSphere::Initialize_pressures( double P_zero )
{
RTriangulation& Tri = T[currentTes].Triangulation();
- // Finite_cells_iterator cell_end = Tri.finite_cells_end();
- //
- // for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- // if (cell->info().fictious()) cell->info().p() = boundary(cell->info().id()).value;
- // else cell->info().p() =0;//FIXME : assign better values for faster convergence?
- // }
+ Finite_cells_iterator cell_end = Tri.finite_cells_end();
+
+ for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++)
+ {if (!cell->info().fictious()) cell->info().p() = P_zero;}
for (int bound=0; bound<6;bound++) {
int& id = *boundsIds[bound];
@@ -1858,7 +1858,8 @@
boundary(y_min_id).value=0;
boundary(y_max_id).value=1;
- Initialize_pressures();
+ double P_zero = abs((boundary(y_min_id).value-boundary(y_max_id).value)/2);
+ Initialize_pressures( P_zero );
GaussSeidel();
=== modified file 'lib/triangulation/FlowBoundingSphere.h'
--- lib/triangulation/FlowBoundingSphere.h 2010-03-01 17:55:59 +0000
+++ lib/triangulation/FlowBoundingSphere.h 2010-03-11 14:48:26 +0000
@@ -37,6 +37,7 @@
int* boundsIds [6];
bool currentTes;
bool SLIP_ON_LATERALS;
+ double TOLERANCE;
Boundary boundaries [6];
short id_offset;
@@ -75,7 +76,7 @@
// void Analytical_Consolidation ( );
void Boundary_Conditions ( RTriangulation& Tri );
- void Initialize_pressures ( );
+ void Initialize_pressures ( double P_zero );
/// Define forces using the same averaging volumes as for permeability
void ComputeTetrahedralForces();
void save_vtk_file ( RTriangulation &T );
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.cpp' (properties changed: +x to -x)
--- pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-03-01 17:55:59 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.cpp 2010-03-11 14:48:26 +0000
@@ -23,148 +23,80 @@
{
}
-std::ofstream plot ("pressurees.txt", std::ios::out);
-// std::ofstream cons_NONDAMP ("cons_NONDAMP_SLIP", std::ios::out);
-// std::ofstream settle_DAMP ("settle_DAMP_SLIP", std::ios::out);
-// std::ofstream settle_NONDAMP ("settle_NONDAMP_SLIP", std::ios::out);
-
void FlowEngine::applyCondition ( Scene* ncb )
{
- if (!flow) {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;}
+ if (!flow) {flow = shared_ptr<CGT::FlowBoundingSphere> (new CGT::FlowBoundingSphere);first=true;Update_Triangulation=false;}
if ( !isActivated ) return;
else
{
+ timingDeltas->start();
+
if ( !triaxialCompressionEngine )
- {
- vector<shared_ptr<Engine> >::iterator itFirst = ncb->engines.begin();
- vector<shared_ptr<Engine> >::iterator itLast = ncb->engines.end();
-
- for ( ;itFirst!=itLast; ++itFirst )
- {
- if ( ( *itFirst )->getClassName() == "TriaxialCompressionEngine" )
- {
- cout << "stress controller engine found - FlowEngine" << endl;
- triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> ( *itFirst );
- }
- }
- if ( !triaxialCompressionEngine ) cout << "stress controller engine NOT found" << endl;
- }
+ {vector<shared_ptr<Engine> >::iterator itFirst = ncb->engines.begin();vector<shared_ptr<Engine> >::iterator itLast = ncb->engines.end();for ( ;itFirst!=itLast; ++itFirst ){
+ if ( ( *itFirst )->getClassName() == "TriaxialCompressionEngine" ){
+ cout << "stress controller engine found - FlowEngine" << endl;
+ triaxialCompressionEngine = YADE_PTR_CAST<TriaxialCompressionEngine> ( *itFirst );}}
+ if ( !triaxialCompressionEngine ) cout << "stress controller engine NOT found" << endl;}
currentStress = triaxialCompressionEngine->stress[triaxialCompressionEngine->wall_top][1];
current_state = triaxialCompressionEngine->currentState;
- if ( !first && current_state==3 )
+ if ( current_state==3 )
{
- if (unload) {triaxialCompressionEngine->sigma_iso=(triaxialCompressionEngine->sigma_iso)/loadFactor;}
-
- UpdateVolumes ( ncb );
-
- cout << "simulation time = " << ncb->simulationTime << endl;
+ if ( first || Update_Triangulation ) { Build_Triangulation( ncb, P_zero );}
+// else
+// {
+ timingDeltas->checkpoint("Triangulating");
+
+ UpdateVolumes ( ncb );
+
+ timingDeltas->checkpoint("Update_Volumes");
///Compute flow and and forces here
- flow->GaussSeidel ( );
+ flow->GaussSeidel ( );
+
+ timingDeltas->checkpoint("Gauss-Seidel");
// flow->MGPost(flow->T[flow->currentTes].Triangulation());
- flow->tess_based_force=tess_based_force;
- flow->Compute_Forces ( );
+ flow->tess_based_force=tess_based_force;
+ flow->Compute_Forces ( );
+
+ timingDeltas->checkpoint("Compute_Forces");
///End Compute flow and forces
- CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end ();
-
- Vector3r f;
- int id;
-
- for ( CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin (); V_it != vertices_end; V_it++ )
- {
- id = V_it->info().id();
- for ( int y=0;y<3;y++ ) f[y] = ( V_it->info().forces ) [y];
-
- ncb->forces.addForce ( id, f );
+ CGT::Finite_vertices_iterator vertices_end = flow->T[flow->currentTes].Triangulation().finite_vertices_end ();
+ Vector3r f;int id;
+ for ( CGT::Finite_vertices_iterator V_it = flow->T[flow->currentTes].Triangulation().finite_vertices_begin (); V_it != vertices_end; V_it++ )
+ {
+ id = V_it->info().id();
+ for ( int y=0;y<3;y++ ) f[y] = ( V_it->info().forces ) [y];
+ ncb->forces.addForce ( id, f );
//ncb->forces.addTorque(id,t);
- }
-
- Real time = Omega::instance().getSimulationTime();
-
- int j = Omega::instance().getCurrentIteration();
-// int j = Omega::instance().getSimulationTime();
- char file [50];
- string consol = flow->key+"%d_Consol";
- const char* keyconsol = consol.c_str();
- sprintf (file, keyconsol, j);
- char *g = file;
-
- MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time);
- plot << j << " " << MaxPressure << endl;
-
-// if (damped) {cons_DAMP << j << " " << time << " " << flow->Pressures[cons] << endl; cons++;}
-// if (!damped){cons_NONDAMP << j << " " << time << " " << flow->Pressures[cons] << endl; cons++;}
-// if (damped) {settle_DAMP << j << " " << time << " " << triaxialCompressionEngine->uniaxialEpsilonCurr << endl;}
-// if (!damped) {settle_NONDAMP << j << " " << time << " " << triaxialCompressionEngine->uniaxialEpsilonCurr << endl;}
-
- if ( Omega::instance().getCurrentIteration() % PermuteInterval == 0 )
- {
-// K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max );
-
- cout << endl << "---NEW TRIANGULATION---" << endl << endl;
-
- NewTriangulation ( ncb );
-
- cout << "Vtotalissimo = " << flow->Vtotalissimo << " Vsolid_tot = " << flow->Vsolid_tot << " Vporale2 = " << flow->Vporale << " Ssolid_tot = " << flow->Ssolid_tot << endl << endl;
-
- flow->DisplayStatistics ();
-
- std::ofstream PFile ( "NewTriangulation_Pressures",std::ios::out );
-
- CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
- int j=0;
-
- for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
- {
- PFile << ++j << " " << cell->info().p() << endl;
}
-
- Initialize_volumes ( ncb );
- }
- }
- else if ( current_state==3 )
- {
- Initialize ( ncb, P_zero );
-
- flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
-
- flow->SLIP_ON_LATERALS=slip_boundary;
-
- flow->k_factor = permeability_factor;
- flow->Compute_Permeability ();
-
- flow->DisplayStatistics ();
-
- CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
- int y=0;
- for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
- {
- cell->info().dv() = 0; cell->info().p() = 0;
- y++;
- }
- cout << y << " deltaV initialised -----------------" << endl;
-
- flow->key = triaxialCompressionEngine->Key;
-
- if (compute_K) {K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
-
- Oedometer_Boundary_Conditions();
- flow->Initialize_pressures();
-
- flow->GaussSeidel ( );
-
-// plotFile << "unset key" << endl;
-
-
-// flow->Analytical_Consolidation();
- first = false; cons=0;
+ timingDeltas->checkpoint("Applying Forces");
+
+ Real time = Omega::instance().getSimulationTime();
+
+ int j = Omega::instance().getCurrentIteration();
+// int j = Omega::instance().getSimulationTime();
+ char file [50];
+ string consol = flow->key+"%d_Consol";
+ const char* keyconsol = consol.c_str();
+ sprintf (file, keyconsol, j);
+ char *g = file;
+ timingDeltas->checkpoint("Writing cons_files");
+
+ MaxPressure = flow->PermeameterCurve(flow->T[flow->currentTes].Triangulation(), g, time);
+
+ if ( Omega::instance().getCurrentIteration() % PermuteInterval == 0 )
+ { Update_Triangulation = true; }
+
+ timingDeltas->checkpoint("Storing Max Pressure");
+
+// }
}
}
}
@@ -186,11 +118,22 @@
triaxialCompressionEngine->sigma_iso=(triaxialCompressionEngine->sigma_iso)*loadFactor;
}
-void FlowEngine::Initialize ( Scene* ncb, double P_zero )
+void FlowEngine::Build_Triangulation ( Scene* ncb, double P_zero )
{
- flow->currentTes=0;
+ if (first)
+ {
+ flow->currentTes=0;
+ flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
+ flow->SLIP_ON_LATERALS=slip_boundary;
+ flow->key = triaxialCompressionEngine->Key;
+ flow->TOLERANCE=Tolerance;
+ }
+ else
+ {
+ if (compute_K) {K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
+ flow->currentTes=!flow->currentTes;
+ }
currentTes=flow->currentTes;
-
flow->x_min = 1000.0, flow->x_max = -10000.0, flow->y_min = 1000.0, flow->y_max = -10000.0, flow->z_min = 1000.0, flow->z_max = -10000.0;
AddBoundary ( ncb );
@@ -201,7 +144,26 @@
flow->Localize ();
+ flow->k_factor = permeability_factor;
+ flow->Compute_Permeability ();
+
+ if (first)
+ {
+ CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
+ for ( CGT::Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){cell->info().dv() = 0; cell->info().p() = 0;}
+ if (compute_K) { K = flow->Sample_Permeability ( flow->T[flow->currentTes].Triangulation(), flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max, flow->key );}
+ Oedometer_Boundary_Conditions();
+ first=!first;
+ }
+ else
+ {
+ flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
+ Update_Triangulation=!Update_Triangulation;
+ }
+
+ flow->Initialize_pressures( P_zero );
Initialize_volumes ( ncb );
+ flow->DisplayStatistics ();
}
void FlowEngine::AddBoundary ( Scene* ncb )
@@ -217,7 +179,7 @@
if ( b->shape->getClassIndex() == Bx_Index )
{
Box* w = YADE_CAST<Box*> ( b->shape.get() );
- const body_id_t& id = b->getId();
+// const body_id_t& id = b->getId();
Real center [3], Extent[3];
for ( int h=0;h<3;h++ ) {center[h] = b->state->pos[h]; Extent[h] = w->extents[h];}
@@ -236,8 +198,6 @@
void FlowEngine::Triangulate ( Scene* ncb )
{
- cout << "Triangulating------" << endl;
-
shared_ptr<Sphere> sph ( new Sphere );
int Sph_Index = sph->getClassIndexStatic();
@@ -260,8 +220,6 @@
contator+=1;
}
}
- cout << contator << "spheres inserted " << endl;
-
double SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min );
cout << "Section area = " << SectionArea << endl;
@@ -274,40 +232,6 @@
cout << "z_max = " << flow->z_max << endl;
}
-void FlowEngine::NewTriangulation ( Scene* ncb )
-{
- flow->x_min = 1000.0, flow->x_max = -10000.0, flow->y_min = 1000.0, flow->y_max = -10000.0, flow->z_min = 1000.0, flow->z_max = -10000.0;
-
- flow->currentTes=!flow->currentTes;
- currentTes = flow->currentTes;
-
- flow->T[flow->currentTes].Clear();
-
- Triangulate ( ncb );
-
- AddBoundary ( ncb );
-
-// flow->Tesselate();
- flow->T[flow->currentTes].Compute();
-
-// flow->Fictious_cells();
-
- flow->Localize ();
-
- flow->k_factor=permeability_factor;
-
- flow->Compute_Permeability ( );
-
-// flow->Sample_Permeability ( t2.Triangulation(), x_min, x_max, z_min, z_max, y_max, y_min );
-
- flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
-
-// flow->currentTes = !flow->currentTes;
-
-// flow->T[flow->currentTes] = flow->T[flow->currentTes];
-// t2 = flow->T[!flow->currentTes];
-}
-
void FlowEngine::Initialize_volumes ( Scene* ncb )
{
CGT::Finite_cells_iterator cell_end = flow->T[flow->currentTes].Triangulation().finite_cells_end();
=== modified file 'pkg/dem/Engine/PartialEngine/FlowEngine.hpp'
--- pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-03-01 17:55:59 +0000
+++ pkg/dem/Engine/PartialEngine/FlowEngine.hpp 2010-03-11 14:48:26 +0000
@@ -19,50 +19,44 @@
private:
shared_ptr<TriaxialCompressionEngine> triaxialCompressionEngine;
shared_ptr<CGT::FlowBoundingSphere> flow;
- public :
-
+ public :
Vector3r gravity;
-// bool first;
-
- int current_state
- ,previous_state
- ,cons;
-
+ int current_state;
Real wall_thickness;
+ bool Update_Triangulation;
void Triangulate ( Scene* ncb );
void AddBoundary ( Scene* ncb );
- void Initialize ( Scene* ncb, double P_zero );
+ void Build_Triangulation ( Scene* ncb, double P_zero );
void UpdateVolumes ( Scene* ncb );
void Initialize_volumes ( Scene* ncb );
Real Volume_cell_single_fictious (CGT::Cell_handle cell, Scene* ncb);
Real Volume_cell_double_fictious (CGT::Cell_handle cell, Scene* ncb);
Real Volume_cell_triple_fictious (CGT::Cell_handle cell, Scene* ncb);
Real Volume_cell (CGT::Cell_handle cell, Scene* ncb);
- void NewTriangulation ( Scene* ncb );
void Oedometer_Boundary_Conditions();
virtual ~FlowEngine();
virtual void applyCondition(Scene*);
- YADE_CLASS_BASE_DOC_ATTRS(FlowEngine,PartialEngine,"An engine to solve the flow problem in saturated granular media",
- ((bool,isActivated,true,"Activates Flow Engine "))
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(FlowEngine,PartialEngine,"An engine to solve flow problem in saturated granular media",
+ ((bool,isActivated,true,"Activates Flow Engine"))
((bool,first,true,"Controls the initialization/update phases"))
- ((bool, damped, true, ""))
((bool, slip_boundary, false, "Controls friction condition on lateral walls"))
((bool,currentTes,false,"Identifies the current triangulation/tesselation of pore space"))
((double,P_zero,0,"Initial internal pressure for oedometer test"))
+ ((double,Tolerance,1e-06,"Gauss-Seidel Tolerance"))
((int,PermuteInterval,100000,"Pore space re-triangulation period"))
((bool,compute_K,true,"Activates permeability measure within a granular sample"))
((double,permeability_factor,1.0,"a permability multiplicator"))
((Real,loadFactor,1.5,"Load multiplicator for oedometer test"))
- ((bool,unload,false,"Remove the load in oedometer test"))
((double, K, 0, "Permeability of the sample"))
((bool, ComputeFlow, 1,"If false only triangulation is done, if true flow problem is solved"))
((double, MaxPressure, 0, "Maximal value of water pressure within the sample"))
((double, currentStress, 0, "Current value of normal stress"))
- ((bool,tess_based_force,true,"true=force computation based on tessalation, false=force computation based on triangulation")));
+ ((bool,tess_based_force,true,"true=force computation based on tessalation, false=force computation based on triangulation")),
+ timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas));
DECLARE_LOGGER;
};
