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.
edgebased_levelset.h
Go to the documentation of this file.
1 // | / |
2 // ' / __| _` | __| _ \ __|
3 // . \ | ( | | ( |\__ `
4 // _|\_\_| \__,_|\__|\___/ ____/
5 // Multi-Physics
6 //
7 // License: BSD License
8 // Kratos default license: kratos/license.txt
9 //
10 // Main authors: Antonia Larese
11 //
12 
13 #if !defined(KRATOS_EDGEBASED_LEVELSET_FLUID_SOLVER_H_INCLUDED)
14 #define KRATOS_EDGEBASED_LEVELSET_FLUID_SOLVER_H_INCLUDED
15 
16 // System includes
17 #include <string>
18 #include <iostream>
19 #include <algorithm>
20 
21 // External includes
22 
23 // Project includes
24 #include "includes/define.h"
25 #include "includes/model_part.h"
26 #include "includes/deprecated_variables.h"
27 #include "includes/global_pointer_variables.h"
28 #include "includes/node.h"
29 #include "includes/cfd_variables.h"
30 //#include "geometries/geometry.h"
31 #include "utilities/geometry_utilities.h"
32 #include "free_surface_application.h"
33 
34 namespace Kratos
35 {
36 
37  template <unsigned int TDim, class MatrixContainer, class TSparseSpace, class TLinearSolver>
38  class EdgeBasedLevelSet
39  {
40  public:
41  // name for the self defined structure
42  typedef EdgesStructureType<TDim> CSR_Tuple;
43  typedef vector<CSR_Tuple> EdgesVectorType;
44 
45  // name for row start and column index vectors
46  typedef vector<unsigned int> IndicesVectorType;
47  // defining matrix type for test calculations
48  typedef vector<array_1d<double, TDim>> CalcVectorType;
49  // defining type for local storage of nodal values
50  typedef vector<double> ValuesVectorType;
51 
52  // defining types for matrix operations
53  typedef typename TSparseSpace::MatrixType TSystemMatrixType;
54  typedef typename TSparseSpace::VectorType TSystemVectorType;
55 
56  typedef std::size_t SizeType;
57 
58  // constructor and destructor
59 
60  EdgeBasedLevelSet(MatrixContainer &mr_matrix_container,
61  ModelPart &mr_model_part,
62  const double viscosity,
63  const double density,
64  bool use_mass_correction,
65  double stabdt_pressure_factor,
66  double stabdt_convection_factor,
67  double tau2_factor,
68  bool assume_constant_dp)
69  : mr_matrix_container(mr_matrix_container),
70  mr_model_part(mr_model_part),
71  mstabdt_pressure_factor(stabdt_pressure_factor),
72  mstabdt_convection_factor(stabdt_convection_factor),
73  mtau2_factor(tau2_factor),
74  massume_constant_dp(assume_constant_dp)
75 
76  {
77  for (ModelPart::NodesContainerType::iterator it = mr_model_part.NodesBegin(); it != mr_model_part.NodesEnd(); it++)
78  it->FastGetSolutionStepValue(VISCOSITY) = viscosity;
79 
80  mRho = density;
81 
82  mdelta_t_avg = 1000.0;
83 
84  max_dt = 1.0;
85 
86  muse_mass_correction = use_mass_correction;
87 
88  mshock_coeff = 0.7;
89  mWallLawIsActive = false;
90  };
91 
92  ~EdgeBasedLevelSet(){};
93 
94  //***********************************
95  // function to initialize fluid solver
96 
97  void Initialize()
98  {
99  KRATOS_TRY
100 
101  // get number of nodes
102  unsigned int n_nodes = mr_model_part.Nodes().size();
103  unsigned int n_edges = mr_matrix_container.GetNumberEdges();
104  // size data vectors
105  mBodyForce.resize(n_nodes);
106  mr_matrix_container.SetToZero(mBodyForce);
107  mViscosity.resize(n_nodes);
108  mr_matrix_container.SetToZero(mViscosity);
109  mWork.resize(n_nodes);
110  mr_matrix_container.SetToZero(mWork);
111  mvel_n.resize(n_nodes);
112  mr_matrix_container.SetToZero(mvel_n);
113  mvel_n1.resize(n_nodes);
114  mr_matrix_container.SetToZero(mvel_n1);
115  mPn.resize(n_nodes);
116  mr_matrix_container.SetToZero(mPn);
117  mPn1.resize(n_nodes);
118  mr_matrix_container.SetToZero(mPn1);
119  mHmin.resize(n_nodes);
120  mr_matrix_container.SetToZero(mHmin);
121  mHavg.resize(n_nodes);
122  mr_matrix_container.SetToZero(mHavg);
123  mNodalFlag.resize(n_nodes);
124  mr_matrix_container.SetToZero(mNodalFlag);
125  mdistances.resize(n_nodes);
126  mr_matrix_container.SetToZero(mdistances);
127 
128  mTauPressure.resize(n_nodes);
129  mr_matrix_container.SetToZero(mTauPressure);
130  mTauConvection.resize(n_nodes);
131  mr_matrix_container.SetToZero(mTauConvection);
132  mTau2.resize(n_nodes);
133  mr_matrix_container.SetToZero(mTau2);
134  mPi.resize(n_nodes);
135  mr_matrix_container.SetToZero(mPi);
136  mXi.resize(n_nodes);
137  mr_matrix_container.SetToZero(mXi);
138  mx.resize(n_nodes);
139  mr_matrix_container.SetToZero(mx);
140 
141  mEdgeDimensions.resize(n_edges);
142  mr_matrix_container.SetToZero(mEdgeDimensions);
143 
144  // convection variables
145  mBeta.resize(n_nodes);
146  mr_matrix_container.SetToZero(mBeta);
147  mPiConvection.resize(n_nodes);
148  mr_matrix_container.SetToZero(mPiConvection);
149  mphi_n.resize(n_nodes);
150  mr_matrix_container.SetToZero(mphi_n);
151  mphi_n1.resize(n_nodes);
152  mr_matrix_container.SetToZero(mphi_n1);
153 
154  mEps.resize(n_nodes);
155  mr_matrix_container.SetToZero(mEps);
156  // mD.resize(n_nodes); mr_matrix_container.SetToZero(mD);
157  mA.resize(n_nodes);
158  mr_matrix_container.SetToZero(mA);
159  mB.resize(n_nodes);
160  mr_matrix_container.SetToZero(mB);
161  mStrVel.resize(n_nodes);
162  mr_matrix_container.SetToZero(mStrVel);
163 
164  mdiv_error.resize(n_nodes);
165  mr_matrix_container.SetToZero(mdiv_error);
166  mdiag_stiffness.resize(n_nodes);
167  mr_matrix_container.SetToZero(mdiag_stiffness);
168  mis_slip.resize(n_nodes);
169 
170  // read velocity and pressure data from Kratos
171  mr_matrix_container.FillVectorFromDatabase(BODY_FORCE, mBodyForce, mr_model_part.Nodes());
172  mr_matrix_container.FillScalarFromDatabase(VISCOSITY, mViscosity, mr_model_part.Nodes());
173  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, mr_model_part.Nodes());
174  mr_matrix_container.FillScalarFromDatabase(PRESSURE, mPn1, mr_model_part.Nodes());
175  mr_matrix_container.FillOldScalarFromDatabase(PRESSURE, mPn, mr_model_part.Nodes());
176  mr_matrix_container.FillOldVectorFromDatabase(VELOCITY, mvel_n, mr_model_part.Nodes());
177  mr_matrix_container.FillCoordinatesFromDatabase(mx, mr_model_part.Nodes());
178 
179  // set flag for first time step
180  mFirstStep = true;
181 
182  // loop to categorize boundary nodes
183  std::vector<unsigned int> tempFixedVelocities;
184  std::vector<array_1d<double, TDim>> tempFixedVelocitiesValues;
185  std::vector<unsigned int> tempPressureOutletList;
186  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
187  inode != mr_model_part.NodesEnd();
188  inode++)
189  {
190  int index = inode->FastGetSolutionStepValue(AUX_INDEX);
191 
192  if (inode->Is(INLET))
193  {
194  tempFixedVelocities.push_back(index);
195  tempFixedVelocitiesValues.push_back(mvel_n1[index]);
196  }
197 
198  if (inode->Is(OUTLET) || inode->IsFixed(PRESSURE))
199  {
200  tempPressureOutletList.push_back(index);
201  }
202  }
203  mFixedVelocities.resize(tempFixedVelocities.size(), false);
204  mFixedVelocitiesValues.resize(tempFixedVelocitiesValues.size(), false);
205  mPressureOutletList.resize(tempPressureOutletList.size(), false);
206 
207  IndexPartition<unsigned int>(tempFixedVelocities.size()).for_each([&](unsigned int i){
208  mFixedVelocities[i] = tempFixedVelocities[i];
209  mFixedVelocitiesValues[i] = tempFixedVelocitiesValues[i];
210  });
211 
212  IndexPartition<unsigned int>(tempPressureOutletList.size()).for_each([&](unsigned int i){
213  mPressureOutletList[i] = tempPressureOutletList[i];
214  });
215 
216  // compute slip normals and fill SlipList
217  CalculateNormals(mr_model_part.Conditions());
218  mr_matrix_container.WriteVectorToDatabase(NORMAL, mSlipNormal, mr_model_part.Nodes());
219 
220  if constexpr (TDim == 3)
221  DetectEdges3D(mr_model_part.Conditions());
222 
223  // allocate memory for variables
224  mL.resize(n_nodes, n_nodes, false);
225 
226  int number_of_threads = ParallelUtilities::GetNumThreads();
227  std::vector<int> row_partition(number_of_threads);
228  OpenMPUtils::DivideInPartitions(n_nodes, number_of_threads, row_partition);
229 
230  for (int k = 0; k < number_of_threads; k++)
231  {
232  #pragma omp parallel
233  if (OpenMPUtils::ThisThread() == k)
234  {
235  for (int i_node = static_cast<int>(row_partition[k]); i_node < static_cast<int>(row_partition[k + 1]); i_node++)
236  {
237  // loop over all nodes
238  // flag for considering diagonal matrix elements
239  bool flag = 0;
240 
241  // loop over all neighbours
242  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
243  {
244  // get global index of neighbouring node j
245  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
246  // define matrix structure row by row (the order does matter!)
247  if ((static_cast<int>(j_neighbour) > i_node) && (flag == 0))
248  {
249  // add diagonal/nodal contribution
250  mL.push_back(i_node, i_node, 0.0);
251  flag = 1;
252  }
253  // add non-diagonal/edge contribution
254  mL.push_back(i_node, j_neighbour, 0.0);
255  }
256  // if diagonal element is the last non-zero element of the row
257  if (flag == 0)
258  mL.push_back(i_node, i_node, 0.0);
259  }
260  }
261  }
262 
263  // compute minimum length of the surrounding edges
264  CalculateEdgeLengths(mr_model_part.Nodes());
265 
266  // set the pressure projection to the body force value
267 
268  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
269  inode != mr_model_part.NodesEnd();
270  inode++)
271  {
272  array_1d<double, 3> temp = ZeroVector(3);
273  array_1d<double, 3> body_force = inode->FastGetSolutionStepValue(BODY_FORCE);
274  for (unsigned int i = 0; i < TDim; i++)
275  temp[i] = mRho * body_force[i];
276  array_1d<double, 3> &press_proj = inode->FastGetSolutionStepValue(PRESS_PROJ);
277  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
278  press_proj[l_comp] = temp[l_comp];
279  }
280 
281  KRATOS_CATCH("")
282  }
283 
284  void SetShockCapturingCoefficient(double coeff)
285  {
286  mshock_coeff = coeff;
287  }
288 
289  //***************************************
290  // function to set adequate time step size
291 
292  double ComputeTimeStep(const double CFLNumber, const double MaxDt)
293  {
294  KRATOS_TRY
295 
296  // save the maximum time step
297  max_dt = MaxDt;
298 
299  // local variable for time step size
300  double delta_t = 1e10; // max_dt;
301 
302  mdelta_t_avg = 1e10; // max_dt;
303 
304  // getting value of current velocity and of viscosity
305  mr_matrix_container.FillVectorFromDatabase(BODY_FORCE, mBodyForce, mr_model_part.Nodes());
306  mr_matrix_container.FillScalarFromDatabase(VISCOSITY, mViscosity, mr_model_part.Nodes());
307  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, mr_model_part.Nodes());
308  mr_matrix_container.FillScalarFromDatabase(POROSITY, mEps, mr_model_part.Nodes());
309  mr_matrix_container.FillScalarFromDatabase(LIN_DARCY_COEF, mA, mr_model_part.Nodes());
310  mr_matrix_container.FillScalarFromDatabase(NONLIN_DARCY_COEF, mB, mr_model_part.Nodes());
311  mr_matrix_container.FillVectorFromDatabase(STRUCTURE_VELOCITY, mStrVel, mr_model_part.Nodes());
312 
313  //*******************
314  // loop over all nodes
315  unsigned int n_nodes = mvel_n1.size();
316  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
317  {
318  const array_1d<double, TDim> &v_i = mvel_n1[i_node];
319  const double havg_i = mHavg[i_node];
320  const double hmin_i = mHmin[i_node];
321  const double eps_i = mEps[i_node];
322  const double nu = mViscosity[i_node];
323 
324  double vel_norm = 0.0;
325  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
326  {
327  vel_norm += v_i[l_comp] * v_i[l_comp];
328  }
329  vel_norm = sqrt(vel_norm);
330 
331  vel_norm /= eps_i;
332 
333  // use CFL condition to compute time step size
334  double delta_t_i = CFLNumber * 1.0 / (2.0 * vel_norm / hmin_i + 4.0 * nu / (hmin_i * hmin_i));
335  double delta_t_i_avg = 1.0 / (2.0 * vel_norm / havg_i + 4.0 * nu / (havg_i * havg_i));
336 
337  // considering the most restrictive case of neighbor's velocities with similar direction but opposite sense.
338  // loop over all neighbours
339  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
340  {
341  // get global index of neighbouring node j
342  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
343 
344  const array_1d<double, TDim> &v_j = mvel_n1[j_neighbour];
345 
346  double v_diff_norm = 0.0;
347  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
348  {
349  double temp = v_i[l_comp] - v_j[l_comp];
350  v_diff_norm += temp * temp;
351  }
352  v_diff_norm = sqrt(v_diff_norm);
353  v_diff_norm /= eps_i;
354  double delta_t_j = CFLNumber * 1.0 / (2.0 * v_diff_norm / hmin_i + 4.0 * nu / (hmin_i * hmin_i));
355 
356  if (delta_t_j < delta_t_i)
357  delta_t_i = delta_t_j;
358  }
359 
360  // choose the overall minimum of delta_t_i
361  if (delta_t_i < delta_t)
362  delta_t = delta_t_i;
363 
364  if (delta_t_i_avg < mdelta_t_avg)
365  mdelta_t_avg = delta_t_i_avg;
366  }
367 
368  return delta_t;
369 
370  KRATOS_CATCH("")
371  }
372 
373  void ApplySmagorinsky(double MolecularViscosity, double Cs)
374  {
375  if (Cs != 0)
376  {
377  if constexpr (TDim == 3)
378  ApplySmagorinsky3D(MolecularViscosity, Cs);
379  else
380  ApplySmagorinsky2D(MolecularViscosity, Cs);
381  }
382  }
383 
384  void UpdateFixedVelocityValues()
385  {
386  KRATOS_TRY
387 
388  // read velocity and pressure data from Kratos
389  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
390  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes);
391 
392  int fixed_size = mFixedVelocities.size();
393 
394  #pragma omp parallel for firstprivate(fixed_size)
395  for (int i_velocity = 0; i_velocity < fixed_size; i_velocity++)
396  {
397  unsigned int i_node = mFixedVelocities[i_velocity];
398  array_1d<double, TDim> &u_i_fix = mFixedVelocitiesValues[i_velocity];
399  const array_1d<double, TDim> &u_i = mvel_n1[i_node];
400 
401  for (unsigned int comp = 0; comp < TDim; comp++)
402  u_i_fix[comp] = u_i[comp];
403  }
404 
405  KRATOS_CATCH("");
406  }
407 
408  //**********************************************************************************
409  // function to solve fluid equations - fractional step 1: compute fractional momentum
410 
411  void SolveStep1()
412  {
413  KRATOS_TRY
414 
415  // PREREQUISITES
416 
417  // variables for node based data handling
418  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
419  int n_nodes = rNodes.size();
420  // storage of nodal values in local variables
421  CalcVectorType rhs;
422  rhs.resize(n_nodes);
423 
424  // read velocity and pressure data from Kratos
425  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes);
426  mr_matrix_container.FillOldVectorFromDatabase(VELOCITY, mvel_n, rNodes);
427  mr_matrix_container.FillVectorFromDatabase(BODY_FORCE, mBodyForce, rNodes);
428  mr_matrix_container.FillScalarFromDatabase(VISCOSITY, mViscosity, rNodes);
429  mr_matrix_container.FillScalarFromDatabase(PRESSURE, mPn1, rNodes);
430  mr_matrix_container.FillOldScalarFromDatabase(PRESSURE, mPn, rNodes);
431 
432  mr_matrix_container.FillScalarFromDatabase(DISTANCE, mdistances, mr_model_part.Nodes());
433  mr_matrix_container.FillScalarFromDatabase(POROSITY, mEps, mr_model_part.Nodes());
434  mr_matrix_container.FillScalarFromDatabase(LIN_DARCY_COEF, mA, mr_model_part.Nodes());
435  mr_matrix_container.FillScalarFromDatabase(NONLIN_DARCY_COEF, mB, mr_model_part.Nodes());
436  mr_matrix_container.FillVectorFromDatabase(STRUCTURE_VELOCITY, mStrVel, rNodes);
437 
438  // read time step size from Kratos
439  ProcessInfo &CurrentProcessInfo = mr_model_part.GetProcessInfo();
440  double delta_t = CurrentProcessInfo[DELTA_TIME];
441 
442  // compute intrinsic time
443  double time_inv_avg = 1.0 / mdelta_t_avg;
444 
445  double stabdt_pressure_factor = mstabdt_pressure_factor;
446  double stabdt_convection_factor = mstabdt_convection_factor;
447  double tau2_factor = mtau2_factor;
448 
449  #pragma omp parallel for firstprivate(time_inv_avg, stabdt_pressure_factor, stabdt_convection_factor, tau2_factor)
450  for (int i_node = 0; i_node < n_nodes; i_node++)
451  {
452  double &h_avg_i = mHavg[i_node];
453 
454  array_1d<double, TDim> &a_i = mvel_n1[i_node];
455  const double nu_i = mViscosity[i_node];
456  const double eps_i = mEps[i_node];
457  const double lindarcy_i = mA[i_node];
458  const double nonlindarcy_i = mB[i_node];
459 
460  double vel_norm = 0.0;
461  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
462  {
463  vel_norm += a_i[l_comp] * a_i[l_comp];
464  }
465  vel_norm = sqrt(vel_norm);
466 
467  const array_1d<double, TDim> &str_v_i = mStrVel[i_node];
468  array_1d<double, TDim> rel_vel_i;
469 
470  double rel_vel_norm = 0.0;
471  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
472  {
473  rel_vel_i[l_comp] = a_i[l_comp] - str_v_i[l_comp];
474  rel_vel_norm += rel_vel_i[l_comp] * rel_vel_i[l_comp];
475  }
476  rel_vel_norm = sqrt(rel_vel_norm);
477 
478  double porosity_coefficient = ComputePorosityCoefficient(rel_vel_norm, eps_i, lindarcy_i, nonlindarcy_i);
479 
480  vel_norm /= eps_i;
481 
482  double tau = 1.0 / (2.0 * vel_norm / h_avg_i + stabdt_pressure_factor * time_inv_avg + (4.0 * nu_i) / (h_avg_i * h_avg_i) + porosity_coefficient);
483  double tau_conv = 1.0 / (2.0 * vel_norm / h_avg_i + stabdt_convection_factor * time_inv_avg + (4.0 * nu_i) / (h_avg_i * h_avg_i) + porosity_coefficient);
484  mTauPressure[i_node] = tau;
485  mTauConvection[i_node] = tau_conv;
486 
487  mTau2[i_node] = (nu_i + h_avg_i * vel_norm * 0.5) * tau2_factor;
488  }
489 
490  // calculating the convective projection
491  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
492  array_1d<double, TDim> &pi_i = mPi[i_node];
493 
494  // setting to zero
495  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
496  pi_i[l_comp] = 0.0;
497 
498  array_1d<double, TDim> a_i = mvel_n1[i_node];
499  const array_1d<double, TDim> &U_i = mvel_n1[i_node];
500  const double &eps_i = mEps[i_node];
501 
502  a_i /= eps_i;
503 
504  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
505  {
506  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
507  array_1d<double, TDim> a_j = mvel_n1[j_neighbour];
508  const array_1d<double, TDim> &U_j = mvel_n1[j_neighbour];
509  const double &eps_j = mEps[j_neighbour];
510 
511  a_j /= eps_j;
512 
513  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
514 
515  edge_ij.Add_ConvectiveContribution(pi_i, a_i, U_i, a_j, U_j);
516  }
517 
518  const double m_inv = mr_matrix_container.GetInvertedMass()[i_node];
519 
520  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
521  pi_i[l_comp] *= m_inv;
522  });
523 
524  mr_matrix_container.AssignVectorToVector(mvel_n, mWork); // mWork = mvel_n
525  // first step of Runge Kutta
526  mr_matrix_container.AssignVectorToVector(mvel_n, mvel_n1); // mvel_n1 = mvel_n
527 
528  mr_matrix_container.SetToZero(rhs);
529  CalculateRHS(mvel_n1, mPn, mvel_n1, rhs, mdiag_stiffness);
530  Add_Effective_Inverse_Multiply(mWork, mWork, delta_t / 6.0, mr_matrix_container.GetLumpedMass(), mdiag_stiffness, rhs);
531  Add_Effective_Inverse_Multiply(mvel_n1, mvel_n, 0.5 * delta_t, mr_matrix_container.GetLumpedMass(), mdiag_stiffness, rhs);
532  ApplyVelocityBC(mvel_n1);
533 
534  // second step
535  mr_matrix_container.SetToZero(rhs);
536  CalculateRHS(mvel_n1, mPn, mvel_n1, rhs, mdiag_stiffness);
537  Add_Effective_Inverse_Multiply(mWork, mWork, delta_t / 3.0, mr_matrix_container.GetLumpedMass(), mdiag_stiffness, rhs);
538  Add_Effective_Inverse_Multiply(mvel_n1, mvel_n, 0.5 * delta_t, mr_matrix_container.GetLumpedMass(), mdiag_stiffness, rhs);
539  ApplyVelocityBC(mvel_n1);
540 
541  // third step
542  mr_matrix_container.SetToZero(rhs);
543  CalculateRHS(mvel_n1, mPn, mvel_n1, rhs, mdiag_stiffness);
544  Add_Effective_Inverse_Multiply(mWork, mWork, delta_t / 3.0, mr_matrix_container.GetLumpedMass(), mdiag_stiffness, rhs);
545  Add_Effective_Inverse_Multiply(mvel_n1, mvel_n, delta_t, mr_matrix_container.GetLumpedMass(), mdiag_stiffness, rhs);
546  ApplyVelocityBC(mvel_n1);
547 
548  // fourth step
549  mr_matrix_container.SetToZero(rhs);
550  CalculateRHS(mvel_n1, mPn, mvel_n1, rhs, mdiag_stiffness);
551  Add_Effective_Inverse_Multiply(mWork, mWork, delta_t / 6.0, mr_matrix_container.GetLumpedMass(), mdiag_stiffness, rhs);
552  // compute right-hand side
553  mr_matrix_container.AssignVectorToVector(mWork, mvel_n1);
554  ApplyVelocityBC(mvel_n1);
555 
556  KRATOS_CATCH("")
557  }
558 
559  //*********************************************************************
560  // function to calculate right-hand side of fractional momentum equation
561 
562  void CalculateRHS(
563  const CalcVectorType &vel,
564  const ValuesVectorType &pressure,
565  const CalcVectorType &convective_velocity,
566  CalcVectorType &rhs,
567  ValuesVectorType &diag_stiffness)
568  {
569  KRATOS_TRY
570 
571  int n_nodes = vel.size();
572 
573  // calculating the RHS
574  array_1d<double, TDim> stab_low;
575  array_1d<double, TDim> stab_high;
576 
577  double inverse_rho = 1.0 / mRho;
578 
579  #pragma omp parallel for private(stab_low, stab_high)
580  for (int i_node = 0; i_node < n_nodes; i_node++)
581  {
582  double dist = mdistances[i_node];
583  if (dist <= 0.0) // node is inside domain ---- if outside do nothing
584  {
585  const double nu_i = mViscosity[i_node];
586  const double nu_j = nu_i;
587  array_1d<double, TDim> &rhs_i = rhs[i_node];
588  const array_1d<double, TDim> &f_i = mBodyForce[i_node];
589  array_1d<double, TDim> a_i = convective_velocity[i_node];
590  const array_1d<double, TDim> &U_i = vel[i_node];
591  const array_1d<double, TDim> &pi_i = mPi[i_node];
592  const double &p_i = pressure[i_node];
593  const double &eps_i = mEps[i_node];
594  const double lindarcy_i = mA[i_node];
595  const double nonlindarcy_i = mB[i_node];
596 
597  const array_1d<double, TDim> &str_v_i = mStrVel[i_node];
598  array_1d<double, TDim> rel_vel_i;
599 
600  double rel_vel_norm = 0.0;
601  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
602  {
603  rel_vel_i[l_comp] = U_i[l_comp] - str_v_i[l_comp];
604  rel_vel_norm += rel_vel_i[l_comp] * rel_vel_i[l_comp];
605  }
606  rel_vel_norm = sqrt(rel_vel_norm);
607 
608  double edge_tau = mTauConvection[i_node];
609 
610  a_i /= eps_i;
611 
612  // initializing with the external forces (e.g. gravity)
613  double &m_i = mr_matrix_container.GetLumpedMass()[i_node];
614  for (unsigned int comp = 0; comp < TDim; comp++)
615  rhs_i[comp] = m_i * eps_i * f_i[comp];
616 
617  // applying the effect of the porosity
618  double porosity_coefficient = ComputePorosityCoefficient(rel_vel_norm, eps_i, lindarcy_i, nonlindarcy_i);
619  diag_stiffness[i_node] = m_i * porosity_coefficient;
620 
621  // /**************************************************rel_vel_modifications_b*/
622  for (unsigned int comp = 0; comp < TDim; comp++)
623  {
624  rhs_i[comp] += m_i * porosity_coefficient * str_v_i[comp];
625  }
626  // /*************************************************rel_vel_modifications_e*/
627 
628  // convective term
629  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
630  {
631  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
632  array_1d<double, TDim> a_j = convective_velocity[j_neighbour];
633  const array_1d<double, TDim> &U_j = vel[j_neighbour];
634  const array_1d<double, TDim> &pi_j = mPi[j_neighbour];
635  const double &p_j = pressure[j_neighbour];
636  const double &eps_j = mEps[j_neighbour];
637 
638  a_j /= eps_j;
639 
640  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
641 
642  edge_ij.Sub_ConvectiveContribution(rhs_i, a_i, U_i, a_j, U_j);
643 
644  // take care! we miss including a B.C. for the external pressure
645  edge_ij.Sub_grad_p(rhs_i, p_i * inverse_rho * eps_i, p_j * inverse_rho * eps_i);
646 
647  edge_ij.Sub_ViscousContribution(rhs_i, U_i, nu_i, U_j, nu_j);
648 
649  // add stabilization
650  edge_ij.CalculateConvectionStabilization_LOW(stab_low, a_i, U_i, a_j, U_j);
651  edge_ij.CalculateConvectionStabilization_HIGH(stab_high, a_i, pi_i, a_j, pi_j);
652  edge_ij.Sub_StabContribution(rhs_i, edge_tau, 1.0, stab_low, stab_high);
653  }
654  }
655  }
656 
657  // apply wall resistance
658  if (mWallLawIsActive == true)
659  ComputeWallResistance(vel, diag_stiffness);
660  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
661  mr_matrix_container.WriteVectorToDatabase(VELOCITY, mvel_n1, rNodes);
662 
663  KRATOS_CATCH("")
664  }
665 
666  //*************************************************************************
667  // function to solve fluid equations - fractional step 2: calculate pressure
668 
669  void SolveStep2(typename TLinearSolver::Pointer pLinearSolver)
670  {
671  KRATOS_TRY
672 
673  typedef Node PointType;
674  typedef GlobalPointersVector<PointType> PointVector;
675  typedef PointVector::iterator PointIterator;
676 
677  // reset is visited flag
678  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
679  inode != mr_model_part.NodesEnd();
680  inode++)
681  {
682  inode->GetValue(IS_VISITED) = 0.0;
683  }
684  // Re-generate a container with LAYER 0 and LAYER 1 after convection of the free surface
685  std::vector<PointVector> layers(2);
686 
687  // detect the nodes inside the fluid surface LAYER_0
688  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
689  inode != mr_model_part.NodesEnd();
690  inode++)
691  {
692  if (inode->FastGetSolutionStepValue(DISTANCE) < 0.0) // candidates are only the ones inside the fluid domain
693  {
694  GlobalPointersVector<Node> &neighb_nodes = inode->GetValue(NEIGHBOUR_NODES);
695  for (GlobalPointersVector<Node>::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
696  {
697  if (i->FastGetSolutionStepValue(DISTANCE) >= 0.0) // add the node as free surface if one of its neighb is outside
698  {
699  if (inode->GetValue(IS_VISITED) == 0.0)
700  {
701  layers[0].push_back(*(inode.base()));
702  inode->GetValue(IS_VISITED) = 1.0;
703  }
704  }
705  }
706  }
707  else
708  inode->FastGetSolutionStepValue(PRESSURE) = 0.0;
709  }
710 
711  // fill layer 1 by neighbour relationships
712  for (PointIterator iii = (layers[0]).begin(); iii != (layers[0]).end(); iii++)
713  {
714  GlobalPointersVector<Node> &neighb_nodes = iii->GetValue(NEIGHBOUR_NODES);
715  for (GlobalPointersVector<Node>::iterator jjj = neighb_nodes.begin(); jjj != neighb_nodes.end(); jjj++) // destination = origin1 + value * Minv*origin
716  {
717 
718  if (jjj->FastGetSolutionStepValue(DISTANCE) >= 0 &&
719  jjj->GetValue(IS_VISITED) == 0.0)
720  {
721  layers[1].push_back(Node::WeakPointer(*jjj.base()));
722  jjj->GetValue(IS_VISITED) = 2.0;
723  }
724  }
725  }
726 
727  // on the first layer outside the pressure is set to a value such that on the free surface the pressure is approx 0
728  for (PointIterator iii = layers[1].begin(); iii != layers[1].end(); iii++)
729  {
730  // get the node
731  unsigned int i_node = iii->FastGetSolutionStepValue(AUX_INDEX);
732 
733  array_1d<double, TDim> grad_d;
734  for (unsigned int comp = 0; comp < TDim; comp++)
735  grad_d[comp] = 0.0;
736 
737  double dist_i = mdistances[i_node];
738 
739  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
740  {
741  // get global index of neighbouring node j
742  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
743 
744  const double &dist_j = mdistances[j_neighbour];
745 
746  // projection of pressure gradients
747  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
748 
749  edge_ij.Add_grad_p(grad_d, dist_i, dist_j);
750  }
751 
752  const double &m_inv = mr_matrix_container.GetInvertedMass()[i_node];
753  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
754  grad_d[l_comp] *= m_inv;
755 
756  double norm_grad = norm_2(grad_d);
757 
758  if (norm_grad < 100.0)
759  {
760  grad_d /= norm_grad; // this is the direction of the gradient of the distances
761 
762  grad_d *= dist_i; // this is the vector with the distance of node_i from the closest point on the free surface
763 
764  double pestimate = 0.0;
765 
766  const array_1d<double, 3> &r_press_proj = iii->FastGetSolutionStepValue(PRESS_PROJ);
767  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
768  pestimate += r_press_proj[l_comp] * grad_d[l_comp];
769 
770  iii->FastGetSolutionStepValue(PRESSURE) = pestimate;
771  }
772  else
773  {
774  std::cout << "attention gradient of distance much greater than 1 on node:" << i_node << std::endl;
775  double avg_number = 0.0;
776 
777  double pavg = 0.0;
778 
779  GlobalPointersVector<Node> &neighb_nodes = iii->GetValue(NEIGHBOUR_NODES);
780  for (GlobalPointersVector<Node>::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
781  {
782  if (i->GetValue(IS_VISITED) == 1.0)
783  {
784  pavg += i->FastGetSolutionStepValue(PRESSURE);
785  avg_number += 1.0;
786  }
787  }
788 
789  KRATOS_ERROR_IF(avg_number == 0) << "can not happen that the extrapolation node has no neighbours" << std::endl;
790 
791  iii->FastGetSolutionStepValue(PRESSURE) = pavg / avg_number;
792  }
793  }
794 
795  // PREREQUISITES
796 
797  // allocate memory for variables
798  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
799  int n_nodes = rNodes.size();
800  // unknown and right-hand side vector
801  TSystemVectorType dp, rhs;
802  dp.resize(n_nodes, false);
803  rhs.resize(n_nodes, false);
804  array_1d<double, TDim> dU_i, dU_j, work_array;
805  // read time step size from Kratos
806  ProcessInfo &CurrentProcessInfo = mr_model_part.GetProcessInfo();
807  double delta_t = CurrentProcessInfo[DELTA_TIME];
808 
809  mr_matrix_container.FillOldScalarFromDatabase(PRESSURE, mPn, mr_model_part.Nodes());
810  mr_matrix_container.FillScalarFromDatabase(PRESSURE, mPn1, mr_model_part.Nodes());
811  mr_matrix_container.FillVectorFromDatabase(PRESS_PROJ, mXi, rNodes);
812  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes);
813 
814  // loop over all nodes
815  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
816 
817  double &rhs_i = rhs[i_node];
818  rhs_i = 0.0;
819  const double &p_i = mPn1[i_node];
820  const double &p_old_i = mPn[i_node];
821  const array_1d<double, TDim> &U_i_curr = mvel_n1[i_node];
822 
823  array_1d<double, TDim> &xi_i = mXi[i_node];
824 
825  double l_ii = 0.0;
826 
827  // loop over all neighbours
828  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
829  {
830  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
831  const double &p_j = mPn1[j_neighbour];
832  const double &p_old_j = mPn[j_neighbour];
833  const array_1d<double, TDim> &U_j_curr = mvel_n1[j_neighbour];
834  const array_1d<double, TDim> &xi_j = mXi[j_neighbour];
835 
836  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
837 
838 #ifdef SYMM_PRESS
839  double edge_tau = 0.25 * (mTauPressure[i_node] + mTauPressure[j_neighbour]);
840 #else
841  double edge_tau = 0.5 * mTauPressure[i_node];
842 #endif
843 
844  if (edge_tau < delta_t)
845  edge_tau = delta_t;
846 
847  // compute laplacian operator
848  double sum_l_ikjk;
849  edge_ij.CalculateScalarLaplacian(sum_l_ikjk);
850 
851  double sum_l_ikjk_onlydt = sum_l_ikjk * (delta_t);
852  sum_l_ikjk *= (delta_t + edge_tau);
853 
854  // assemble right-hand side
855  // pressure contribution
856  rhs_i -= sum_l_ikjk * (p_j - p_i);
857  rhs_i += sum_l_ikjk_onlydt * (p_old_j - p_old_i);
858 
859  // calculating the divergence of the fract vel
860  edge_ij.Sub_D_v(rhs_i, U_i_curr * mRho, U_j_curr * mRho);
861 
862  // high order stabilizing term
863  double temp = 0.0;
864 
865  edge_ij.Add_div_v(temp, xi_i, xi_j);
866  rhs_i += edge_tau * temp;
867 
868  // assemble laplacian matrix
869  mL(i_node, j_neighbour) = sum_l_ikjk;
870  l_ii -= sum_l_ikjk;
871  }
872 
873  mL(i_node, i_node) = l_ii;
874  });
875 
876  if (muse_mass_correction == true)
877  {
878  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
879  double &rhs_i = rhs[i_node];
880  rhs_i -= mdiv_error[i_node];
881  });
882  }
883 
884  // find the max diagonal term
885  double max_diag = 0.0;
886  for (int i_node = 0; i_node < n_nodes; i_node++)
887  {
888  double L_diag = mL(i_node, i_node);
889  if (std::abs(L_diag) > std::abs(max_diag))
890  max_diag = L_diag;
891  }
892  if (max_diag < 1e15)
893  max_diag = 1e15;
894 
895  // respect pressure boundary conditions by penalization
896  for (unsigned int i_pressure = 0; i_pressure < mPressureOutletList.size(); i_pressure++)
897  {
898  unsigned int i_node = mPressureOutletList[i_pressure];
899  mL(i_node, i_node) = max_diag;
900  rhs[i_node] = 0.0;
901  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
902  {
903  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
904  mL(i_node, j_neighbour) = 0.0;
905  }
906  }
907 
908  // modification for level_set
909  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
910  if (mdistances[i_node] >= 0)
911  {
912  mL(i_node, i_node) = max_diag;
913  rhs[i_node] = 0.0;
914  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
915  {
916  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
917  mL(i_node, j_neighbour) = 0.0;
918  }
919  }
920  else
921  {
922  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
923  {
924  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
925  if (mdistances[j_neighbour] >= 0)
926  mL(i_node, j_neighbour) = 0.0;
927  }
928  }
929  });
930 
931  // compute row scaling factors
932  TSystemVectorType scaling_factors(n_nodes);
933  double *Lvalues = mL.value_data().begin();
934  SizeType *Lrow_indices = mL.index1_data().begin();
935  SizeType *Lcol_indices = mL.index2_data().begin();
936 
937  IndexPartition<unsigned int>(mL.size1()).for_each([&](unsigned int k){
938  double t = 0.0;
939  SizeType col_begin = Lrow_indices[k];
940  SizeType col_end = Lrow_indices[k + 1];
941 
942  for (SizeType j = col_begin; j < col_end; j++){
943  if (static_cast<unsigned int>(Lcol_indices[j]) == k)
944  {
945  t = fabs(Lvalues[j]);
946  break;
947  }
948  }
949  scaling_factors[k] = 1.0 / sqrt(t);
950  });
951 
952  IndexPartition<unsigned int>(mL.size1()).for_each([&](unsigned int k){
953  SizeType col_begin = Lrow_indices[k];
954  SizeType col_end = Lrow_indices[k + 1];
955  double k_factor = scaling_factors[k];
956 
957  rhs[k] *= k_factor;
958 
959  for (SizeType j = col_begin; j < col_end; j++)
960  {
961  Lvalues[j] *= scaling_factors[Lcol_indices[j]] * k_factor;
962  }
963  });
964 
965  // set starting vector for iterative solvers
966  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
967  dp[i_node] = 0.0;
968  });
969 
970  pLinearSolver->Solve(mL, dp, rhs);
971 
972  // update pressure
973  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
974  mPn1[i_node] += dp[i_node] * scaling_factors[i_node];
975  });
976 
977  // write pressure to Kratos
978  mr_matrix_container.WriteScalarToDatabase(PRESSURE, mPn1, rNodes);
979 
980  // compute pressure proj for the next step
981  #pragma omp parallel for private(work_array)
982  for (int i_node = 0; i_node < n_nodes; i_node++)
983  {
984  array_1d<double, TDim> &xi_i = mXi[i_node];
985  for (unsigned int comp = 0; comp < TDim; comp++)
986  xi_i[comp] = 0.0;
987 
988  double dist = mdistances[i_node];
989  if (dist <= 0.0) // node is inside domain ---- if outside do nothing
990  {
991 
992  const double &p_i = mPn1[i_node];
993 
994  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
995  {
996  // get global index of neighbouring node j
997  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
998 
999  const double &p_j = mPn1[j_neighbour];
1000 
1001  // projection of pressure gradients
1002  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
1003 
1004  edge_ij.Add_grad_p(xi_i, p_i, p_j);
1005  }
1006 
1007  const double &m_inv = mr_matrix_container.GetInvertedMass()[i_node];
1008  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
1009  xi_i[l_comp] *= m_inv;
1010  }
1011  }
1012 
1013  mr_matrix_container.WriteVectorToDatabase(PRESS_PROJ, mXi, rNodes);
1014 
1015  KRATOS_CATCH("")
1016  }
1017 
1018  //**********************************************************************************
1019  // function to solve fluid equations - fractional step 3: correct fractional momentum
1020 
1021  void SolveStep3()
1022  {
1023  KRATOS_TRY
1024  // get number of nodes
1025  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
1026 
1027  int n_nodes = rNodes.size();
1028 
1029  // define work array
1030  array_1d<double, TDim> correction;
1031  // read time step size from Kratos
1032  ProcessInfo &CurrentProcessInfo = mr_model_part.GetProcessInfo();
1033  double delta_t = CurrentProcessInfo[DELTA_TIME];
1034 
1035  double factor = 0.5;
1036  if (massume_constant_dp == true)
1037  factor = 1.0;
1038 
1039  // compute end of step momentum
1040  double rho_inv = 1.0 / mRho;
1041  #pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv, factor)
1042  for (int i_node = 0; i_node < n_nodes; i_node++)
1043  {
1044  double dist = mdistances[i_node];
1045  if (dist < 0.0) // node is inside domain ---- if outside do nothing
1046  {
1047  array_1d<double, TDim> &U_i_curr = mvel_n1[i_node];
1048  double delta_p_i = (mPn1[i_node] - mPn[i_node]) * rho_inv * factor;
1049 
1050  // setting to zero
1051  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
1052  correction[l_comp] = 0.0;
1053 
1054  // compute edge contributions dt*M^(-1)Gp
1055  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
1056  {
1057  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
1058  double delta_p_j = (mPn1[j_neighbour] - mPn[j_neighbour]) * rho_inv * factor;
1059 
1060  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
1061 
1062  edge_ij.Sub_grad_p(correction, delta_p_i, delta_p_j);
1063  }
1064  // compute prefactor
1065  const double m = mr_matrix_container.GetLumpedMass()[i_node];
1066  const double &d = mdiag_stiffness[i_node];
1067 
1068  // correct fractional momentum
1069  for (unsigned int comp = 0; comp < TDim; comp++)
1070  {
1071  U_i_curr[comp] += delta_t / (m + delta_t * d) * correction[comp];
1072  }
1073  }
1074  }
1075 
1076  ApplyVelocityBC(mvel_n1);
1077 
1078  // write velocity of time step n+1 to Kratos
1079  mr_matrix_container.WriteVectorToDatabase(VELOCITY, mvel_n1, rNodes);
1080 
1081  // calculate the error on the divergence
1082  if (muse_mass_correction == true)
1083  {
1084  #pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv)
1085  for (int i_node = 0; i_node < n_nodes; i_node++)
1086  {
1087  const double dist = mdistances[i_node];
1088  double &div_i_err = mdiv_error[i_node];
1089  div_i_err = 0.0;
1090  if (dist < 0.0) // node is inside domain ---- if outside do nothing
1091  {
1092  const array_1d<double, TDim> &U_i_curr = mvel_n1[i_node];
1093 
1094  // compute edge contributions dt*M^(-1)Gp
1095  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
1096  {
1097  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
1098  array_1d<double, TDim> &U_j_curr = mvel_n1[j_neighbour];
1099 
1100  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
1101 
1102  edge_ij.Add_D_v(div_i_err, U_i_curr * mRho, U_j_curr * mRho);
1103  }
1104  }
1105  }
1106  }
1107 
1108  KRATOS_CATCH("")
1109  }
1110 
1111  //************************************
1112 
1113  void ApplyVelocityBC(CalcVectorType &VelArray)
1114  {
1115  KRATOS_TRY
1116 
1117  if (mWallLawIsActive == false)
1118  {
1119  // apply conditions on corner edges
1120  int edge_size = medge_nodes_direction.size();
1121 
1122  #pragma omp parallel for firstprivate(edge_size)
1123  for (int i = 0; i < edge_size; i++)
1124  {
1125  int i_node = medge_nodes[i];
1126  const array_1d<double, TDim> &direction = medge_nodes_direction[i];
1127  double dist = mdistances[i_node];
1128 
1129  if (dist <= 0.0)
1130  {
1131  array_1d<double, TDim> &U_i = VelArray[i_node];
1132  double temp = 0.0;
1133  for (unsigned int comp = 0; comp < TDim; comp++)
1134  temp += U_i[comp] * direction[comp];
1135 
1136  for (unsigned int comp = 0; comp < TDim; comp++)
1137  U_i[comp] = direction[comp] * temp;
1138  }
1139  }
1140 
1141  // apply conditions on corners
1142  int corner_size = mcorner_nodes.size();
1143  for (int i = 0; i < corner_size; i++)
1144  {
1145  int i_node = mcorner_nodes[i];
1146 
1147  array_1d<double, TDim> &U_i = VelArray[i_node];
1148  for (unsigned int comp = 0; comp < TDim; comp++)
1149  U_i[comp] = 0.0;
1150  }
1151  }
1152 
1153  // slip condition
1154  int slip_size = mSlipBoundaryList.size();
1155 
1156  #pragma omp parallel for firstprivate(slip_size)
1157  for (int i_slip = 0; i_slip < slip_size; i_slip++)
1158  {
1159  unsigned int i_node = mSlipBoundaryList[i_slip];
1160  double dist = mdistances[i_node];
1161  if (dist <= 0.0)
1162  {
1163  array_1d<double, TDim> &U_i = VelArray[i_node];
1164  array_1d<double, TDim> &an_i = mSlipNormal[i_node];
1165  double projection_length = 0.0;
1166  double normalization = 0.0;
1167  for (unsigned int comp = 0; comp < TDim; comp++)
1168  {
1169  projection_length += U_i[comp] * an_i[comp];
1170  normalization += an_i[comp] * an_i[comp];
1171  }
1172  projection_length /= normalization;
1173  // tangential momentum as difference between original and normal momentum
1174  for (unsigned int comp = 0; comp < TDim; comp++)
1175  U_i[comp] -= projection_length * an_i[comp];
1176  }
1177  }
1178 
1179  // fixed condition
1180  int fixed_size = mFixedVelocities.size();
1181 
1182  #pragma omp parallel for firstprivate(fixed_size)
1183  for (int i_velocity = 0; i_velocity < fixed_size; i_velocity++)
1184  {
1185  unsigned int i_node = mFixedVelocities[i_velocity];
1186  double dist = mdistances[i_node];
1187  if (dist <= 0.0)
1188  {
1189  const array_1d<double, TDim> &u_i_fix = mFixedVelocitiesValues[i_velocity];
1190  array_1d<double, TDim> &u_i = VelArray[i_node];
1191 
1192  for (unsigned int comp = 0; comp < TDim; comp++)
1193  u_i[comp] = u_i_fix[comp];
1194  }
1195  }
1196 
1197  KRATOS_CATCH("")
1198  }
1199 
1200  //********************************
1201  // function to compute coefficients
1202 
1203  void ExtrapolateValues(unsigned int extrapolation_layers)
1204  {
1205  KRATOS_TRY
1206 
1207  // ensure that corner nodes are wet if all of the nodes around them have a negative distance
1208 
1209  typedef Node PointType;
1210  typedef GlobalPointersVector<PointType> PointVector;
1211  typedef PointVector::iterator PointIterator;
1212 
1213  mr_matrix_container.FillScalarFromDatabase(DISTANCE, mdistances, mr_model_part.Nodes());
1214 
1215  // reset is visited flag
1216  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1217  inode != mr_model_part.NodesEnd();
1218  inode++)
1219  {
1220  inode->GetValue(IS_VISITED) = 0.0;
1221  }
1222 
1223  // generate a container with the layers to be extrapolated
1224  std::vector<PointVector> layers(extrapolation_layers);
1225 
1226  // detect the nodes inside the fluid surface
1227  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1228  inode != mr_model_part.NodesEnd();
1229  inode++)
1230  {
1231  if (inode->FastGetSolutionStepValue(DISTANCE) < 0.0) // candidates are only the ones inside the fluid domain
1232  {
1233  GlobalPointersVector<Node> &neighb_nodes = inode->GetValue(NEIGHBOUR_NODES);
1234  for (GlobalPointersVector<Node>::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
1235  {
1236  if (i->FastGetSolutionStepValue(DISTANCE) >= 0.0) // add the node as free surface if one of its neighb is outside
1237  {
1238  if (inode->GetValue(IS_VISITED) == 0.0)
1239  {
1240  layers[0].push_back(*(inode.base()));
1241  inode->GetValue(IS_VISITED) = 1.0;
1242  }
1243  }
1244  }
1245  }
1246  else
1247  {
1248  // set everything to zero
1249  noalias(inode->FastGetSolutionStepValue(VELOCITY)) = ZeroVector(3);
1250  inode->FastGetSolutionStepValue(PRESSURE) = 0.0;
1251  noalias(inode->FastGetSolutionStepValue(VELOCITY, 1)) = ZeroVector(3);
1252  inode->FastGetSolutionStepValue(PRESSURE, 1) = 0.0;
1253 
1254  noalias(inode->FastGetSolutionStepValue(PRESS_PROJ)) = ZeroVector(3);
1255  noalias(inode->FastGetSolutionStepValue(PRESS_PROJ, 1)) = ZeroVector(3);
1256  }
1257  }
1258 
1259  // fill the following layers by neighbour relationships
1260  // each layer fills the following
1261  for (unsigned int il = 0; il < extrapolation_layers - 1; il++)
1262  {
1263  for (PointIterator iii = (layers[il]).begin(); iii != (layers[il]).end(); iii++)
1264  {
1265  GlobalPointersVector<Node> &neighb_nodes = iii->GetValue(NEIGHBOUR_NODES);
1266  for (GlobalPointersVector<Node>::iterator jjj = neighb_nodes.begin(); jjj != neighb_nodes.end(); jjj++) // destination = origin1 + value * Minv*origin
1267  {
1268 
1269  if (jjj->FastGetSolutionStepValue(DISTANCE) >= 0 &&
1270  jjj->GetValue(IS_VISITED) == 0.0)
1271  {
1272  layers[il + 1].push_back(Node::WeakPointer(*jjj.base()));
1273  jjj->GetValue(IS_VISITED) = double(il + 2.0);
1274  }
1275  }
1276  }
1277  }
1278 
1279  array_1d<double, 3> aux, aux_proj;
1280 
1281  // TESTING!!!
1282  // fill the pressure projection on the first layer inside the fluid
1283  // by extrapolating from the pressure projection on the layer -1 (the first layer completely inside the domain)
1284  for (PointIterator iii = (layers[0]).begin(); iii != (layers[0]).end(); iii++)
1285  {
1286  noalias(aux_proj) = ZeroVector(3);
1287  double avg_number = 0.0;
1288 
1289  GlobalPointersVector<Node> &neighb_nodes = iii->GetValue(NEIGHBOUR_NODES);
1290  for (GlobalPointersVector<Node>::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
1291  {
1292  if (i->GetValue(IS_VISITED) == 0.0) // the node will be considered for extrapolation only if completely inside
1293  {
1294  const array_1d<double, 3> &inside_press_grad = i->FastGetSolutionStepValue(PRESS_PROJ);
1295  noalias(aux_proj) += inside_press_grad;
1296  avg_number += 1.0;
1297  }
1298  }
1299 
1300  if (avg_number != 0.0) // this case means that it has some neighbours that are completely internal
1301  {
1302  aux_proj /= avg_number;
1303  noalias(iii->FastGetSolutionStepValue(PRESS_PROJ)) = aux_proj;
1304  }
1305  else // case in which there is not a layer of nodes completely internal
1306  {
1307  array_1d<double, 3> body_force = iii->FastGetSolutionStepValue(BODY_FORCE);
1308  array_1d<double, 3> &pproj = iii->FastGetSolutionStepValue(PRESS_PROJ);
1309  for (unsigned int i = 0; i < TDim; i++)
1310  pproj[i] = mRho * body_force[i];
1311  }
1312  }
1313 
1314  // perform extrapolation layer by layer by making an average
1315  // of the neighbours of lower order
1316  for (unsigned int il = 1; il < extrapolation_layers; il++)
1317  {
1318  for (PointIterator iii = layers[il].begin(); iii != layers[il].end(); iii++)
1319  {
1320 
1321  const array_1d<double, 3> &coords_top = iii->Coordinates();
1322 
1323  // extrapolate the average velocity
1324  noalias(aux) = ZeroVector(3);
1325  noalias(aux_proj) = ZeroVector(3);
1326  double avg_number = 0.0;
1327 
1328  double pavg = 0.0;
1329 
1330  GlobalPointersVector<Node> &neighb_nodes = iii->GetValue(NEIGHBOUR_NODES);
1331  for (GlobalPointersVector<Node>::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
1332  {
1333  if (i->GetValue(IS_VISITED) < (il + 1) && i->GetValue(IS_VISITED) != 0.0)
1334  {
1335  const array_1d<double, 3> &coords_bottom = i->Coordinates();
1336  array_1d<double, 3> direction_vec = coords_top;
1337  noalias(direction_vec) -= coords_bottom;
1338  const array_1d<double, 3> &press_grad = i->FastGetSolutionStepValue(PRESS_PROJ);
1339  double temp = inner_prod(direction_vec, press_grad);
1340  double pestimate = i->FastGetSolutionStepValue(PRESSURE, 1) + temp;
1341  pavg += pestimate;
1342  noalias(aux_proj) += press_grad;
1343 
1344  noalias(aux) += i->FastGetSolutionStepValue(VELOCITY);
1345  avg_number += 1.0;
1346  }
1347  }
1348 
1349  if (avg_number != 0.0)
1350  {
1351  aux /= avg_number;
1352  pavg /= avg_number;
1353  aux_proj /= avg_number;
1354  }
1355  else
1356  {
1357  KRATOS_THROW_ERROR(std::runtime_error, "error in extrapolation:: no neighbours find on a extrapolation layer -- impossible", "");
1358  }
1359 
1360  noalias(iii->FastGetSolutionStepValue(VELOCITY)) = aux;
1361  noalias(iii->FastGetSolutionStepValue(VELOCITY, 1)) = aux;
1362 
1363  iii->FastGetSolutionStepValue(PRESSURE, 1) = pavg;
1364 
1365  noalias(iii->FastGetSolutionStepValue(PRESS_PROJ)) = aux_proj;
1366  noalias(iii->FastGetSolutionStepValue(PRESS_PROJ, 1)) = aux_proj;
1367  }
1368  }
1369 
1370  mr_matrix_container.FillVectorFromDatabase(PRESS_PROJ, mXi, mr_model_part.Nodes());
1371 
1372  // mark nodes on which we will have to solve for convection
1373  // mark all of internal nodes
1374  ModelPart::NodesContainerType::iterator it_begin = mr_model_part.NodesBegin();
1375  for (unsigned int i_node = 0; i_node < mr_model_part.Nodes().size(); i_node++)
1376  {
1377  ModelPart::NodesContainerType::iterator it = it_begin + i_node;
1378  if (it->FastGetSolutionStepValue(DISTANCE) <= 0.0)
1379  it->GetValue(IS_VISITED) = 1.0;
1380  else
1381  it->GetValue(IS_VISITED) = 0.0;
1382  }
1383 
1384  // now mark all of the nodes up to the extrapolation layers - 1
1385  for (unsigned int il = 0; il < extrapolation_layers - 1; il++)
1386  for (PointIterator iii = layers[il].begin(); iii != layers[il].end(); iii++)
1387  iii->GetValue(IS_VISITED) = 1.0;
1388 
1389  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, mr_model_part.Nodes());
1390  ApplyVelocityBC(mvel_n1);
1391  mr_matrix_container.WriteVectorToDatabase(VELOCITY, mvel_n1, mr_model_part.Nodes());
1392 
1393  KRATOS_CATCH("")
1394  }
1395 
1396  void ChangeSignToDistance()
1397  {
1398  KRATOS_TRY
1399 
1400  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1401  inode != mr_model_part.NodesEnd();
1402  inode++)
1403  {
1404  double dist = inode->FastGetSolutionStepValue(DISTANCE);
1405  inode->FastGetSolutionStepValue(DISTANCE) = -dist;
1406  }
1407 
1408  KRATOS_CATCH("")
1409  }
1410 
1411  void MarkNodesByDistance(double min, double max)
1412  {
1413  KRATOS_TRY
1414 
1415  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1416  inode != mr_model_part.NodesEnd();
1417  inode++)
1418  {
1419  double dist = inode->FastGetSolutionStepValue(DISTANCE);
1420  if (dist > min && dist < max)
1421  inode->GetValue(IS_VISITED) = 1.0;
1422  else
1423  inode->GetValue(IS_VISITED) = 0.0;
1424  }
1425 
1426  KRATOS_CATCH("")
1427  }
1428 
1429  void SaveScalarVariableToOldStep(Variable<double> &rVar)
1430  {
1431  KRATOS_TRY
1432 
1433  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1434  inode != mr_model_part.NodesEnd();
1435  inode++)
1436  {
1437  inode->FastGetSolutionStepValue(rVar, 1) = inode->FastGetSolutionStepValue(rVar);
1438  }
1439 
1440  KRATOS_CATCH("")
1441  }
1442 
1443  void MarkExternalAndMixedNodes()
1444  {
1445  KRATOS_TRY
1446 
1447  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1448  inode != mr_model_part.NodesEnd();
1449  inode++)
1450  {
1451  inode->GetValue(IS_VISITED) = 0.0;
1452  }
1453 
1454  // detect the nodes inside the fluid surface
1455  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1456  inode != mr_model_part.NodesEnd();
1457  inode++)
1458  {
1459  if (inode->FastGetSolutionStepValue(DISTANCE) > 0.0) // candidates are only the ones inside the fluid domain
1460  {
1461  inode->GetValue(IS_VISITED) = 1.0;
1462  GlobalPointersVector<Node> &neighb_nodes = inode->GetValue(NEIGHBOUR_NODES);
1463  for (GlobalPointersVector<Node>::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
1464  {
1465  i->GetValue(IS_VISITED) = 1.0;
1466  }
1467  }
1468  }
1469  KRATOS_CATCH("")
1470  }
1471 
1472  void MarkInternalAndMixedNodes()
1473  {
1474  KRATOS_TRY
1475 
1476  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1477  inode != mr_model_part.NodesEnd();
1478  inode++)
1479  {
1480  inode->GetValue(IS_VISITED) = 0.0;
1481  }
1482 
1483  // detect the nodes inside the fluid surface
1484  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1485  inode != mr_model_part.NodesEnd();
1486  inode++)
1487  {
1488  if (inode->FastGetSolutionStepValue(DISTANCE) <= 0.0) // candidates are only the ones inside the fluid domain
1489  {
1490  inode->GetValue(IS_VISITED) = 1.0;
1491  GlobalPointersVector<Node> &neighb_nodes = inode->GetValue(NEIGHBOUR_NODES);
1492  for (GlobalPointersVector<Node>::iterator i = neighb_nodes.begin(); i != neighb_nodes.end(); i++)
1493  {
1494  i->GetValue(IS_VISITED) = 1.0;
1495  }
1496  }
1497  }
1498  KRATOS_CATCH("")
1499  }
1500 
1501  void MarkInternalNodes()
1502  {
1503  KRATOS_TRY
1504 
1505  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1506  inode != mr_model_part.NodesEnd();
1507  inode++)
1508  {
1509  inode->GetValue(IS_VISITED) = 0.0;
1510  }
1511 
1512  // detect the nodes inside the fluid surface
1513  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
1514  inode != mr_model_part.NodesEnd();
1515  inode++)
1516  {
1517  if (inode->FastGetSolutionStepValue(DISTANCE) <= 0.0) // candidates are only the ones inside the fluid domain
1518  {
1519  inode->GetValue(IS_VISITED) = 1.0;
1520  }
1521  }
1522  KRATOS_CATCH("")
1523  }
1524 
1525  //**************************************
1526  // function to calculate the area normals
1527 
1528  void CalculateNormals(ModelPart::ConditionsContainerType &rConditions)
1529  {
1530  KRATOS_TRY
1531 
1532  // calculate area normals face-by-face
1533  array_1d<double, 3> area_normal;
1534  // 2D case
1535  if constexpr (TDim == 2)
1536  {
1537  for (ModelPart::ConditionsContainerType::iterator cond_it = rConditions.begin(); cond_it != rConditions.end(); cond_it++)
1538  CalculateNormal2D(cond_it, area_normal);
1539  } // 3D case
1540  else if constexpr (TDim == 3)
1541  {
1542  // help vectors for cross product
1543  array_1d<double, 3> v1;
1544  array_1d<double, 3> v2;
1545  for (ModelPart::ConditionsContainerType::iterator cond_it = rConditions.begin(); cond_it != rConditions.end(); cond_it++)
1546  CalculateNormal3D(cond_it, area_normal, v1, v2);
1547  }
1548 
1549  //(re)initialize normals
1550  unsigned int n_nodes = mNodalFlag.size();
1551  mInOutNormal.resize(n_nodes);
1552  mSlipNormal.resize(n_nodes);
1553 
1554  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
1555  {
1556  noalias(mSlipNormal[i_node]) = ZeroVector(TDim);
1557  mis_slip[i_node] = false;
1558 
1559  noalias(mInOutNormal[i_node]) = ZeroVector(TDim);
1560  }
1561 
1562  // loop over all faces
1563  const double node_factor = 1.0 / TDim;
1564  for (ModelPart::ConditionsContainerType::iterator cond_it = rConditions.begin(); cond_it != rConditions.end(); cond_it++)
1565  {
1566  // get geometry data of the face
1567  Geometry<Node> &face_geometry = cond_it->GetGeometry();
1568 
1569  // reference for area normal of the face
1570  array_1d<double, 3> &face_normal = cond_it->GetValue(NORMAL);
1571 
1572  // slip condition
1573  if (static_cast<bool>(cond_it->Is(SLIP)))
1574  for (unsigned int if_node = 0; if_node < TDim; if_node++)
1575  {
1576  unsigned int i_node = static_cast<unsigned int>(face_geometry[if_node].FastGetSolutionStepValue(AUX_INDEX));
1577  array_1d<double, TDim> &slip_normal = mSlipNormal[i_node];
1578  mis_slip[i_node] = true;
1579  for (unsigned int comp = 0; comp < TDim; comp++)
1580  {
1581  slip_normal[comp] += node_factor * face_normal[comp];
1582  }
1583  }
1584  }
1585 
1586  // fill the list of slip nodes
1587  std::vector<unsigned int> tempmSlipBoundaryList;
1588  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
1589  {
1590  if (mis_slip[i_node] == true)
1591  tempmSlipBoundaryList.push_back(i_node);
1592 
1593  mis_slip[i_node] = false;
1594  }
1595  mSlipBoundaryList.resize(tempmSlipBoundaryList.size(), false);
1596 
1597  IndexPartition<unsigned int>(tempmSlipBoundaryList.size()).for_each([&](unsigned int i){
1598  mSlipBoundaryList[i] = tempmSlipBoundaryList[i];
1599  });
1600 
1601  // loop over all faces to fill inlet outlet
1602 
1603  for (ModelPart::ConditionsContainerType::iterator cond_it = rConditions.begin(); cond_it != rConditions.end(); cond_it++)
1604  {
1605  // get geometry data of the face
1606  Geometry<Node> &face_geometry = cond_it->GetGeometry();
1607 
1608  // reference for area normal of the face
1609  array_1d<double, 3> &face_normal = cond_it->GetValue(NORMAL);
1610 
1611  // inlet or outlet condition
1612  bool is_inlet_or_outlet = false;
1613  if (cond_it->IsNot(SLIP))
1614  is_inlet_or_outlet = true;
1615  else
1616  {
1617  for (unsigned int if_node = 0; if_node < TDim; if_node++)
1618  if (face_geometry[if_node].Is(INLET) || face_geometry[if_node].Is(OUTLET) || face_geometry[if_node].IsFixed(PRESSURE))
1619  is_inlet_or_outlet = true;
1620  }
1621  // slip condition
1622  if (is_inlet_or_outlet) // the opposite of the loop before
1623  for (unsigned int if_node = 0; if_node < TDim; if_node++)
1624  {
1625  unsigned int i_node = static_cast<unsigned int>(face_geometry[if_node].FastGetSolutionStepValue(AUX_INDEX));
1626  array_1d<double, TDim> &inout_normal = mInOutNormal[i_node];
1627  mis_slip[i_node] = true; // reutilize it!
1628  for (unsigned int comp = 0; comp < TDim; comp++)
1629  {
1630  inout_normal[comp] += node_factor * face_normal[comp];
1631  }
1632  }
1633  }
1634 
1635  // fill the list of inlet outlet nodes nodes
1636  std::vector<unsigned int> tempmInOutBoundaryList;
1637  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
1638  {
1639  if (mis_slip[i_node] == true)
1640  tempmInOutBoundaryList.push_back(i_node);
1641  }
1642  mInOutBoundaryList.resize(tempmInOutBoundaryList.size(), false);
1643 
1644  IndexPartition<unsigned int>(tempmInOutBoundaryList.size()).for_each([&](unsigned int i){
1645  mInOutBoundaryList[i] = tempmInOutBoundaryList[i];
1646  });
1647 
1648  KRATOS_CATCH("")
1649  }
1650 
1651  //*******************************
1652  // function to free dynamic memory
1653 
1654  void Clear()
1655  {
1656  KRATOS_TRY
1657  mBodyForce.clear();
1658  mViscosity.clear();
1659  mWork.clear();
1660  mvel_n.clear();
1661  mvel_n1.clear();
1662  mPn.clear();
1663  mPn1.clear();
1664  mHmin.clear();
1665  mHavg.clear();
1666  mSlipNormal.clear();
1667  mNodalFlag.clear();
1668  mFixedVelocities.clear();
1669  mFixedVelocitiesValues.clear();
1670  mPressureOutletList.clear();
1671  mSlipBoundaryList.clear();
1672  mL.clear();
1673  mTauPressure.clear();
1674  mTauConvection.clear();
1675  mTau2.clear();
1676 
1677  mBeta.clear();
1678  mPiConvection.clear();
1679  mphi_n.clear();
1680  mphi_n1.clear();
1681 
1682  mEps.clear();
1683  mA.clear();
1684  mB.clear();
1685  mStrVel.clear();
1686 
1687  mdiv_error.clear();
1688  mdiag_stiffness.clear();
1689  mis_slip.clear();
1690  KRATOS_CATCH("")
1691  }
1692 
1693  void ConvectDistance()
1694  {
1695  KRATOS_TRY
1696 
1697  // variables for node based data handling
1698  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
1699  int n_nodes = rNodes.size();
1700 
1701  // storage of nodal values in local variables
1702  ValuesVectorType rhs, WorkConvection;
1703  rhs.resize(n_nodes);
1704  WorkConvection.resize(n_nodes);
1705 
1706  ValuesVectorType active_nodes;
1707  active_nodes.resize(n_nodes);
1708 
1709  mr_matrix_container.FillScalarFromDatabase(POROSITY, mEps, mr_model_part.Nodes());
1710 
1711  // read variables from Kratos
1712  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, mr_model_part.Nodes());
1713  mr_matrix_container.FillOldVectorFromDatabase(VELOCITY, mvel_n, mr_model_part.Nodes());
1714 
1715  mr_matrix_container.FillScalarFromDatabase(DISTANCE, mphi_n1, mr_model_part.Nodes());
1716  mr_matrix_container.FillOldScalarFromDatabase(DISTANCE, mphi_n, mr_model_part.Nodes());
1717 
1718  // create and fill a vector of nodes for which we want to convect the velocity
1719  for (int i_node = 0; i_node < n_nodes; i_node++)
1720  {
1721  ModelPart::NodesContainerType::iterator it_begin = mr_model_part.NodesBegin();
1722  active_nodes[i_node] = (it_begin + i_node)->GetValue(IS_VISITED);
1723  }
1724 
1725  // calculating the convective projection
1726  array_1d<double, TDim> a_i;
1727  array_1d<double, TDim> a_j;
1728 
1729  #pragma omp parallel for private(a_i, a_j)
1730  for (int i_node = 0; i_node < n_nodes; i_node++)
1731  {
1732  array_1d<double, TDim> &pi_i = mPiConvection[i_node];
1733 
1734  // setting to zero the projection
1735  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
1736  pi_i[l_comp] = 0.0;
1737 
1738  const double &phi_i = mphi_n1[i_node];
1739 
1740  noalias(a_i) = mvel_n1[i_node];
1741  a_i /= mEps[i_node];
1742 
1743  // loop to all the edges surrounding node I
1744  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
1745  {
1746  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
1747 
1748  noalias(a_j) = mvel_n1[j_neighbour];
1749  a_j /= mEps[j_neighbour];
1750 
1751  const double &phi_j = mphi_n1[j_neighbour];
1752 
1753  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
1754 
1755  edge_ij.Add_grad_p(pi_i, phi_i, phi_j);
1756  }
1757 
1758  // apply inverted mass matrix
1759  const double m_inv = mr_matrix_container.GetInvertedMass()[i_node];
1760 
1761  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
1762  pi_i[l_comp] *= m_inv;
1763  }
1764 
1765  // calculating limitor
1766  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
1767 
1768  const array_1d<double, TDim> &pi_i = mPiConvection[i_node];
1769  const double &p_i = mphi_n1[i_node];
1770  double &beta_i = mBeta[i_node];
1771  beta_i = 0.0;
1772  double n = 0.0;
1773  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
1774  {
1775  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
1776 
1777  const double &p_j = mphi_n1[j_neighbour];
1778  const array_1d<double, TDim> &l_k = mEdgeDimensions[csr_index];
1779  const array_1d<double, TDim> &pi_j = mPiConvection[j_neighbour];
1780 
1781  double proj = 0.0;
1782  for (unsigned int comp = 0; comp < TDim; comp++)
1783  proj += 0.5 * l_k[comp] * (pi_i[comp] + pi_j[comp]);
1784 
1785  double numerator = fabs(fabs(p_j - p_i) - fabs(proj));
1786  double denom = fabs(fabs(p_j - p_i) + 1e-6);
1787 
1788  beta_i += numerator / denom;
1789  n += 1.0;
1790  }
1791 
1792  beta_i /= n;
1793  if (beta_i > 1.0)
1794  beta_i = 1.0;
1795  });
1796 
1797  // read time step size from Kratos
1798  ProcessInfo &CurrentProcessInfo = mr_model_part.GetProcessInfo();
1799  double delta_t = CurrentProcessInfo[DELTA_TIME];
1800 
1801  mr_matrix_container.AssignVectorToVector(mphi_n, WorkConvection); // mWork = mphi_n
1802 
1803  // first step of Runge Kutta
1804  mr_matrix_container.SetToZero(rhs);
1805  CalculateRHS_convection(mphi_n1, mvel_n1, rhs, active_nodes);
1806  mr_matrix_container.Add_Minv_value(WorkConvection, WorkConvection, delta_t / 6.0, mr_matrix_container.GetInvertedMass(), rhs);
1807  mr_matrix_container.Add_Minv_value(mphi_n1, mphi_n, 0.5 * delta_t, mr_matrix_container.GetInvertedMass(), rhs);
1808 
1809  // second step
1810  mr_matrix_container.SetToZero(rhs);
1811  CalculateRHS_convection(mphi_n1, mvel_n1, rhs, active_nodes);
1812  mr_matrix_container.Add_Minv_value(WorkConvection, WorkConvection, delta_t / 3.0, mr_matrix_container.GetInvertedMass(), rhs);
1813  mr_matrix_container.Add_Minv_value(mphi_n1, mphi_n, 0.5 * delta_t, mr_matrix_container.GetInvertedMass(), rhs);
1814 
1815  // third step
1816  mr_matrix_container.SetToZero(rhs);
1817  CalculateRHS_convection(mphi_n1, mvel_n1, rhs, active_nodes);
1818  mr_matrix_container.Add_Minv_value(WorkConvection, WorkConvection, delta_t / 3.0, mr_matrix_container.GetInvertedMass(), rhs);
1819  mr_matrix_container.Add_Minv_value(mphi_n1, mphi_n, delta_t, mr_matrix_container.GetInvertedMass(), rhs);
1820 
1821  // fourth step
1822  mr_matrix_container.SetToZero(rhs);
1823  CalculateRHS_convection(mphi_n1, mvel_n1, rhs, active_nodes);
1824  mr_matrix_container.Add_Minv_value(WorkConvection, WorkConvection, delta_t / 6.0, mr_matrix_container.GetInvertedMass(), rhs);
1825 
1826  // compute right-hand side
1827  mr_matrix_container.AssignVectorToVector(WorkConvection, mphi_n1);
1828 
1829  // wetten corner nodes if needed
1830  int corner_size = mcorner_nodes.size();
1831  for (int i = 0; i < corner_size; i++)
1832  {
1833  int i_node = mcorner_nodes[i];
1834  bool to_be_wettened = true;
1835  double min_dist = 0.0;
1836  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
1837  {
1838  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
1839  double neighb_dist = mphi_n1[j_neighbour];
1840  if (min_dist > neighb_dist)
1841  min_dist = neighb_dist;
1842  if (neighb_dist >= 0.0)
1843  {
1844  to_be_wettened = false;
1845  }
1846  }
1847  if (to_be_wettened == true)
1848  mphi_n1[i_node] = min_dist;
1849  }
1850 
1851  mr_matrix_container.WriteScalarToDatabase(DISTANCE, mphi_n1, mr_model_part.Nodes());
1852 
1853  KRATOS_CATCH("")
1854  }
1855 
1856  void ReduceTimeStep(ModelPart &rModelPart, double NewTime)
1857  {
1858  KRATOS_TRY
1859 
1860  rModelPart.OverwriteSolutionStepData(1, 0);
1861  rModelPart.GetProcessInfo().SetCurrentTime(NewTime);
1862 
1863  KRATOS_CATCH("error in reducing the time step")
1864  }
1865 
1866  bool CheckDistanceConvection()
1867  {
1868  int n_large_distance_gradient = 0;
1869  array_1d<double, TDim> grad_d;
1870 
1871  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
1872  int n_nodes = rNodes.size();
1873 
1874  // calculate gradient of distance on the nodes and count occurrences of large gradients (that indicate a failure)
1875  for (int i_node = 0; i_node < n_nodes; i_node++)
1876  {
1877  double dist = mdistances[i_node];
1878 
1879  if (dist <= 0.0)
1880  {
1881  for (unsigned int comp = 0; comp < TDim; comp++)
1882  grad_d[comp] = 0.0;
1883 
1884  double dist_i = mdistances[i_node];
1885 
1886  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
1887  {
1888  // get global index of neighbouring node j
1889  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
1890 
1891  const double &dist_j = mdistances[j_neighbour];
1892 
1893  // projection of pressure gradients
1894  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
1895 
1896  edge_ij.Add_grad_p(grad_d, dist_i, dist_j);
1897  }
1898 
1899  const double &m_inv = mr_matrix_container.GetInvertedMass()[i_node];
1900  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
1901  grad_d[l_comp] *= m_inv;
1902 
1903  double norm_grad = norm_2(grad_d);
1904 
1905  if (norm_grad > 1.5) // large gradient found
1906  n_large_distance_gradient += 1;
1907  }
1908  }
1909 
1910  if (n_large_distance_gradient != 0)
1911  {
1912  bool success = false;
1913  return success;
1914  }
1915  else
1916  {
1917  bool success = true;
1918  return success;
1919  }
1920  }
1921 
1922  void ActivateWallResistance(double Ywall)
1923  {
1924  mWallLawIsActive = true;
1925  mY_wall = Ywall;
1926  }
1927 
1928  double ComputeVolumeVariation()
1929  {
1930  ProcessInfo &CurrentProcessInfo = mr_model_part.GetProcessInfo();
1931  double dt = CurrentProcessInfo[DELTA_TIME];
1932 
1933  // slip condition
1934  int inout_size = mInOutBoundaryList.size();
1935  double vol_var = 0.0;
1936 
1937  for (int i = 0; i < inout_size; i++)
1938  {
1939  unsigned int i_node = mInOutBoundaryList[i];
1940  double dist = mdistances[i_node];
1941  if (dist <= 0.0)
1942  {
1943  const array_1d<double, TDim> &U_i = mvel_n1[i_node];
1944  const array_1d<double, TDim> &an_i = mInOutNormal[i_node];
1945  double projection_length = 0.0;
1946  for (unsigned int comp = 0; comp < TDim; comp++)
1947  {
1948  projection_length += U_i[comp] * an_i[comp];
1949  }
1950  vol_var += projection_length;
1951  }
1952  }
1953  return vol_var * dt;
1954  }
1955 
1956  double ComputeWetVolume()
1957  {
1958  KRATOS_TRY
1959 
1960  mr_matrix_container.FillScalarFromDatabase(DISTANCE, mdistances, mr_model_part.Nodes());
1961 
1962  // slip condition
1963 
1964  double wet_volume = 0.0;
1965 
1966  for (int i = 0; i < static_cast<int>(mdistances.size()); i++)
1967  {
1968  double dist = mdistances[i];
1969  const double m_inv = mr_matrix_container.GetInvertedMass()[i];
1970 
1971  if (dist <= 0.0)
1972  {
1973  wet_volume += 1.0 / m_inv;
1974  }
1975  }
1976 
1977  return wet_volume;
1978 
1979  KRATOS_CATCH("");
1980  }
1981 
1982  void DiscreteVolumeCorrection(double expected_volume, double measured_volume)
1983  {
1984  double volume_error = expected_volume - measured_volume;
1985  if (measured_volume < expected_volume)
1986  {
1987  double layer_volume = 0.0;
1988  std::vector<unsigned int> first_outside;
1989  int n_nodes = mdistances.size();
1990 
1991  // find list of the first nodes outside of the fluid and compute their volume
1992  for (int i_node = 0; i_node < n_nodes; i_node++)
1993  {
1994  double dist = mdistances[i_node];
1995  if (dist > 0.0) // node is outside domain
1996  {
1997  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
1998  {
1999  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2000  if (mdistances[j_neighbour] <= 0.0)
2001  {
2002  const double nodal_mass = 1.0 / mr_matrix_container.GetInvertedMass()[i_node];
2003 
2004  if (nodal_mass < volume_error - layer_volume)
2005  {
2006  first_outside.push_back(i_node);
2007  layer_volume += nodal_mass;
2008  }
2009  }
2010  }
2011  }
2012  // mark the nodes in the outside layer with a small negative distance
2013  for (unsigned int i = 0; i < first_outside.size(); i++)
2014  {
2015  unsigned int i_node = first_outside[i];
2016  mdistances[i_node] = -mHavg[i_node];
2017  }
2018  }
2019  }
2020 
2021  mr_matrix_container.WriteScalarToDatabase(DISTANCE, mdistances, mr_model_part.Nodes());
2022  }
2023 
2024  void PushFreeSurface()
2025  {
2026 
2027  // double layer_volume = 0.0;
2028  std::vector<unsigned int> first_outside;
2029  int n_nodes = mdistances.size();
2030 
2031  // find list of the first nodes outside of the fluid and compute their volume
2032  for (int i_node = 0; i_node < n_nodes; i_node++)
2033  {
2034  double dist = mdistances[i_node];
2035  if (dist > 0.0) // node is outside domain
2036  {
2037  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
2038  {
2039  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2040  if (mdistances[j_neighbour] <= 0.0)
2041  {
2042  // mark the nodes in the outside layer with a small negative distance
2043  mdistances[i_node] = -mHavg[i_node];
2044  }
2045  }
2046  }
2047  }
2048 
2049  mr_matrix_container.WriteScalarToDatabase(DISTANCE, mdistances, mr_model_part.Nodes());
2050  }
2051 
2052  //***************************************
2053  // function to set adequate time step size
2054 
2055  double ComputeBoundedTimeStep(const double CFLNumber, const double MaxDt)
2056  {
2057  KRATOS_TRY
2058 
2059  // save the maximum time step
2060  max_dt = MaxDt;
2061 
2062  // local variable for time step size
2063  double delta_t = 1e10; // max_dt;
2064 
2065  mdelta_t_avg = 1e10; // max_dt;
2066 
2067  // getting value of current velocity and of viscosity
2068  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, mr_model_part.Nodes());
2069  mr_matrix_container.FillVectorFromDatabase(BODY_FORCE, mBodyForce, mr_model_part.Nodes());
2070  mr_matrix_container.FillScalarFromDatabase(VISCOSITY, mViscosity, mr_model_part.Nodes());
2071 
2072  mr_matrix_container.FillScalarFromDatabase(POROSITY, mEps, mr_model_part.Nodes());
2073  mr_matrix_container.FillScalarFromDatabase(LIN_DARCY_COEF, mA, mr_model_part.Nodes());
2074  mr_matrix_container.FillScalarFromDatabase(NONLIN_DARCY_COEF, mB, mr_model_part.Nodes());
2075  mr_matrix_container.FillVectorFromDatabase(STRUCTURE_VELOCITY, mStrVel, mr_model_part.Nodes());
2076 
2077  // loop over all nodes
2078  double n_nodes = mvel_n1.size();
2079  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
2080  {
2081  array_1d<double, TDim> &v_i = mvel_n1[i_node];
2082  const double havg_i = mHavg[i_node];
2083  const double hmin_i = mHmin[i_node];
2084  const double eps_i = mEps[i_node];
2085  const double nu_i = mViscosity[i_node];
2086 
2087  double vel_norm = 0.0;
2088  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
2089  {
2090  vel_norm += v_i[l_comp] * v_i[l_comp];
2091  }
2092  vel_norm = sqrt(vel_norm);
2093 
2094  vel_norm /= eps_i;
2095 
2096  // use CFL condition to compute time step size
2097  double delta_t_i = CFLNumber * 1.0 / (2.0 * vel_norm / hmin_i + 4.0 * nu_i / (hmin_i * hmin_i));
2098  double delta_t_i_avg = 1.0 / (2.0 * vel_norm / havg_i + 4.0 * nu_i / (havg_i * havg_i));
2099 
2100  if (delta_t_i < 10e-8) // NO PHYSICS AT ALL!!!!! bounding the delata_t to 10e-08 by reducing the velocity!!
2101  {
2102  v_i *= delta_t_i / 10e-8;
2103  delta_t_i = 10e-8;
2104  }
2105 
2106  if (delta_t_i_avg < 10e-8) // NO PHYSICS AT ALL!!!!! bounding the delta_t_i_avg to 10e-08 by reducing the velocity!!
2107  {
2108  v_i *= delta_t_i_avg / 10e-8;
2109  delta_t_i_avg = 10e-8;
2110  }
2111 
2112  // considering the most restrictive case of neighbor's velocities with similar direction but opposite sense.
2113  // loop over all neighbours
2114  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
2115  {
2116  // get global index of neighbouring node j
2117  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2118 
2119  array_1d<double, TDim> &v_j = mvel_n1[j_neighbour];
2120 
2121  double v_diff_norm = 0.0;
2122  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
2123  {
2124  double temp = v_i[l_comp] - v_j[l_comp];
2125  v_diff_norm += temp * temp;
2126  }
2127  v_diff_norm = sqrt(v_diff_norm);
2128  v_diff_norm /= eps_i;
2129  double delta_t_j = CFLNumber * 1.0 / (2.0 * v_diff_norm / hmin_i + 4.0 * nu_i / (hmin_i * hmin_i));
2130 
2131  if (delta_t_j < 10e-8) // NO PHYSICS AT ALL!!!!! bounding the delata_t to 10e-08 by reducing the velocity!!
2132  {
2133  v_j *= delta_t_j / 10e-8;
2134  delta_t_j = 10e-8;
2135  }
2136 
2137  if (delta_t_j < delta_t_i)
2138  delta_t_i = delta_t_j;
2139  }
2140 
2141  // choose the overall minimum of delta_t_i
2142  if (delta_t_i < delta_t)
2143  delta_t = delta_t_i;
2144 
2145  if (delta_t_i_avg < mdelta_t_avg)
2146  mdelta_t_avg = delta_t_i_avg;
2147  }
2148  //*******************
2149  // perform MPI syncronization of the dt (minimum should be kept)
2150 
2151  if (delta_t <= 10 - 7) // writing back the changed velocities
2152  mr_matrix_container.WriteVectorToDatabase(VELOCITY, mvel_n1, mr_model_part.Nodes());
2153 
2154  return delta_t;
2155 
2156  KRATOS_CATCH("")
2157  }
2158 
2159  void CalculatePorousResistanceLaw(unsigned int res_law)
2160  {
2161  if (res_law == 1)
2162  {
2163  /* if the chosen resistance law is ERGUN calculate Ergun A and B*/
2164  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
2165  inode != mr_model_part.NodesEnd();
2166  inode++)
2167  {
2168  const double eps = inode->FastGetSolutionStepValue(POROSITY); /*reading from kratos database*/
2169  const double d = inode->FastGetSolutionStepValue(DIAMETER); /*reading from kratos database*/
2170  const double nu = inode->FastGetSolutionStepValue(VISCOSITY); /*reading from kratos database*/
2171  double &a = inode->FastGetSolutionStepValue(LIN_DARCY_COEF); /*changing kratos database*/
2172  double &b = inode->FastGetSolutionStepValue(NONLIN_DARCY_COEF); /*changing kratos database*/
2173  if (eps < 1.0)
2174  {
2175  double k_inv = 150.0 * (1.0 - eps) * (1.0 - eps) / (eps * eps * eps * d * d);
2176  a = nu * k_inv;
2177  b = (1.75 / eps) * sqrt(k_inv / (150.0 * eps));
2178  }
2179  else
2180  {
2181  a = 0.0;
2182  b = 0.0;
2183  }
2184  }
2185  }
2186  else
2187  {
2188  /* whether it is a Custom Resistance law or NO resistance law is present ---> set to zero A and B for non porous nodes*/
2189  for (ModelPart::NodesContainerType::iterator inode = mr_model_part.NodesBegin();
2190  inode != mr_model_part.NodesEnd();
2191  inode++)
2192  {
2193  const double eps = inode->FastGetSolutionStepValue(POROSITY); /*reading from kratos database*/
2194 
2195  double &a = inode->FastGetSolutionStepValue(LIN_DARCY_COEF); /*changing kratos database*/
2196  double &b = inode->FastGetSolutionStepValue(NONLIN_DARCY_COEF); /*changing kratos database*/
2197  if (eps == 1.0)
2198  {
2199  a = 0.0;
2200  b = 0.0;
2201  }
2202  }
2203  }
2204 
2205  mr_matrix_container.FillScalarFromDatabase(LIN_DARCY_COEF, mA, mr_model_part.Nodes()); /*filling edgebased database reading from kratos database*/
2206  mr_matrix_container.FillScalarFromDatabase(NONLIN_DARCY_COEF, mB, mr_model_part.Nodes()); /*filling edgebased database reading from kratos database*/
2207  }
2208 
2209  private:
2210  MatrixContainer &mr_matrix_container;
2211  ModelPart &mr_model_part;
2212 
2213  bool muse_mass_correction;
2214 
2215  // parameters controlling the wall law
2216  bool mWallLawIsActive;
2217  double mY_wall;
2218 
2219  // parameters for controlling the usage of the delta time in the stabilization
2220  double mstabdt_pressure_factor;
2221  double mstabdt_convection_factor;
2222  double mtau2_factor;
2223  bool massume_constant_dp;
2224 
2225  // nodal values
2226  ValuesVectorType mViscosity;
2227  CalcVectorType mBodyForce;
2228  // velocity vector U at time steps n and n+1
2229  CalcVectorType mWork, mvel_n, mvel_n1, mx;
2230  // pressure vector p at time steps n and n+1
2231  ValuesVectorType mPn, mPn1;
2232  // coefficients
2233  ValuesVectorType mdistances;
2234  // minimum length of the edges surrounding edges surrounding each nodal point
2235  ValuesVectorType mHmin;
2236  ValuesVectorType mHavg;
2237  CalcVectorType mEdgeDimensions;
2238 
2239  // area normal
2240  CalcVectorType mSlipNormal;
2241  CalcVectorType mInOutNormal;
2242 
2243  // projection terms
2244  CalcVectorType mPi, mXi;
2245 
2246  // flag for first time step
2247  bool mFirstStep;
2248 
2249  // flag to differentiate interior and boundary nodes
2250  ValuesVectorType mNodalFlag;
2251  // lists of nodes with different types of boundary conditions
2252  IndicesVectorType mSlipBoundaryList, mPressureOutletList, mFixedVelocities, mInOutBoundaryList;
2253  CalcVectorType mFixedVelocitiesValues;
2254  // ValuesVectorType mPressureOutlet;
2255 
2256  // intrinsic time step size
2257  ValuesVectorType mTauPressure;
2258  ValuesVectorType mTauConvection;
2259  ValuesVectorType mTau2;
2260 
2261  ValuesVectorType mdiv_error;
2262  std::vector<bool> mis_slip;
2263  // variables for resolving pressure equation
2264  // laplacian matrix
2265  TSystemMatrixType mL;
2266 
2267  // constant variables
2268  double mRho;
2269 
2270  // variables for convection
2271  ValuesVectorType mphi_n;
2272  ValuesVectorType mphi_n1;
2273  CalcVectorType mPiConvection;
2274  ValuesVectorType mBeta;
2275 
2276  // variables for edge BCs
2277  IndicesVectorType medge_nodes;
2278  CalcVectorType medge_nodes_direction;
2279  IndicesVectorType mcorner_nodes;
2280 
2281  ValuesVectorType mEps;
2282  ValuesVectorType mdiag_stiffness;
2283  ValuesVectorType mA;
2284  ValuesVectorType mB;
2285  CalcVectorType mStrVel;
2286 
2287  double mdelta_t_avg;
2288  double max_dt;
2289 
2290  double mshock_coeff;
2291 
2292  //***********************************************************
2293  // functions to calculate area normals for boundary conditions
2294 
2295  void CalculateNormal2D(ModelPart::ConditionsContainerType::iterator cond_it, array_1d<double, 3> &area_normal)
2296  {
2297  Geometry<Node> &face_geometry = (cond_it)->GetGeometry();
2298 
2299  area_normal[0] = face_geometry[1].Y() - face_geometry[0].Y();
2300  area_normal[1] = -(face_geometry[1].X() - face_geometry[0].X());
2301  area_normal[2] = 0.00;
2302 
2303  noalias((cond_it)->GetValue(NORMAL)) = area_normal;
2304  }
2305 
2306  void CalculateNormal3D(ModelPart::ConditionsContainerType::iterator cond_it, array_1d<double, 3> &area_normal, array_1d<double, 3> &v1, array_1d<double, 3> &v2)
2307  {
2308  Geometry<Node> &face_geometry = (cond_it)->GetGeometry();
2309 
2310  v1[0] = face_geometry[1].X() - face_geometry[0].X();
2311  v1[1] = face_geometry[1].Y() - face_geometry[0].Y();
2312  v1[2] = face_geometry[1].Z() - face_geometry[0].Z();
2313 
2314  v2[0] = face_geometry[2].X() - face_geometry[0].X();
2315  v2[1] = face_geometry[2].Y() - face_geometry[0].Y();
2316  v2[2] = face_geometry[2].Z() - face_geometry[0].Z();
2317 
2318  MathUtils<double>::CrossProduct(area_normal, v1, v2);
2319  area_normal *= -0.5;
2320 
2321  noalias((cond_it)->GetValue(NORMAL)) = area_normal;
2322  }
2323 
2324  //*********************************************************
2325  // function to calculate minimum length of surrounding edges
2326 
2327  void CalculateEdgeLengths(ModelPart::NodesContainerType &rNodes)
2328  {
2329  KRATOS_TRY
2330 
2331  // get number of nodes
2332  unsigned int n_nodes = rNodes.size();
2333  // reserve memory for storage of nodal coordinates
2334  std::vector<array_1d<double, 3>> position;
2335  position.resize(n_nodes);
2336 
2337  // get position of all nodes
2338  for (typename ModelPart::NodesContainerType::iterator node_it = rNodes.begin(); node_it != rNodes.end(); node_it++)
2339  {
2340  // get the global index of the node
2341  unsigned int i_node = static_cast<unsigned int>(node_it->FastGetSolutionStepValue(AUX_INDEX));
2342  // save its coordinates locally
2343  noalias(position[i_node]) = node_it->Coordinates();
2344 
2345  // initialize minimum edge length with relatively big values
2346  mHmin[i_node] = 1e10;
2347  }
2348 
2349  ValuesVectorType &aaa = mr_matrix_container.GetHmin();
2350  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
2351  {
2352  mHmin[i_node] = aaa[i_node];
2353  }
2354 
2355  // take unstructured meshes into account
2356  if constexpr (TDim == 2)
2357  {
2358  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
2359  {
2360  double &h_i = mHavg[i_node];
2361  double &m_i = mr_matrix_container.GetLumpedMass()[i_node];
2362 
2363  h_i = sqrt(2.0 * m_i);
2364  }
2365  }
2366  else if constexpr (TDim == 3)
2367  {
2368  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
2369  {
2370  double &h_i = mHavg[i_node];
2371  double &m_i = mr_matrix_container.GetLumpedMass()[i_node];
2372 
2373  h_i = pow(6.0 * m_i, 1.0 / 3.0);
2374  }
2375  }
2376 
2377  // compute edge coordinates
2378  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
2379  {
2380  array_1d<double, 3> &pos_i = position[i_node];
2381 
2382  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
2383  {
2384  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2385  array_1d<double, 3> &pos_j = position[j_neighbour];
2386 
2387  array_1d<double, TDim> &l_k = mEdgeDimensions[csr_index];
2388  for (unsigned int comp = 0; comp < TDim; comp++)
2389  l_k[comp] = pos_i[comp] - pos_j[comp];
2390  }
2391  }
2392 
2393  KRATOS_CATCH("")
2394  }
2395 
2396  //*********************************************************************
2397  // function to calculate right-hand side of fractional momentum equation
2398 
2399  void CalculateRHS_convection(
2400  const ValuesVectorType &mphi,
2401  const CalcVectorType &convective_velocity,
2402  ValuesVectorType &rhs,
2403  ValuesVectorType &active_nodes)
2404  {
2405  KRATOS_TRY
2406 
2407  int n_nodes = mphi.size();
2408 
2409  // calculating the RHS
2410  double stab_low;
2411  double stab_high;
2412  array_1d<double, TDim> a_i;
2413  array_1d<double, TDim> a_j;
2414 
2415  #pragma omp parallel for private(stab_low, stab_high, a_i, a_j)
2416  for (int i_node = 0; i_node < n_nodes; i_node++)
2417  {
2418  double &rhs_i = rhs[i_node];
2419  const double &h_i = mHavg[i_node];
2420  const double &phi_i = mphi[i_node];
2421  noalias(a_i) = convective_velocity[i_node];
2422  a_i /= mEps[i_node];
2423 
2424  const array_1d<double, TDim> &proj_i = mPiConvection[i_node];
2425  double pi_i = proj_i[0] * a_i[0];
2426  for (unsigned int l_comp = 1; l_comp < TDim; l_comp++)
2427  pi_i += proj_i[l_comp] * a_i[l_comp];
2428 
2429  rhs_i = 0.0;
2430 
2431  if (active_nodes[i_node] != 0.0)
2432  {
2433  const double &beta = mBeta[i_node];
2434 
2435  double norm_a = a_i[0] * a_i[0];
2436  for (unsigned int l_comp = 1; l_comp < TDim; l_comp++)
2437  norm_a += a_i[l_comp] * a_i[l_comp];
2438  norm_a = sqrt(norm_a);
2439 
2440  // loop to all the edges surrounding node I
2441  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
2442  {
2443  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2444 
2445  if (active_nodes[j_neighbour] != 0.0)
2446  {
2447  const double &phi_j = mphi[j_neighbour];
2448  noalias(a_j) = convective_velocity[j_neighbour];
2449  a_j /= mEps[j_neighbour];
2450 
2451  const array_1d<double, TDim> &proj_j = mPiConvection[j_neighbour];
2452  double pi_j = proj_j[0] * a_i[0];
2453  for (unsigned int l_comp = 1; l_comp < TDim; l_comp++)
2454  pi_j += proj_j[l_comp] * a_i[l_comp];
2455 
2456  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
2457 
2458  // convection operator
2459  edge_ij.Sub_ConvectiveContribution(rhs_i, a_i, phi_i, a_j, phi_j);
2460 
2461  // calculate stabilization part
2462  edge_ij.CalculateConvectionStabilization_LOW(stab_low, a_i, phi_i, a_j, phi_j);
2463 
2464  double edge_tau = mTauConvection[i_node];
2465 
2466  edge_ij.CalculateConvectionStabilization_HIGH(stab_high, a_i, pi_i, a_j, pi_j);
2467 
2468  edge_ij.Sub_StabContribution(rhs_i, edge_tau, 1.0, stab_low, stab_high);
2469 
2470  double coeff = 0.5 * mshock_coeff; //=0.7*0.5;
2471  double laplacian_ij = 0.0;
2472  edge_ij.CalculateScalarLaplacian(laplacian_ij);
2473  double capturing = laplacian_ij * (phi_j - phi_i);
2474 
2475  double aaa = 0.0;
2476  for (unsigned int k_comp = 0; k_comp < TDim; k_comp++)
2477  for (unsigned int m_comp = 0; m_comp < TDim; m_comp++)
2478  aaa += a_i[k_comp] * a_i[m_comp] * edge_ij.LaplacianIJ(k_comp, m_comp);
2479 
2480  if (norm_a > 1e-10)
2481  {
2482  aaa /= (norm_a * norm_a);
2483  double capturing2 = aaa * (phi_j - phi_i);
2484 
2485  if (fabs(capturing) > fabs(capturing2))
2486  rhs_i -= coeff * (capturing - capturing2) * beta * norm_a * h_i;
2487  }
2488  }
2489  }
2490  }
2491  }
2492 
2493  KRATOS_CATCH("")
2494  }
2495 
2496  //**************************************
2497 
2498  void CornerDectectionHelper(Geometry<Node> &face_geometry,
2499  const array_1d<double, 3> &face_normal,
2500  const double An,
2501  const GlobalPointersVector<Condition> &neighb,
2502  const unsigned int i1,
2503  const unsigned int i2,
2504  const unsigned int neighb_index,
2505  std::vector<unsigned int> &edge_nodes,
2506  CalcVectorType &cornern_list)
2507  {
2508  double acceptable_angle = 45.0 / 180.0 * 3.1; // angles of less than 45 deg will be accepted
2509  double acceptable_cos = cos(acceptable_angle);
2510 
2511  if (face_geometry[i1].Id() < face_geometry[i2].Id()) // we do this to add the face ones
2512  {
2513  const array_1d<double, 3> &neighb_normal = neighb[neighb_index].GetValue(NORMAL);
2514  double neighb_An = norm_2(neighb_normal);
2515 
2516  double cos_normal = 1.0 / (An * neighb_An) * inner_prod(face_normal, neighb_normal);
2517 
2518  // if the angle is too big between the two normals then the edge in the middle is a corner
2519  if (cos_normal < acceptable_cos)
2520  {
2521  array_1d<double, 3> edge = face_geometry[i2].Coordinates() - face_geometry[i1].Coordinates();
2522  double temp = norm_2(edge);
2523  edge /= temp;
2524 
2525  int index1 = face_geometry[i1].FastGetSolutionStepValue(AUX_INDEX);
2526  int index2 = face_geometry[i2].FastGetSolutionStepValue(AUX_INDEX);
2527 
2528  edge_nodes[index1] += 1;
2529  edge_nodes[index2] += 1;
2530 
2531  double sign1 = 0.0;
2532  for (unsigned int i = 0; i < edge.size(); i++)
2533  {
2534  sign1 += cornern_list[index1][i] * edge[i];
2535  }
2536 
2537  if (sign1 >= 0)
2538  {
2539  for (unsigned int i = 0; i < edge.size(); i++)
2540  cornern_list[index1][i] += edge[i];
2541  }
2542  else
2543  {
2544  for (unsigned int i = 0; i < edge.size(); i++)
2545  cornern_list[index1][i] -= edge[i];
2546  }
2547 
2548  double sign2 = inner_prod(cornern_list[index2], edge);
2549  if (sign2 >= 0)
2550  {
2551  for (unsigned int i = 0; i < edge.size(); i++)
2552  cornern_list[index2][i] += edge[i];
2553  }
2554  else
2555  {
2556  for (unsigned int i = 0; i < edge.size(); i++)
2557  cornern_list[index2][i] -= edge[i];
2558  }
2559  }
2560  }
2561  }
2562 
2563  // function to calculate the area normals
2564 
2565  void DetectEdges3D(ModelPart::ConditionsContainerType &rConditions)
2566  {
2567  KRATOS_TRY
2568 
2569  // calculate area normals face-by-face
2570  array_1d<double, 3> area_normal;
2571 
2572  //(re)initialize normals
2573  unsigned int n_nodes = mNodalFlag.size();
2574  std::vector<unsigned int> temp_edge_nodes(n_nodes);
2575  CalcVectorType temp_cornern_list(n_nodes);
2576  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
2577  {
2578  temp_edge_nodes[i_node] = 0.0;
2579  noalias(temp_cornern_list[i_node]) = ZeroVector(TDim);
2580  }
2581 
2582  // loop over all faces
2583  for (ModelPart::ConditionsContainerType::iterator cond_it = rConditions.begin(); cond_it != rConditions.end(); cond_it++)
2584  {
2585  // get geometry data of the face
2586  Geometry<Node> &face_geometry = cond_it->GetGeometry();
2587 
2588  // reference for area normal of the face
2589  const array_1d<double, 3> &face_normal = cond_it->GetValue(NORMAL);
2590  double An = norm_2(face_normal);
2591 
2592  unsigned int current_id = cond_it->Id();
2593 
2594  // slip condition
2595  if (cond_it->Is(SLIP)) // this is a slip face --> now look for its neighbours
2596  {
2597  const GlobalPointersVector<Condition> &neighb = cond_it->GetValue(NEIGHBOUR_CONDITIONS);
2598 
2599  // check for neighbour zero
2600  if (neighb[0].Id() != current_id) // check if the neighbour exists
2601  CornerDectectionHelper(face_geometry, face_normal, An, neighb, 1, 2, 0, temp_edge_nodes, temp_cornern_list);
2602 
2603  // check for neighbour one
2604  if (neighb[1].Id() != current_id) // check if the neighbour exists
2605  CornerDectectionHelper(face_geometry, face_normal, An, neighb, 2, 0, 1, temp_edge_nodes, temp_cornern_list);
2606 
2607  // check for neighbour two
2608  if (neighb[2].Id() != current_id) // check if the neighbour exists
2609  CornerDectectionHelper(face_geometry, face_normal, An, neighb, 0, 1, 2, temp_edge_nodes, temp_cornern_list);
2610  }
2611  }
2612 
2613  // fill the list of edge_nodes
2614  std::vector<unsigned int> tempmedge_nodes;
2615  std::vector<array_1d<double, TDim>> tempmedge_nodes_direction;
2616  std::vector<unsigned int> tempmcorner_nodes;
2617  for (unsigned int i_node = 0; i_node < n_nodes; i_node++)
2618  {
2619  if (temp_edge_nodes[i_node] == 2) // node is a edge_node
2620  {
2621  tempmedge_nodes.push_back(i_node);
2622  array_1d<double, TDim> &node_edge = temp_cornern_list[i_node];
2623 
2624  node_edge /= norm_2(node_edge);
2625  tempmedge_nodes_direction.push_back(node_edge);
2626  }
2627  else if (temp_edge_nodes[i_node] > 2)
2628  tempmcorner_nodes.push_back(i_node);
2629  }
2630  medge_nodes.resize(tempmedge_nodes.size(), false);
2631  medge_nodes_direction.resize(tempmedge_nodes_direction.size(), false);
2632  mcorner_nodes.resize(tempmcorner_nodes.size(), false);
2633 
2634  IndexPartition<unsigned int>(tempmedge_nodes.size()).for_each([&](unsigned int i){
2635  medge_nodes[i] = tempmedge_nodes[i];
2636  medge_nodes_direction[i] = tempmedge_nodes_direction[i];
2637  });
2638 
2639  IndexPartition<unsigned int>(tempmcorner_nodes.size()).for_each([&](unsigned int i){
2640  mcorner_nodes[i] = tempmcorner_nodes[i];
2641  });
2642 
2643  for (int i = 0; i < static_cast<int>(mcorner_nodes.size()); i++)
2644  {
2645  KRATOS_WATCH(mcorner_nodes[i]);
2646  }
2647 
2648  KRATOS_CATCH("")
2649  }
2650 
2651  double ComputePorosityCoefficient(const double &vel_norm, const double &eps, const double &a, const double &b)
2652  {
2653  double linear;
2654  double non_linear;
2655 
2656  linear = eps * a;
2657  non_linear = eps * b * vel_norm;
2658 
2659  return linear + non_linear;
2660  }
2661 
2662  void LaplacianSmooth(ValuesVectorType &to_be_smoothed, ValuesVectorType &aux)
2663  {
2664  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
2665  int n_nodes = rNodes.size();
2666 
2667  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
2668  double dist = mdistances[i_node];
2669 
2670  double correction = 0.0;
2671  const double &origin_i = to_be_smoothed[i_node];
2672 
2673  if (dist <= 0.0) // node is inside domain ---- if outside do nothing
2674  {
2675 
2676  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
2677  {
2678  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2679  const double &origin_j = to_be_smoothed[j_neighbour];
2680 
2681  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
2682 
2683  double l_ikjk;
2684  edge_ij.CalculateScalarLaplacian(l_ikjk);
2685 
2686  correction += l_ikjk * (origin_j - origin_i);
2687  }
2688  }
2689 
2690  aux[i_node] = origin_i - correction;
2691  });
2692 
2693  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
2694  to_be_smoothed[i_node] = aux[i_node];
2695  });
2696  }
2697 
2698  void ComputeWallResistance(
2699  const CalcVectorType &vel,
2700  ValuesVectorType &diag_stiffness
2701  )
2702  {
2703  // parameters:
2704  double k = 0.41;
2705  double B = 5.1;
2706  double toll = 1e-6;
2707  double ym = mY_wall;
2708  double y_plus_incercept = 10.9931899;
2709  unsigned int itmax = 100;
2710 
2711  // slip condition
2712  int slip_size = mSlipBoundaryList.size();
2713 
2714  #pragma omp parallel for firstprivate(slip_size, B, toll, ym, y_plus_incercept, itmax)
2715  for (int i_slip = 0; i_slip < slip_size; i_slip++)
2716  {
2717  unsigned int i_node = mSlipBoundaryList[i_slip];
2718  double dist = mdistances[i_node];
2719  if (dist <= 0.0)
2720  {
2721 
2722  KRATOS_ERROR_IF(mViscosity[i_node] == 0) << "it is not possible to use the wall law with 0 viscosity" << std::endl;
2723 
2724  const double nu = mViscosity[i_node];
2725 
2726  // array_1d<double, TDim>& rhs_i = rhs[i_node];
2727  const array_1d<double, TDim> &U_i = vel[i_node];
2728  const array_1d<double, TDim> &an_i = mSlipNormal[i_node];
2729 
2730  // compute the modulus of the velocity
2731  double mod_vel = 0.0;
2732  double area = 0.0;
2733  for (unsigned int comp = 0; comp < TDim; comp++)
2734  {
2735  mod_vel += U_i[comp] * U_i[comp];
2736  area += an_i[comp] * an_i[comp];
2737  }
2738  mod_vel = sqrt(mod_vel);
2739  area = sqrt(area);
2740  diag_stiffness[i_node] += area * mod_vel / pow(1.0 / k * log(100.00) + B, 2); /* * mWallReductionFactor[ i_node ];*/
2741  // now compute the skin friction
2742  double mod_uthaw = sqrt(mod_vel * nu / ym);
2743  const double y_plus = ym * mod_uthaw / nu;
2744 
2745  if (y_plus > y_plus_incercept)
2746  {
2747  // begin cicle to calculate the real u_thaw's module:
2748  unsigned int it = 0;
2749  double dx = 1e10;
2750 
2751  while (fabs(dx) > toll * mod_uthaw && it < itmax)
2752  {
2753  double a = 1.0 / k;
2754  double temp = a * log(ym * mod_uthaw / nu) + B;
2755  double y = mod_uthaw * (temp)-mod_vel;
2756  double y1 = temp + a;
2757  dx = y / y1;
2758  mod_uthaw -= dx;
2759  it = it + 1;
2760  }
2761 
2762  if (it == itmax)
2763  std::cout << "attention max number of iterations exceeded in wall law computation" << std::endl;
2764  }
2765  }
2766  else
2767  diag_stiffness[i_node] += 0.0;
2768  }
2769  }
2770 
2771  void ApplySmagorinsky3D(double MolecularViscosity, double Cs)
2772  {
2773  KRATOS_TRY
2774  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
2775  // calculating the RHS
2776  array_1d<double, TDim> grad_vx;
2777  array_1d<double, TDim> grad_vy;
2778  array_1d<double, TDim> grad_vz;
2779  int n_nodes = rNodes.size();
2780  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes);
2781  array_1d<double, TDim> stab_high;
2782 
2783  #pragma omp parallel for private(grad_vx, grad_vy, grad_vz)
2784  for (int i_node = 0; i_node < n_nodes; i_node++)
2785  {
2786  // set to zero the gradients
2787  for (unsigned int comp = 0; comp < TDim; comp++)
2788  {
2789  grad_vx[comp] = 0.0;
2790  grad_vy[comp] = 0.0;
2791  grad_vz[comp] = 0.0;
2792  }
2793  // compute node by node the gradients
2794  const array_1d<double, TDim> &U_i = mvel_n1[i_node];
2795  const double h = mHmin[i_node];
2796  const double m_inv = mr_matrix_container.GetInvertedMass()[i_node];
2797  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
2798  {
2799  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2800  const array_1d<double, TDim> &U_j = mvel_n1[j_neighbour];
2801  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
2802  edge_ij.Add_grad_p(grad_vx, U_i[0], U_j[0]);
2803  edge_ij.Add_grad_p(grad_vy, U_i[1], U_j[1]);
2804  edge_ij.Add_grad_p(grad_vz, U_i[2], U_j[2]);
2805  }
2806  // finalize computation of the gradients
2807  // set to zero the gradients
2808  for (unsigned int comp = 0; comp < TDim; comp++)
2809  {
2810  grad_vx[comp] *= m_inv;
2811  grad_vy[comp] *= m_inv;
2812  grad_vz[comp] *= m_inv;
2813  }
2814  // symmetrize and multiply by 2
2815  grad_vx[0] *= 2.0;
2816  grad_vy[1] *= 2.0;
2817  grad_vz[2] *= 2.0;
2818  grad_vx[1] += grad_vy[0];
2819  grad_vx[2] += grad_vz[0];
2820  grad_vy[2] += grad_vz[1];
2821  grad_vy[0] += grad_vx[1];
2822  grad_vz[0] += grad_vx[2];
2823  grad_vz[1] += grad_vy[2];
2824 
2825  // compute smagorinsky term
2826  double aux = 0.0;
2827  for (unsigned int comp = 0; comp < TDim; comp++)
2828  {
2829  aux += grad_vx[comp] * grad_vx[comp];
2830  aux += grad_vy[comp] * grad_vy[comp];
2831  aux += grad_vz[comp] * grad_vz[comp];
2832  }
2833  aux *= 0.5;
2834  if (aux < 0.0)
2835  aux = 0.0;
2836  double turbulent_viscosity = Cs * h * h * sqrt(aux) ;
2837 
2838  mViscosity[i_node] = turbulent_viscosity + MolecularViscosity;
2839  }
2840  mr_matrix_container.WriteScalarToDatabase(VISCOSITY, mViscosity, rNodes);
2841  KRATOS_CATCH("");
2842  }
2843 
2844  void ApplySmagorinsky2D(double MolecularViscosity, double Cs)
2845  {
2846  KRATOS_TRY
2847  ModelPart::NodesContainerType &rNodes = mr_model_part.Nodes();
2848  // calculating the RHS
2849  array_1d<double, TDim> grad_vx;
2850  array_1d<double, TDim> grad_vy;
2851 
2852  int n_nodes = rNodes.size();
2853  mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes);
2854  array_1d<double, TDim> stab_high;
2855 
2856  #pragma omp parallel for private(grad_vx, grad_vy)
2857  for (int i_node = 0; i_node < n_nodes; i_node++)
2858  {
2859  // set to zero the gradients
2860  for (unsigned int comp = 0; comp < TDim; comp++)
2861  {
2862  grad_vx[comp] = 0.0;
2863  grad_vy[comp] = 0.0;
2864  }
2865  // compute node by node the gradients
2866  const array_1d<double, TDim> &U_i = mvel_n1[i_node];
2867  const double h = mHmin[i_node];
2868  const double m_inv = mr_matrix_container.GetInvertedMass()[i_node];
2869  for (unsigned int csr_index = mr_matrix_container.GetRowStartIndex()[i_node]; csr_index != mr_matrix_container.GetRowStartIndex()[i_node + 1]; csr_index++)
2870  {
2871  unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index];
2872  const array_1d<double, TDim> &U_j = mvel_n1[j_neighbour];
2873  CSR_Tuple &edge_ij = mr_matrix_container.GetEdgeValues()[csr_index];
2874  edge_ij.Add_grad_p(grad_vx, U_i[0], U_j[0]);
2875  edge_ij.Add_grad_p(grad_vy, U_i[1], U_j[1]);
2876  }
2877  // finalize computation of the gradients
2878  // set to zero the gradients
2879  for (unsigned int comp = 0; comp < TDim; comp++)
2880  {
2881  grad_vx[comp] *= m_inv;
2882  grad_vy[comp] *= m_inv;
2883  }
2884  // symmetrize and multiply by 2
2885  grad_vx[0] *= 2.0;
2886  grad_vy[1] *= 2.0;
2887  grad_vx[1] += grad_vy[0];
2888  grad_vy[0] += grad_vx[1];
2889  // compute smagorinsky term
2890  double aux = 0.0;
2891  for (unsigned int comp = 0; comp < TDim; comp++)
2892  {
2893  aux += grad_vx[comp] * grad_vx[comp];
2894  aux += grad_vy[comp] * grad_vy[comp];
2895  }
2896  aux *= 0.5;
2897  if (aux < 0.0)
2898  aux = 0.0;
2899  double turbulent_viscosity = Cs * h * h * sqrt(aux) ;
2900 
2901  mViscosity[i_node] = turbulent_viscosity + MolecularViscosity;
2902  }
2903  mr_matrix_container.WriteScalarToDatabase(VISCOSITY, mViscosity, rNodes);
2904  KRATOS_CATCH("");
2905  }
2906 
2907  void Add_Effective_Inverse_Multiply(
2908  CalcVectorType &destination,
2909  const CalcVectorType &origin1,
2910  const double value,
2911  const ValuesVectorType &mass,
2912  const ValuesVectorType &diag_stiffness,
2913  const CalcVectorType &origin)
2914  {
2915  KRATOS_TRY
2916  int loop_size = destination.size();
2917 
2918  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
2919  array_1d<double, TDim> &dest = destination[i_node];
2920  const double m = mass[i_node];
2921  const double d = diag_stiffness[i_node];
2922  const array_1d<double, TDim> &origin_vec1 = origin1[i_node];
2923  const array_1d<double, TDim> &origin_value = origin[i_node];
2924 
2925  for (unsigned int comp = 0; comp < TDim; comp++)
2926  dest[comp] = value / (m + value * d) * (m / value * origin_vec1[comp] + origin_value[comp]);
2927  });
2928 
2929  KRATOS_CATCH("")
2930  }
2931  };
2932 } // namespace Kratos
2933 #undef SYMM_PRESS
2934 #endif // KRATOS_EDGEBASED_LEVELSET_FLUID_SOLVER_H_INCLUDED defined
cfd_variables.h
Kratos::EdgeBasedLevelSet
Definition: edgebased_levelset.h:39
Kratos::EdgeBasedLevelSet::TSystemVectorType
TSparseSpace::VectorType TSystemVectorType
Definition: edgebased_levelset.h:54
Kratos::EdgeBasedLevelSet::EdgeBasedLevelSet
EdgeBasedLevelSet(MatrixContainer &mr_matrix_container, ModelPart &mr_model_part, const double viscosity, const double density, bool use_mass_correction, double stabdt_pressure_factor, double stabdt_convection_factor, double tau2_factor, bool assume_constant_dp)
Definition: edgebased_levelset.h:60
Kratos::EdgeBasedLevelSet::UpdateFixedVelocityValues
void UpdateFixedVelocityValues()
Definition: edgebased_levelset.h:384
Kratos::EdgeBasedLevelSet::PushFreeSurface
void PushFreeSurface()
Definition: edgebased_levelset.h:2024
Kratos::EdgeBasedLevelSet::ComputeWetVolume
double ComputeWetVolume()
Definition: edgebased_levelset.h:1956
Kratos::EdgeBasedLevelSet::ComputeTimeStep
double ComputeTimeStep(const double CFLNumber, const double MaxDt)
Definition: edgebased_levelset.h:292
Kratos::EdgeBasedLevelSet::Initialize
void Initialize()
Definition: edgebased_levelset.h:97
Kratos::EdgeBasedLevelSet::DiscreteVolumeCorrection
void DiscreteVolumeCorrection(double expected_volume, double measured_volume)
Definition: edgebased_levelset.h:1982
Kratos::EdgeBasedLevelSet::IndicesVectorType
vector< unsigned int > IndicesVectorType
Definition: edgebased_levelset.h:46
Kratos::EdgeBasedLevelSet::~EdgeBasedLevelSet
~EdgeBasedLevelSet()
Definition: edgebased_levelset.h:92
Kratos::EdgeBasedLevelSet::CalculatePorousResistanceLaw
void CalculatePorousResistanceLaw(unsigned int res_law)
Definition: edgebased_levelset.h:2159
Kratos::EdgeBasedLevelSet::Clear
void Clear()
Definition: edgebased_levelset.h:1654
Kratos::EdgeBasedLevelSet::SizeType
std::size_t SizeType
Definition: edgebased_levelset.h:56
Kratos::EdgeBasedLevelSet::ConvectDistance
void ConvectDistance()
Definition: edgebased_levelset.h:1693
Kratos::EdgeBasedLevelSet::ApplyVelocityBC
void ApplyVelocityBC(CalcVectorType &VelArray)
Definition: edgebased_levelset.h:1113
Kratos::EdgeBasedLevelSet::EdgesVectorType
vector< CSR_Tuple > EdgesVectorType
Definition: edgebased_levelset.h:43
Kratos::EdgeBasedLevelSet::MarkInternalNodes
void MarkInternalNodes()
Definition: edgebased_levelset.h:1501
Kratos::EdgeBasedLevelSet::ComputeVolumeVariation
double ComputeVolumeVariation()
Definition: edgebased_levelset.h:1928
Kratos::EdgeBasedLevelSet::MarkInternalAndMixedNodes
void MarkInternalAndMixedNodes()
Definition: edgebased_levelset.h:1472
Kratos::EdgeBasedLevelSet::ChangeSignToDistance
void ChangeSignToDistance()
Definition: edgebased_levelset.h:1396
Kratos::EdgeBasedLevelSet::SolveStep2
void SolveStep2(typename TLinearSolver::Pointer pLinearSolver)
Definition: edgebased_levelset.h:669
Kratos::EdgeBasedLevelSet::ComputeBoundedTimeStep
double ComputeBoundedTimeStep(const double CFLNumber, const double MaxDt)
Definition: edgebased_levelset.h:2055
Kratos::EdgeBasedLevelSet::ActivateWallResistance
void ActivateWallResistance(double Ywall)
Definition: edgebased_levelset.h:1922
Kratos::EdgeBasedLevelSet::CalcVectorType
vector< array_1d< double, TDim > > CalcVectorType
Definition: edgebased_levelset.h:48
Kratos::EdgeBasedLevelSet::SolveStep3
void SolveStep3()
Definition: edgebased_levelset.h:1021
Kratos::EdgeBasedLevelSet::CheckDistanceConvection
bool CheckDistanceConvection()
Definition: edgebased_levelset.h:1866
Kratos::EdgeBasedLevelSet::SetShockCapturingCoefficient
void SetShockCapturingCoefficient(double coeff)
Definition: edgebased_levelset.h:284
Kratos::EdgeBasedLevelSet::SolveStep1
void SolveStep1()
Definition: edgebased_levelset.h:411
Kratos::EdgeBasedLevelSet::ApplySmagorinsky
void ApplySmagorinsky(double MolecularViscosity, double Cs)
Definition: edgebased_levelset.h:373
Kratos::EdgeBasedLevelSet::TSystemMatrixType
TSparseSpace::MatrixType TSystemMatrixType
Definition: edgebased_levelset.h:53
Kratos::EdgeBasedLevelSet::MarkExternalAndMixedNodes
void MarkExternalAndMixedNodes()
Definition: edgebased_levelset.h:1443
Kratos::EdgeBasedLevelSet::ValuesVectorType
vector< double > ValuesVectorType
Definition: edgebased_levelset.h:50
Kratos::EdgeBasedLevelSet::ExtrapolateValues
void ExtrapolateValues(unsigned int extrapolation_layers)
Definition: edgebased_levelset.h:1203
Kratos::EdgeBasedLevelSet::ReduceTimeStep
void ReduceTimeStep(ModelPart &rModelPart, double NewTime)
Definition: edgebased_levelset.h:1856
Kratos::EdgeBasedLevelSet::MarkNodesByDistance
void MarkNodesByDistance(double min, double max)
Definition: edgebased_levelset.h:1411
Kratos::EdgeBasedLevelSet::CSR_Tuple
EdgesStructureType< TDim > CSR_Tuple
Definition: edgebased_levelset.h:42
Kratos::EdgeBasedLevelSet::CalculateNormals
void CalculateNormals(ModelPart::ConditionsContainerType &rConditions)
Definition: edgebased_levelset.h:1528
Kratos::EdgeBasedLevelSet::SaveScalarVariableToOldStep
void SaveScalarVariableToOldStep(Variable< double > &rVar)
Definition: edgebased_levelset.h:1429
Kratos::EdgeBasedLevelSet::CalculateRHS
void CalculateRHS(const CalcVectorType &vel, const ValuesVectorType &pressure, const CalcVectorType &convective_velocity, CalcVectorType &rhs, ValuesVectorType &diag_stiffness)
Definition: edgebased_levelset.h:562
Kratos::EdgesStructureType
Definition: edge_data.h:42
Kratos::EdgesStructureType::Sub_ViscousContribution
void Sub_ViscousContribution(array_1d< double, TDim > &destination, const array_1d< double, TDim > &U_i, const double &nu_i, const array_1d< double, TDim > &U_j, const double &nu_j)
Definition: edge_data.h:362
Kratos::EdgesStructureType::Sub_ConvectiveContribution
void Sub_ConvectiveContribution(array_1d< double, TDim > &destination, const array_1d< double, TDim > &a_i, const array_1d< double, TDim > &U_i, const array_1d< double, TDim > &a_j, const array_1d< double, TDim > &U_j)
Definition: edge_data.h:167
Kratos::EdgesStructureType::Add_grad_p
void Add_grad_p(array_1d< double, TDim > &destination, const double &p_i, const double &p_j)
Definition: edge_data.h:99
Kratos::EdgesStructureType::Sub_StabContribution
void Sub_StabContribution(array_1d< double, TDim > &destination, const double tau, const double beta, const array_1d< double, TDim > &stab_low, const array_1d< double, TDim > &stab_high)
Definition: edge_data.h:330
Kratos::EdgesStructureType::CalculateConvectionStabilization_HIGH
void CalculateConvectionStabilization_HIGH(array_1d< double, TDim > &stab_high, const array_1d< double, TDim > &a_i, const array_1d< double, TDim > &pi_i, const array_1d< double, TDim > &a_j, const array_1d< double, TDim > &pi_j)
Definition: edge_data.h:265
Kratos::EdgesStructureType::Add_D_v
void Add_D_v(double &destination, const array_1d< double, TDim > &v_i, const array_1d< double, TDim > &v_j)
Definition: edge_data.h:78
Kratos::EdgesStructureType::Sub_grad_p
void Sub_grad_p(array_1d< double, TDim > &destination, const double &p_i, const double &p_j)
Definition: edge_data.h:105
Kratos::EdgesStructureType::CalculateConvectionStabilization_LOW
void CalculateConvectionStabilization_LOW(array_1d< double, TDim > &stab_low, const array_1d< double, TDim > &a_i, const array_1d< double, TDim > &U_i, const array_1d< double, TDim > &a_j, const array_1d< double, TDim > &U_j)
Definition: edge_data.h:240
Kratos::Geometry
Geometry base class.
Definition: geometry.h:71
Kratos::GlobalPointersVector
This class is a vector which stores global pointers.
Definition: global_pointers_vector.h:61
Kratos::GlobalPointersVector::iterator
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: global_pointers_vector.h:79
Kratos::GlobalPointersVector::push_back
void push_back(TPointerType x)
Definition: global_pointers_vector.h:322
Kratos::GlobalPointersVector::begin
iterator begin()
Definition: global_pointers_vector.h:221
Kratos::GlobalPointersVector::end
iterator end()
Definition: global_pointers_vector.h:229
Kratos::IndexPartition
This class is useful for index iteration over containers.
Definition: parallel_utilities.h:451
Kratos::IndexPartition::for_each
void for_each(TUnaryFunction &&f)
Definition: parallel_utilities.h:514
Kratos::MathUtils::CrossProduct
static T CrossProduct(const T &a, const T &b)
Performs the vector product of the two input vectors a,b.
Definition: math_utils.h:762
Kratos::MatrixContainer
Definition: edge_data.h:381
Kratos::MatrixContainer::WriteVectorToDatabase
void WriteVectorToDatabase(Variable< array_1d< double, 3 >> &rVariable, CalcVectorType &rOrigin, ModelPart::NodesContainerType &rNodes)
Definition: edge_data.h:905
Kratos::MatrixContainer::AssignVectorToVector
void AssignVectorToVector(const CalcVectorType &origin, CalcVectorType &destination)
Definition: edge_data.h:1082
Kratos::MatrixContainer::GetLumpedMass
ValuesVectorType & GetLumpedMass()
Definition: edge_data.h:420
Kratos::MatrixContainer::FillOldScalarFromDatabase
void FillOldScalarFromDatabase(Variable< double > &rVariable, ValuesVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data.h:879
Kratos::MatrixContainer::SetToZero
void SetToZero(EdgesVectorType &data_vector)
Definition: edge_data.h:1038
Kratos::MatrixContainer::GetColumnIndex
IndicesVectorType & GetColumnIndex()
Definition: edge_data.h:410
Kratos::MatrixContainer::GetInvertedMass
ValuesVectorType & GetInvertedMass()
Definition: edge_data.h:425
Kratos::MatrixContainer::GetRowStartIndex
IndicesVectorType & GetRowStartIndex()
Definition: edge_data.h:415
Kratos::MatrixContainer::FillVectorFromDatabase
void FillVectorFromDatabase(Variable< array_1d< double, 3 >> &rVariable, CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data.h:790
Kratos::MatrixContainer::FillScalarFromDatabase
void FillScalarFromDatabase(Variable< double > &rVariable, ValuesVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data.h:851
Kratos::MatrixContainer::WriteScalarToDatabase
void WriteScalarToDatabase(Variable< double > &rVariable, ValuesVectorType &rOrigin, ModelPart::NodesContainerType &rNodes)
Definition: edge_data.h:933
Kratos::MatrixContainer::Add_Minv_value
void Add_Minv_value(CalcVectorType &destination, const CalcVectorType &origin1, const double value, const ValuesVectorType &Minv_vec, const CalcVectorType &origin)
Definition: edge_data.h:963
Kratos::MatrixContainer::FillOldVectorFromDatabase
void FillOldVectorFromDatabase(Variable< array_1d< double, 3 >> &rVariable, CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data.h:821
Kratos::MatrixContainer::FillCoordinatesFromDatabase
void FillCoordinatesFromDatabase(CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data.h:759
Kratos::MatrixContainer::GetHmin
ValuesVectorType & GetHmin()
Definition: edge_data.h:435
Kratos::MatrixContainer::GetNumberEdges
unsigned int GetNumberEdges()
Definition: edge_data.h:400
Kratos::MatrixContainer::GetEdgeValues
EdgesVectorType & GetEdgeValues()
Definition: edge_data.h:405
Kratos::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::ModelPart::ConditionsContainerType
MeshType::ConditionsContainerType ConditionsContainerType
Condintions container. A vector set of Conditions with their Id's as key.
Definition: model_part.h:183
Kratos::ModelPart::OverwriteSolutionStepData
void OverwriteSolutionStepData(IndexType SourceSolutionStepIndex, IndexType DestinationSourceSolutionStepIndex)
Definition: model_part.cpp:180
Kratos::ModelPart::NodesBegin
NodeIterator NodesBegin(IndexType ThisIndex=0)
Definition: model_part.h:487
Kratos::ModelPart::Conditions
ConditionsContainerType & Conditions(IndexType ThisIndex=0)
Definition: model_part.h:1381
Kratos::ModelPart::GetProcessInfo
ProcessInfo & GetProcessInfo()
Definition: model_part.h:1746
Kratos::ModelPart::NodesContainerType
MeshType::NodesContainerType NodesContainerType
Nodes container. Which is a vector set of nodes with their Id's as key.
Definition: model_part.h:128
Kratos::ModelPart::NodesEnd
NodeIterator NodesEnd(IndexType ThisIndex=0)
Definition: model_part.h:497
Kratos::ModelPart::Nodes
NodesContainerType & Nodes(IndexType ThisIndex=0)
Definition: model_part.h:507
Kratos::Node
This class defines the node.
Definition: node.h:65
Kratos::OpenMPUtils::DivideInPartitions
static void DivideInPartitions(const int NumTerms, const int NumThreads, PartitionVector &Partitions)
Divide an array of length NumTerms between NumThreads threads.
Definition: openmp_utils.h:158
Kratos::OpenMPUtils::ThisThread
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
Kratos::ParallelUtilities::GetNumThreads
static int GetNumThreads()
Returns the current number of threads.
Definition: parallel_utilities.cpp:34
Kratos::ProcessInfo
ProcessInfo holds the current value of different solution parameters.
Definition: process_info.h:59
Kratos::ProcessInfo::SetCurrentTime
void SetCurrentTime(double NewTime)
Definition: process_info.cpp:63
Kratos::array_1d< double, 3 >
KRATOS_THROW_ERROR
#define KRATOS_THROW_ERROR(ExceptionType, ErrorMessage, MoreInfo)
Definition: define.h:77
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_WATCH
#define KRATOS_WATCH(variable)
Definition: define.h:806
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
deprecated_variables.h
free_surface_application.h
global_pointer_variables.h
KRATOS_ERROR_IF
#define KRATOS_ERROR_IF(conditional)
Definition: exception.h:162
model_part.h
DEM_benchmarks.dt
dt
Definition: DEM_benchmarks.py:173
DEM_benchmarks.end
end
Definition: DEM_benchmarks.py:180
Kratos::GeometricalTransformationUtilities::VectorType
Vector VectorType
Definition: geometrical_transformation_utilities.h:56
Kratos::GeometricalTransformationUtilities::MatrixType
Matrix MatrixType
Definition: geometrical_transformation_utilities.h:55
Kratos::GeometryFunctions::max
static double max(double a, double b)
Definition: GeometryFunctions.h:79
Kratos::GeometryFunctions::min
static double min(double a, double b)
Definition: GeometryFunctions.h:71
Kratos::Python::GetValue
Parameters GetValue(Parameters &rParameters, const std::string &rEntry)
Definition: add_kratos_parameters_to_python.cpp:53
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::norm_2
TExpressionType::data_type norm_2(AMatrix::MatrixExpression< TExpressionType, TCategory > const &TheExpression)
Definition: amatrix_interface.h:625
Kratos::inner_prod
TExpression1Type::data_type inner_prod(AMatrix::MatrixExpression< TExpression1Type, TCategory1 > const &First, AMatrix::MatrixExpression< TExpression2Type, TCategory2 > const &Second)
Definition: amatrix_interface.h:592
Kratos::PointVector
std::vector< PointTypePointer > PointVector
Definition: internal_variables_interpolation_process.h:53
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
Kratos::double
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
Kratos::PointIterator
PointVector::iterator PointIterator
Definition: internal_variables_interpolation_process.h:54
Kratos::PointType
Point PointType
Geometric definitions.
Definition: mortar_classes.h:36
bombardelli_test.coeff
list coeff
Definition: bombardelli_test.py:41
edgebased_PureConvection.dist
float dist
Definition: edgebased_PureConvection.py:89
edgebased_var.use_mass_correction
bool use_mass_correction
Definition: edgebased_var.py:13
edgebased_var.viscosity
float viscosity
Definition: edgebased_var.py:8
edgebased_var.extrapolation_layers
int extrapolation_layers
Definition: edgebased_var.py:15
face_heat.density
float density
Definition: face_heat.py:56
generate_axisymmetric_navier_stokes_element.y
y
Other simbols definition.
Definition: generate_axisymmetric_navier_stokes_element.py:54
generate_convection_diffusion_explicit_element.tau
tau
Definition: generate_convection_diffusion_explicit_element.py:115
generate_droplet_dynamics.h
h
Definition: generate_droplet_dynamics.py:91
generate_frictional_mortar_condition.rhs
rhs
Definition: generate_frictional_mortar_condition.py:297
generate_stokes_twofluid_element.a
a
Definition: generate_stokes_twofluid_element.py:77
generate_total_lagrangian_mixed_volumetric_strain_element.b
b
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:31
generate_total_lagrangian_mixed_volumetric_strain_element.n_nodes
int n_nodes
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:15
isotropic_damage_automatic_differentiation.nu
nu
Definition: isotropic_damage_automatic_differentiation.py:135
ode_solve.t
int t
Definition: ode_solve.py:392
ode_solve.d
int d
Definition: ode_solve.py:397
ode_solve.n
int n
manufactured solution and derivatives (u=0 at z=0 dudz=0 at z=domain_height)
Definition: ode_solve.py:402
pure_conduction.vel
vel
Definition: pure_conduction.py:76
quadrature.k
int k
Definition: quadrature.py:595
quadrature.j
int j
Definition: quadrature.py:648
rotating_cone.temp
float temp
Definition: rotating_cone.py:85
rotatingcone_PureConvectionBenchmarking.delta_t
float delta_t
Definition: rotatingcone_PureConvectionBenchmarking.py:129
run_marine_rain_substepping.m
int m
Definition: run_marine_rain_substepping.py:8
script_ELASTIC.body_force
body_force
Definition: script_ELASTIC.py:102
script.current_id
int current_id
Output settings end ####.
Definition: script.py:194
sensitivityMatrix.B
B
Definition: sensitivityMatrix.py:76
tensors::i
integer i
Definition: TensorModule.f:17
tensors::sign1
integer function sign1(a)
Definition: TensorModule.f:872
test_pureconvectionsolver_benchmarking.pressure
float pressure
Definition: test_pureconvectionsolver_benchmarking.py:101
testing.run_cpp_mpi_tests.e
e
Definition: run_cpp_mpi_tests.py:31
ulf_PGLASS.factor
float factor
for node in (self.combined_model_part).Nodes: pold = node.GetSolutionStepValue(PRESSURE,...
Definition: ulf_PGLASS.py:254
Generated on Sat Feb 3 2024 02:01:53 for KratosMultiphysics by doxygen 1.9.1