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