KratosMultiphysics
KRATOS Multiphysics (Kratos) is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface.
pure_convection_tools.h
Go to the documentation of this file.
1 // KRATOS ___ ___ _ ___ __ ___ ___ ___ ___
2 // / __/ _ \| \| \ \ / /__| \_ _| __| __|
3 // | (_| (_) | .` |\ V /___| |) | || _|| _|
4 // \___\___/|_|\_| \_/ |___/___|_| |_| APPLICATION
5 //
6 // License: BSD License
7 // Kratos default license: kratos/license.txt
8 //
9 // Main authors: Riccardo Rossi
10 //
11 
12 #if !defined(KRATOS_PURE_CONVECTION_UTILITIES_INCLUDED )
13 #define KRATOS_PURE_CONVECTION_UTILITIES_INCLUDED
14 
15 
16 
17 // System includes
18 #include <string>
19 #include <iostream>
20 #include <algorithm>
21 
22 // External includes
23 
24 
25 // Project includes
26 #include "includes/define.h"
27 #include "includes/model_part.h"
28 #include "includes/node.h"
29 #include "utilities/geometry_utilities.h"
31 
32 namespace Kratos
33 {
34 
35 template<int TDim,class TSparseSpace,class TLinearSolver>
37 {
38 public:
39 
43 
44 
47 
48 
49 
50  //************************************************************************
51  //************************************************************************
53  {
55 
56  //loop over nodes with neighbours
57  mDofSet.clear(); // = DofsArrayType();
58  //mDofSet.resize(0);
59  mDofSet.reserve(model_part.Nodes().size() );
60 
61  int tot_nnz=0;
62  //int tot_row=0;
63  for (auto it=model_part.NodesBegin(); it!=model_part.NodesEnd(); ++it)
64  {
65  if( (it->GetValue(NEIGHBOUR_NODES)).size() != 0 )
66  {
67  mDofSet.push_back( it->pGetDof(rScalarVar) );
68  //tot_row += 1;
69  tot_nnz += (it->GetValue(NEIGHBOUR_NODES)).size()+1;
70  }
71  }
72 
73  //fills the DofList and give a unique progressive tag to each node
74 
75  int free_id = 0;
76  int fix_id = mDofSet.size();
77 
78  for (typename DofsArrayType::iterator dof_iterator = mDofSet.begin(); dof_iterator != mDofSet.end(); ++dof_iterator)
79  if (dof_iterator->IsFixed())
80  dof_iterator->SetEquationId(--fix_id);
81  else
82  dof_iterator->SetEquationId(free_id++);
83 
84  mEquationSystemSize = fix_id;
85 
86  std::vector< int > work_array;
87  work_array.reserve(1000);
88 
89  //guessing the total size of the matrix
90  mA.resize(mEquationSystemSize, mEquationSystemSize, tot_nnz);
91 
92  //getting the dof position
93  //unsigned int dof_position = (model_part.NodesBegin())->GetDofPosition(rScalarVar);
94 
95  //building up the matrix graph row by row
96  //int total_size = 0;
97  for (auto it=model_part.NodesBegin(); it!=model_part.NodesEnd(); ++it)
98  {
99  unsigned int index_i = it->GetDof(rScalarVar).EquationId();
100 
101  if(index_i < mEquationSystemSize && (it->GetValue(NEIGHBOUR_NODES)).size() != 0)
102  {
103  GlobalPointersVector< Node >& neighb_nodes = it->GetValue(NEIGHBOUR_NODES);
104 
105  //filling the first neighbours list
106  work_array.push_back(index_i);
107  for( GlobalPointersVector< Node >::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
108  {
109  unsigned int index_j = i->GetDof(rScalarVar).EquationId();
110  if(index_j < mEquationSystemSize)
111  work_array.push_back(index_j);
112  }
113 
114  //sorting the indices and elminating the duplicates
115  std::sort(work_array.begin(),work_array.end());
116  typename std::vector<int>::iterator new_end = std::unique(work_array.begin(),work_array.end());
117  unsigned int number_of_entries = new_end - work_array.begin();
118 
119  //filling up the matrix
120  for(unsigned int j=0; j<number_of_entries; j++)
121  {
122  mA.push_back(index_i,work_array[j] , 0.00);
123  }
124 
125  //clearing the array for the next step
126  work_array.erase(work_array.begin(),work_array.end());
127  //total_size += number_of_entries;
128  }
129 
130 
131 
132 
133  }
134 
135  mDx.resize(mA.size1(),false);
136  mb.resize(mA.size1(),false);
137  KRATOS_CATCH("")
138  }
139 
140  //************************************************************************
141  //************************************************************************
143  typename TLinearSolver::Pointer linear_solver,
144  Variable<double>& rScalarVar,
145  const Variable< array_1d<double,3> >& rTransportVel,
146  const Variable< array_1d<double,3> >& rMeshVel,
147  const Variable< double >& rProjVar,
148  unsigned int time_order )
149  {
150  KRATOS_TRY
151 
152  TSparseSpace::SetToZero(mA);
153  TSparseSpace::SetToZero(mb);
154  TSparseSpace::SetToZero(mDx);
155 
156  ProcessInfo& CurrentProcessInfo = model_part.GetProcessInfo();
157  double dt = CurrentProcessInfo[DELTA_TIME];
158 
159  Vector BDFcoeffs(time_order+1);
160 
161  if(time_order == 2)
162  {
163  if(model_part.GetBufferSize() < 3)
164  KRATOS_THROW_ERROR(std::logic_error,"insufficient buffer size for BDF2","")
165 
166  BDFcoeffs[0] = 1.5 / dt; //coefficient for step n+1
167  BDFcoeffs[1] = -2.0 / dt;//coefficient for step n
168  BDFcoeffs[2] = 0.5 / dt;//coefficient for step n-1
169  KRATOS_WATCH("CONVECTION: BDF_2.......................................................................................");
170  }
171  else
172  {
173  BDFcoeffs[0] = 1.0 / dt; //coefficient for step n+1
174  BDFcoeffs[1] = -1.0 / dt;//coefficient for step n
175  KRATOS_WATCH("CONVECTION: BDF_1.......................................................................................");
176  }
177 
178  //**********************************
179  //BUILD PHASE
180  //**********************************
182  BoundedMatrix<double,TDim+1,TDim+1> lhs_contribution;
184  array_1d<unsigned int ,TDim+1> local_indices;
185  array_1d<double,TDim+1> rhs_contribution;
186  array_1d<double,TDim+1> dp_vector;
187  array_1d<double,TDim+1> MassFactors;
188  array_1d<double,TDim> vel_gauss;
191  array_1d<double,TDim+1> temp_vec_np;
192 
193 
194 
195  for(unsigned int i=0; i<TDim+1; i++)
196  MassFactors[i] = 1.0/double(TDim+1);
197 
198  //getting the dof position
199 // unsigned int dof_position = (model_part.NodesBegin())->GetDofPosition(rScalarVar);
200 
201  double lumping_factor = 1.0/double(TDim+1.0);
202  const unsigned int number_of_points = TDim+1;
203 
205  for(ModelPart::ElementsContainerType::iterator i = model_part.ElementsBegin();
206  i!=model_part.ElementsEnd(); i++)
207  {
208  Geometry< Node >& geom = i->GetGeometry();
209 
210  //calculating elemental values
211  double Volume;
212  GeometryUtils::CalculateGeometryData(geom, DN_DX, N, Volume);
213 
214 
215  //finiding local indices
216  for(int ii = 0; ii<TDim+1; ii++)
217  {
218  local_indices[ii] = geom[ii].GetDof(rScalarVar).EquationId();
219 
220  /* dp_vector[ii] = geom[ii].FastGetSolutionStepValue(rScalarVar)
221  - geom[ii].FastGetSolutionStepValue(rScalarVar,1);*/
222  }
223 
224  //calculating convective velocity
225  noalias(vel_gauss) = ZeroVector(TDim);
226  double proj = 0.0;
227  for(unsigned int i = 0; i<number_of_points; i++)
228  {
229  const array_1d<double,3>& v = geom[i].FastGetSolutionStepValue(rTransportVel);
230  const array_1d<double,3>& w = geom[i].FastGetSolutionStepValue(rMeshVel);
231  proj += geom[i].FastGetSolutionStepValue(rProjVar);
232  for(unsigned int j = 0; j<TDim; j++)
233  {
234  vel_gauss[j] += v[j] - w[j];
235  }
236  }
237  vel_gauss *= lumping_factor;
238  proj *= lumping_factor;
239 
240 
241  //CONVECTIVE CONTRIBUTION TO THE STIFFNESS MATRIX
242  noalias(u_DN) = prod(DN_DX , vel_gauss);
243  noalias(lhs_contribution) = outer_prod(N,u_DN);
244 
245  //CONVECTION STABILIZING CONTRIBUTION (Suu)
246  double tau = CalculateTAU(vel_gauss,Volume);
247  noalias(lhs_contribution) += (tau) * outer_prod(u_DN,u_DN);
248 
249  //INERTIA CONTRIBUTION
250  for(unsigned int iii = 0; iii<number_of_points; iii++)
251  lhs_contribution(iii,iii) += BDFcoeffs[0] * MassFactors[iii];
252 
253  //adding projection contribution
254  //RHS += Suy * proj[component]
255  noalias(rhs_contribution) = (tau*proj)*u_DN;
256 
257 
258  //adding the inertia terms
259  // RHS += M*vhistory
260  //calculating the historical velocity
261  for(unsigned int iii = 0; iii<number_of_points; iii++)
262  {
263  temp_vec_np[iii] = BDFcoeffs[1]*geom[iii].FastGetSolutionStepValue(rScalarVar,1);
264  }
265  for(unsigned int step = 2; step<BDFcoeffs.size(); step++)
266  {
267  for(unsigned int iii = 0; iii<number_of_points; iii++)
268  temp_vec_np[iii] += BDFcoeffs[step]*geom[iii].FastGetSolutionStepValue(rScalarVar,step);
269  }
270  for(unsigned int iii = 0; iii<number_of_points; iii++)
271  rhs_contribution[iii] -= MassFactors[iii]*temp_vec_np[iii];
272 // rhs_contribution[iii] = - MassFactors[iii]*temp_vec_np[iii];
273 
274 
275  //subtracting the dirichlet term
276  // RHS -= LHS*rScalarVars
277  for(unsigned int iii = 0; iii<number_of_points; iii++)
278  temp_vec_np[iii] = geom[iii].FastGetSolutionStepValue(rScalarVar);
279  noalias(rhs_contribution) -= prod(lhs_contribution,temp_vec_np);
280 
281  //multiplying by Volume, rho and density
282  rhs_contribution *= Volume;
283  lhs_contribution *= Volume;
284 
285  AssembleLHS(mA,lhs_contribution,local_indices);
286  AssembleRHS(mb,rhs_contribution,local_indices);
287 
288 
289 
290  }
291 
292  //**********************************
293  //SOLVE PHASE
294  //**********************************
295  // KRATOS_WATCH(mA);
296  // KRATOS_WATCH(mb);
297  // KRATOS_WATCH(mDx);
298 // std::cout << "convection residual norm" << TSparseSpace::TwoNorm(mb) << std::endl;
299 
300  linear_solver->Solve(mA,mDx,mb);
301  std::cout << *(linear_solver) << std::endl;
302  // KRATOS_WATCH(mDx);
303 
304  //**********************************
305  //UPDATE rScalarVarS
306  //**********************************
307  for(typename DofsArrayType::iterator i_dof = mDofSet.begin() ; i_dof != mDofSet.end() ; ++i_dof)
308  {
309  if(i_dof->IsFree())
310  {
311  i_dof->GetSolutionStepValue() += mDx[i_dof->EquationId()];
312  }
313  }
314 
315  KRATOS_CATCH("")
316  }
317 
318 
319 
320 
321 
322  //************************************************************************
323  //************************************************************************
325  const Variable<double>& rScalarVar,
326  const Variable<double>& rNodalArea,
327  const Variable< array_1d<double,3> >& rTransportVel,
328  const Variable< array_1d<double,3> >& rMeshVel,
329  Variable< double >& rProjVar )
330  {
331  KRATOS_TRY
332 
333 
334  //**********************************
335  //BUILD PHASE
336  //**********************************
339  array_1d<double,TDim> vel_gauss;
340  array_1d<double,TDim+1> temp_vec_np;
342 
343  double lumping_factor = 1.0/double(TDim+1.0);
344  const unsigned int number_of_points = TDim+1;
345 
346  //set to zero nodal areas
347  for(ModelPart::NodesContainerType::iterator i = model_part.NodesBegin();
348  i!=model_part.NodesEnd(); i++)
349  {
350  i->FastGetSolutionStepValue(rNodalArea) = 0.0;
351  i->FastGetSolutionStepValue(rProjVar) = 0.0;
352  }
353 
354  //build up nodal areas and projections
356  for(ModelPart::ElementsContainerType::iterator i = model_part.ElementsBegin();
357  i!=model_part.ElementsEnd(); i++)
358  {
359 
360  Geometry< Node >& geom = i->GetGeometry();
361 
362  //calculating elemental values
363  double Volume;
364  GeometryUtils::CalculateGeometryData(geom, DN_DX, N, Volume);
365 
366  //calculating convective velocity
367  double nodal_area = Volume*lumping_factor;
368  noalias(vel_gauss) = ZeroVector(TDim);
369  for(unsigned int i = 0; i<number_of_points; i++)
370  {
371  const array_1d<double,3>& v = geom[i].FastGetSolutionStepValue(rTransportVel);
372  const array_1d<double,3>& w = geom[i].FastGetSolutionStepValue(rMeshVel);
373 
374  //adding up nodal area as needed
375  geom[i].FastGetSolutionStepValue(rNodalArea) += nodal_area;
376 
377  for(unsigned int j = 0; j<TDim; j++)
378  {
379  vel_gauss[j] += v[j] - w[j];
380  }
381  }
382  vel_gauss *= lumping_factor;
383 
384  //filling up a vector with all of the nodal values
385  for(unsigned int iii = 0; iii<number_of_points; iii++)
386  temp_vec_np[iii] = geom[iii].FastGetSolutionStepValue(rScalarVar);
387 
388  //
389  noalias(u_DN) = prod(DN_DX , vel_gauss);
390  double conv_proj = inner_prod( u_DN , temp_vec_np);
391  conv_proj *= Volume * lumping_factor;
392 
393  for(unsigned int iii = 0; iii<number_of_points; iii++)
394  geom[iii].FastGetSolutionStepValue(rProjVar) += conv_proj;
395  }
396 
397  //dividing by the overall mass matrix
398  for(ModelPart::NodesContainerType::iterator i = model_part.NodesBegin();
399  i!=model_part.NodesEnd(); i++)
400  {
401  i->FastGetSolutionStepValue(rProjVar) /= i->FastGetSolutionStepValue(rNodalArea);
402  }
403 
404  KRATOS_CATCH("")
405  }
406 
407  void ClearSystem()
408  {
409  KRATOS_TRY
410  mDofSet.clear(); // = DofsArrayType();
411  // mDofSet.resize(0);
412  mDx.resize(0,false);
413  mb.resize(0,false);
414  mA.resize(0,0,false);
415 
416  KRATOS_CATCH("")
417  }
418 
419 
420 
421 
422 private:
423 
424  unsigned int mEquationSystemSize;
425  TSystemVectorType mDx;
428  DofsArrayType mDofSet;
429 
430 
431 
432 
433 
434  //**************************************************************************
435  void AssembleLHS(
437  const BoundedMatrix<double,TDim+1,TDim+1>& LHS_Contribution,
438  const array_1d<unsigned int,TDim+1>& EquationId
439  )
440  {
441  unsigned int local_size = LHS_Contribution.size1();
442 
443  for (unsigned int i_local=0; i_local<local_size; i_local++)
444  {
445  unsigned int i_global=EquationId[i_local];
446  if ( i_global < mEquationSystemSize )
447  {
448  for (unsigned int j_local=0; j_local<local_size; j_local++)
449  {
450  unsigned int j_global=EquationId[j_local];
451  if ( j_global < mEquationSystemSize )
452  A(i_global,j_global) += LHS_Contribution(i_local,j_local);
453  }
454  }
455  }
456  }
457 
458 
459 
460  //**************************************************************************
461  void AssembleRHS(
463  const array_1d<double,TDim+1>& RHS_Contribution,
464  const array_1d<unsigned int,TDim+1>& EquationId
465  )
466  {
467  unsigned int local_size = RHS_Contribution.size();
468 
469  for (unsigned int i_local=0; i_local<local_size; i_local++)
470  {
471  unsigned int i_global=EquationId[i_local];
472  if ( i_global < mEquationSystemSize ) //on "free" DOFs
473  {
474  // ASSEMBLING THE SYSTEM VECTOR
475  b[i_global] += RHS_Contribution[i_local];
476  }
477  }
478 
479  }
480 
481  double CalculateTAU(const array_1d<double,2>& vel, const double& Volume)
482  {
483  //calculating parameter tau
484 // double c1 = 4.00;
485  double c2 = 2.00;
486  double h = sqrt(2.00*Volume);
487  double norm_u =norm_2(vel);
488  if(norm_u > 1e-15)
489  return h / ( c2*norm_u );
490  return 0.0;
491  }
492 
493  double CalculateTAU(const array_1d<double,3>& vel, const double& Volume)
494  {
495  //calculating parameter tau
496 // double c1 = 4.00;
497  double c2 = 2.00;
498  double h = pow(6.00*Volume,0.3333333);
499  double norm_u =norm_2(vel);
500  if(norm_u > 1e-15)
501  return h / ( c2*norm_u );
502  return 0.0;
503  }
504 
505 };
506 
507 } // namespace Kratos.
508 
509 #endif // KRATOS_PURE_CONVECTION_UTILITIES_INCLUDED defined
510 
511 
Geometry base class.
Definition: geometry.h:71
static void CalculateGeometryData(const GeometryType &rGeometry, BoundedMatrix< double, 4, 3 > &rDN_DX, array_1d< double, 4 > &rN, double &rVolume)
This function is designed to compute the shape function derivatives, shape functions and volume in 3D...
Definition: geometry_utilities.h:176
This class is a vector which stores global pointers.
Definition: global_pointers_vector.h:61
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: global_pointers_vector.h:79
void push_back(TPointerType x)
Definition: global_pointers_vector.h:322
iterator begin()
Definition: global_pointers_vector.h:221
iterator end()
Definition: global_pointers_vector.h:229
This object defines an indexed object.
Definition: indexed_object.h:54
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
A sorted associative container similar to an STL set, but uses a vector to store pointers to its data...
Definition: pointer_vector_set.h:72
void clear()
Clear the set, removing all elements.
Definition: pointer_vector_set.h:663
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: pointer_vector_set.h:95
size_type size() const
Returns the number of elements in the container.
Definition: pointer_vector_set.h:502
void push_back(TPointerType x)
Adds a pointer to the end of the set.
Definition: pointer_vector_set.h:544
iterator begin()
Returns an iterator pointing to the beginning of the container.
Definition: pointer_vector_set.h:278
void reserve(int reservedsize)
Reserves memory for a specified number of elements.
Definition: pointer_vector_set.h:733
iterator end()
Returns an iterator pointing to the end of the container.
Definition: pointer_vector_set.h:314
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Definition: pure_convection_tools.h:37
void ClearSystem()
Definition: pure_convection_tools.h:407
PureConvectionUtilities()
Definition: pure_convection_tools.h:45
void ConvectScalarVar(ModelPart &model_part, typename TLinearSolver::Pointer linear_solver, Variable< double > &rScalarVar, const Variable< array_1d< double, 3 > > &rTransportVel, const Variable< array_1d< double, 3 > > &rMeshVel, const Variable< double > &rProjVar, unsigned int time_order)
Definition: pure_convection_tools.h:142
void ConstructSystem(ModelPart &model_part, Variable< double > &rScalarVar, Variable< array_1d< double, 3 > > &rTransportVel, Variable< array_1d< double, 3 > > &rMeshVel)
Definition: pure_convection_tools.h:52
TSparseSpace::VectorType TSystemVectorType
Definition: pure_convection_tools.h:41
PointerVectorSet< Dof< double >, IndexedObject > DofsArrayType
Definition: pure_convection_tools.h:42
TSparseSpace::MatrixType TSystemMatrixType
Definition: pure_convection_tools.h:40
~PureConvectionUtilities()
Definition: pure_convection_tools.h:46
void CalculateProjection(ModelPart &model_part, const Variable< double > &rScalarVar, const Variable< double > &rNodalArea, const Variable< array_1d< double, 3 > > &rTransportVel, const Variable< array_1d< double, 3 > > &rMeshVel, Variable< double > &rProjVar)
Definition: pure_convection_tools.h:324
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
#define KRATOS_WATCH(variable)
Definition: define.h:806
#define KRATOS_TRY
Definition: define.h:109
dt
Definition: DEM_benchmarks.py:173
Vector VectorType
Definition: geometrical_transformation_utilities.h:56
Matrix MatrixType
Definition: geometrical_transformation_utilities.h:55
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
AMatrix::MatrixProductExpression< TExpression1Type, TExpression2Type > prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:568
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
AMatrix::VectorOuterProductExpression< TExpression1Type, TExpression2Type > outer_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:582
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
TABLE_NUMBER_ANGULAR_VELOCITY TABLE_NUMBER_MOMENT I33 BEAM_INERTIA_ROT_UNIT_LENGHT_Y KRATOS_DEFINE_APPLICATION_VARIABLE(DEM_APPLICATION, double, BEAM_INERTIA_ROT_UNIT_LENGHT_Z) typedef std double
Definition: DEM_application_variables.h:182
int time_order
Definition: ProjectParameters.py:51
int step
Definition: face_heat.py:88
model_part
Definition: face_heat.py:14
v
Definition: generate_convection_diffusion_explicit_element.py:114
w
Definition: generate_convection_diffusion_explicit_element.py:108
tau
Definition: generate_convection_diffusion_explicit_element.py:115
h
Definition: generate_droplet_dynamics.py:91
int local_size
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:17
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
vel
Definition: pure_conduction.py:76
int j
Definition: quadrature.py:648
A
Definition: sensitivityMatrix.py:70
N
Definition: sensitivityMatrix.py:29
integer i
Definition: TensorModule.f:17
e
Definition: run_cpp_mpi_tests.py:31