From 93e67f416db2a2f0e1956b12dd3b8a428285930a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Joaqu=C3=ADn=20Iraz=C3=A1bal=20Gonz=C3=A1lez?= Date: Mon, 9 Oct 2023 12:49:24 +0200 Subject: [PATCH] Fix wall law and clean-up (#11670) --- .../custom_utilities/edgebased_levelset.h | 53 ++++--- .../edgebased_levelset_substep.h | 144 ++++++------------ .../edgebased_levelset_solver.py | 14 +- .../python_scripts/free_surface_analysis.py | 7 +- 4 files changed, 86 insertions(+), 132 deletions(-) diff --git a/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset.h b/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset.h index ce24da2e2e5c..d2dd6a5bbbc9 100644 --- a/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset.h +++ b/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset.h @@ -77,8 +77,6 @@ namespace Kratos for (ModelPart::NodesContainerType::iterator it = mr_model_part.NodesBegin(); it != mr_model_part.NodesEnd(); it++) it->FastGetSolutionStepValue(VISCOSITY) = viscosity; - mMolecularViscosity = viscosity; - mRho = density; mdelta_t_avg = 1000.0; @@ -392,7 +390,8 @@ namespace Kratos mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes); int fixed_size = mFixedVelocities.size(); -#pragma omp parallel for firstprivate(fixed_size) + + #pragma omp parallel for firstprivate(fixed_size) for (int i_velocity = 0; i_velocity < fixed_size; i_velocity++) { unsigned int i_node = mFixedVelocities[i_velocity]; @@ -447,7 +446,7 @@ namespace Kratos double stabdt_convection_factor = mstabdt_convection_factor; double tau2_factor = mtau2_factor; -#pragma omp parallel for firstprivate(time_inv_avg, stabdt_pressure_factor, stabdt_convection_factor, tau2_factor) + #pragma omp parallel for firstprivate(time_inv_avg, stabdt_pressure_factor, stabdt_convection_factor, tau2_factor) for (int i_node = 0; i_node < n_nodes; i_node++) { double &h_avg_i = mHavg[i_node]; @@ -576,7 +575,8 @@ namespace Kratos array_1d stab_high; double inverse_rho = 1.0 / mRho; -#pragma omp parallel for private(stab_low, stab_high) + + #pragma omp parallel for private(stab_low, stab_high) for (int i_node = 0; i_node < n_nodes; i_node++) { double dist = mdistances[i_node]; @@ -974,12 +974,11 @@ namespace Kratos mPn1[i_node] += dp[i_node] * scaling_factors[i_node]; }); - // write pressure and density to Kratos + // write pressure to Kratos mr_matrix_container.WriteScalarToDatabase(PRESSURE, mPn1, rNodes); // compute pressure proj for the next step - -#pragma omp parallel for private(work_array) + #pragma omp parallel for private(work_array) for (int i_node = 0; i_node < n_nodes; i_node++) { array_1d &xi_i = mXi[i_node]; @@ -1039,7 +1038,7 @@ namespace Kratos // compute end of step momentum double rho_inv = 1.0 / mRho; -#pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv, factor) + #pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv, factor) for (int i_node = 0; i_node < n_nodes; i_node++) { double dist = mdistances[i_node]; @@ -1082,7 +1081,7 @@ namespace Kratos // calculate the error on the divergence if (muse_mass_correction == true) { -#pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv) + #pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv) for (int i_node = 0; i_node < n_nodes; i_node++) { const double dist = mdistances[i_node]; @@ -1119,7 +1118,8 @@ namespace Kratos { // apply conditions on corner edges int edge_size = medge_nodes_direction.size(); -#pragma omp parallel for firstprivate(edge_size) + + #pragma omp parallel for firstprivate(edge_size) for (int i = 0; i < edge_size; i++) { int i_node = medge_nodes[i]; @@ -1152,7 +1152,8 @@ namespace Kratos // slip condition int slip_size = mSlipBoundaryList.size(); -#pragma omp parallel for firstprivate(slip_size) + + #pragma omp parallel for firstprivate(slip_size) for (int i_slip = 0; i_slip < slip_size; i_slip++) { unsigned int i_node = mSlipBoundaryList[i_slip]; @@ -1177,7 +1178,8 @@ namespace Kratos // fixed condition int fixed_size = mFixedVelocities.size(); -#pragma omp parallel for firstprivate(fixed_size) + + #pragma omp parallel for firstprivate(fixed_size) for (int i_velocity = 0; i_velocity < fixed_size; i_velocity++) { unsigned int i_node = mFixedVelocities[i_velocity]; @@ -1723,7 +1725,8 @@ namespace Kratos // calculating the convective projection array_1d a_i; array_1d a_j; -#pragma omp parallel for private(a_i, a_j) + + #pragma omp parallel for private(a_i, a_j) for (int i_node = 0; i_node < n_nodes; i_node++) { array_1d &pi_i = mPiConvection[i_node]; @@ -2204,7 +2207,6 @@ namespace Kratos } private: - double mMolecularViscosity; MatrixContainer &mr_matrix_container; ModelPart &mr_model_part; @@ -2357,7 +2359,6 @@ namespace Kratos { double &h_i = mHavg[i_node]; double &m_i = mr_matrix_container.GetLumpedMass()[i_node]; - // double& rho_i = mRho[i_node]; h_i = sqrt(2.0 * m_i); } @@ -2368,7 +2369,6 @@ namespace Kratos { double &h_i = mHavg[i_node]; double &m_i = mr_matrix_container.GetLumpedMass()[i_node]; - // double& rho_i = mRho[i_node]; h_i = pow(6.0 * m_i, 1.0 / 3.0); } @@ -2411,7 +2411,8 @@ namespace Kratos double stab_high; array_1d a_i; array_1d a_j; -#pragma omp parallel for private(stab_low, stab_high, a_i, a_j) + + #pragma omp parallel for private(stab_low, stab_high, a_i, a_j) for (int i_node = 0; i_node < n_nodes; i_node++) { double &rhs_i = rhs[i_node]; @@ -2707,18 +2708,21 @@ namespace Kratos double y_plus_incercept = 10.9931899; unsigned int itmax = 100; - KRATOS_ERROR_IF(mViscosity[0] == 0) << "it is not possible to use the wall law with 0 viscosity" << std::endl; - // slip condition int slip_size = mSlipBoundaryList.size(); -#pragma omp parallel for firstprivate(slip_size, B, toll, ym, y_plus_incercept, itmax) + + #pragma omp parallel for firstprivate(slip_size, B, toll, ym, y_plus_incercept, itmax) for (int i_slip = 0; i_slip < slip_size; i_slip++) { unsigned int i_node = mSlipBoundaryList[i_slip]; double dist = mdistances[i_node]; - const double nu = mViscosity[i_node]; if (dist <= 0.0) { + + KRATOS_ERROR_IF(mViscosity[i_node] == 0) << "it is not possible to use the wall law with 0 viscosity" << std::endl; + + const double nu = mViscosity[i_node]; + // array_1d& rhs_i = rhs[i_node]; const array_1d &U_i = vel[i_node]; const array_1d &an_i = mSlipNormal[i_node]; @@ -2775,7 +2779,8 @@ namespace Kratos int n_nodes = rNodes.size(); mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes); array_1d stab_high; -#pragma omp parallel for private(grad_vx, grad_vy, grad_vz) + + #pragma omp parallel for private(grad_vx, grad_vy, grad_vz) for (int i_node = 0; i_node < n_nodes; i_node++) { // set to zero the gradients @@ -2848,7 +2853,7 @@ namespace Kratos mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes); array_1d stab_high; -#pragma omp parallel for private(grad_vx, grad_vy) + #pragma omp parallel for private(grad_vx, grad_vy) for (int i_node = 0; i_node < n_nodes; i_node++) { // set to zero the gradients diff --git a/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset_substep.h b/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset_substep.h index 89d65e385676..07dc8a06a4ca 100644 --- a/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset_substep.h +++ b/applications/FreeSurfaceApplication/custom_utilities/edgebased_levelset_substep.h @@ -85,7 +85,6 @@ namespace Kratos for (ModelPart::NodesContainerType::iterator it = mr_model_part.NodesBegin(); it != mr_model_part.NodesEnd(); it++) it->FastGetSolutionStepValue(VISCOSITY) = viscosity; - mMolecularViscosity = viscosity; mRho = density; mdelta_t_avg = 1000.0; max_dt = 1.0; @@ -242,7 +241,7 @@ namespace Kratos OpenMPUtils::DivideInPartitions(n_nodes, number_of_threads, row_partition); for (int k = 0; k < number_of_threads; k++) { -#pragma omp parallel + #pragma omp parallel if (OpenMPUtils::ThisThread() == k) { for (int i_node = static_cast(row_partition[k]); i_node < static_cast(row_partition[k + 1]); i_node++) @@ -348,7 +347,7 @@ namespace Kratos Vector dt_vec(n_proc, 1e10); Vector dt_avg_novisc_vec(n_proc, 1e10); -#pragma omp parallel for firstprivate(n_nodes) + #pragma omp parallel for firstprivate(n_nodes) for (int i_node = 0; i_node < n_nodes; i_node++) { unsigned int my_id = OpenMPUtils::ThisThread(); @@ -447,7 +446,8 @@ namespace Kratos KRATOS_TRY // read velocity and pressure data from Kratos int fixed_size = mFixedVelocities.size(); -#pragma omp parallel for firstprivate(fixed_size) + + #pragma omp parallel for firstprivate(fixed_size) for (int i_velocity = 0; i_velocity < fixed_size; i_velocity++) { unsigned int i_node = mFixedVelocities[i_velocity]; @@ -487,8 +487,7 @@ namespace Kratos double stabdt_pressure_factor = mstabdt_pressure_factor; double stabdt_convection_factor = mstabdt_convection_factor; -// double tau2_factor = mtau2_factor; -#pragma omp parallel for firstprivate(time_inv_avg, stabdt_pressure_factor, stabdt_convection_factor) + #pragma omp parallel for firstprivate(time_inv_avg, stabdt_pressure_factor, stabdt_convection_factor) for (int i_node = 0; i_node < n_nodes; i_node++) { double &h_avg_i = mHavg[i_node]; @@ -558,26 +557,6 @@ namespace Kratos pi_i[l_comp] *= m_inv; }); -#ifdef DEBUG_OUTPUT - KRATOS_WATCH("before RK of step1 - new") - double aux_v = 0.0; - for (int i_node = 0; i_node < mvel_n1.size(); i_node++) - aux_v += inner_prod(mvel_n1[i_node], mvel_n1[i_node]); - double aux_oldv = 0.0; - for (int i_node = 0; i_node < mvel_n1.size(); i_node++) - aux_oldv += inner_prod(mvel_n[i_node], mvel_n[i_node]); - double aux_pi = 0.0; - for (int i_node = 0; i_node < mvel_n1.size(); i_node++) - aux_pi += inner_prod(mPi[i_node], mPi[i_node]); - - KRATOS_WATCH(inner_prod(mPn, mPn)); - KRATOS_WATCH(aux_v); - KRATOS_WATCH(aux_oldv); - KRATOS_WATCH(aux_pi); - KRATOS_WATCH(inner_prod(mdistances, mdistances)); - KRATOS_WATCH(inner_prod(mViscosity, mViscosity)); -#endif - CalcVectorType auxn = mvel_n; double n_substeps = mnumsubsteps + 1; @@ -662,22 +641,6 @@ namespace Kratos } } -#ifdef DEBUG_OUTPUT - KRATOS_WATCH("end of step1 - new") - aux_v = 0.0; - for (int i_node = 0; i_node < mvel_n1.size(); i_node++) - aux_v += inner_prod(mvel_n1[i_node], mvel_n1[i_node]); - double aux_xi = 0.0; - for (int i_node = 0; i_node < mvel_n1.size(); i_node++) - aux_xi += inner_prod(mXi[i_node], mXi[i_node]); - - KRATOS_WATCH(inner_prod(mPn, mPn)); - KRATOS_WATCH(inner_prod(mdistances, mdistances)); - KRATOS_WATCH(inner_prod(mViscosity, mViscosity)); - - KRATOS_WATCH(aux_v); - KRATOS_WATCH(aux_xi); -#endif KRATOS_CATCH("") } //********************************************************************* @@ -696,7 +659,8 @@ namespace Kratos array_1d stab_low; array_1d stab_high; double inverse_rho = 1.0 / mRho; -#pragma omp parallel for private(stab_low, stab_high) + + #pragma omp parallel for private(stab_low, stab_high) for (int i_node = 0; i_node < n_nodes; i_node++) { double dist = mdistances[i_node]; @@ -789,7 +753,8 @@ namespace Kratos // Re-generate a container with LAYER 0 and LAYER 1 after convection of the free surface layer_limits[0] = 0; -#pragma omp parallel for + + #pragma omp parallel for for (int i_node = 0; i_node < static_cast(mr_model_part.Nodes().size()); i_node++) { if (mdistances[i_node] < 0.0) @@ -800,7 +765,7 @@ namespace Kratos unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index]; if (mdistances[j_neighbour] >= 0.0 && mis_visited[i_node] == 0) { -#pragma omp critical + #pragma omp critical layers[++layer_counter] = i_node; mis_visited[i_node] = 1; @@ -832,7 +797,7 @@ namespace Kratos int return_value = 0; // on the first layer outside the pressure is set to a value such that on the free surface the pressure is approx 0 -#pragma omp parallel for + #pragma omp parallel for for (int iii = static_cast(layer_limits[1]); iii < static_cast(layer_limits[2]); iii++) { unsigned int i_node = layers[iii]; @@ -1065,11 +1030,11 @@ namespace Kratos mPn1[i_node] += dp[i_node] * scaling_factors[i_node]; }); - // write pressure and density to Kratos + // write pressure to Kratos mr_matrix_container.WriteScalarToDatabase(PRESSURE, mPn1, rNodes); -// compute pressure proj for the next step -#pragma omp parallel for private(work_array) + // compute pressure proj for the next step + #pragma omp parallel for private(work_array) for (int i_node = 0; i_node < n_nodes; i_node++) { array_1d &xi_i = mXi[i_node]; @@ -1095,21 +1060,6 @@ namespace Kratos } mr_matrix_container.WriteVectorToDatabase(PRESS_PROJ, mXi, rNodes); -#ifdef DEBUG_OUTPUT - - KRATOS_WATCH("end of step2 - new") - double aux_v = 0.0; - for (int i_node = 0; i_node < mvel_n1.size(); i_node++) - aux_v += inner_prod(mvel_n1[i_node], mvel_n1[i_node]); - double aux_xi = 0.0; - for (int i_node = 0; i_node < mvel_n1.size(); i_node++) - aux_xi += inner_prod(mXi[i_node], mXi[i_node]); - - KRATOS_WATCH(inner_prod(mPn1, mPn1)); - KRATOS_WATCH(aux_v); - KRATOS_WATCH(aux_xi); -#endif - return return_value; KRATOS_CATCH("") } @@ -1132,7 +1082,8 @@ namespace Kratos factor = 1.0; // compute end of step momentum double rho_inv = 1.0 / mRho; -#pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv, factor) + + #pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv, factor) for (int i_node = 0; i_node < n_nodes; i_node++) { double dist = mdistances[i_node]; @@ -1181,7 +1132,7 @@ namespace Kratos // calculate the error on the divergence if (muse_mass_correction == true) { -#pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv) + #pragma omp parallel for private(correction) firstprivate(delta_t, rho_inv) for (int i_node = 0; i_node < n_nodes; i_node++) { const double dist = mdistances[i_node]; @@ -1202,17 +1153,6 @@ namespace Kratos } } -#ifdef DEBUG_OUTPUT - - KRATOS_WATCH("end of step 3") - double aux = 0.0; - for (int i_node = 0; i_node < n_nodes; i_node++) - aux += inner_prod(mvel_n1[i_node], mvel_n1[i_node]); - - KRATOS_WATCH(inner_prod(mPn1, mPn1)); - KRATOS_WATCH(aux); -#endif - mr_matrix_container.WriteVectorToDatabase(VELOCITY, mvel_n1, rNodes); KRATOS_CATCH("") @@ -1222,7 +1162,8 @@ namespace Kratos KRATOS_TRY // slip condition int size = mDistanceBoundaryList.size(); -#pragma omp parallel for firstprivate(size) + + #pragma omp parallel for firstprivate(size) for (int i_dist = 0; i_dist < size; i_dist++) { unsigned int i_node = mDistanceBoundaryList[i_dist]; @@ -1271,7 +1212,8 @@ namespace Kratos // slip condition int slip_size = mSlipBoundaryList.size(); -#pragma omp parallel for firstprivate(slip_size) + + #pragma omp parallel for firstprivate(slip_size) for (int i_slip = 0; i_slip < slip_size; i_slip++) { unsigned int i_node = mSlipBoundaryList[i_slip]; @@ -1296,7 +1238,8 @@ namespace Kratos // fixed condition int fixed_size = mFixedVelocities.size(); -#pragma omp parallel for firstprivate(fixed_size) + + #pragma omp parallel for firstprivate(fixed_size) for (int i_velocity = 0; i_velocity < fixed_size; i_velocity++) { unsigned int i_node = mFixedVelocities[i_velocity]; @@ -1330,7 +1273,7 @@ namespace Kratos layer_limits[0] = 0; int layer_counter = -1; -#pragma omp parallel for + #pragma omp parallel for for (int i_node = 0; i_node < static_cast(mr_model_part.Nodes().size()); i_node++) { if (mdistances[i_node] < 0.0) @@ -1341,7 +1284,7 @@ namespace Kratos unsigned int j_neighbour = mr_matrix_container.GetColumnIndex()[csr_index]; if (mdistances[j_neighbour] >= 0.0 && mis_visited[i_node] == 0) { -#pragma omp critical + #pragma omp critical layers[++layer_counter] = i_node; mis_visited[i_node] = 1; @@ -1383,12 +1326,9 @@ namespace Kratos array_1d aux, aux_proj; -// ProcessInfo& CurrentProcessInfo = mr_model_part.GetProcessInfo(); -// double delta_t = CurrentProcessInfo[DELTA_TIME]; - -// fill the pressure projection on the first layer inside the fluid -// by extrapolating from the pressure projection on the layer -1 (the first layer completely inside the domain) -#pragma omp parallel for + // fill the pressure projection on the first layer inside the fluid + // by extrapolating from the pressure projection on the layer -1 (the first layer completely inside the domain) + #pragma omp parallel for for (int i = layer_limits[0]; i < layer_limits[1]; i++) { unsigned int i_node = layers[i]; @@ -1475,7 +1415,7 @@ namespace Kratos // now mark all of the nodes up to the extrapolation layers - 1 for (unsigned int il = 0; il < extrapolation_layers - 1; il++) { -#pragma omp parallel for + #pragma omp parallel for for (int iii = static_cast(layer_limits[il]); iii < static_cast(layer_limits[il + 1]); iii++) { unsigned int i_node = layers[iii]; @@ -1922,7 +1862,8 @@ namespace Kratos // compute wall reduction factor // slip condition int slip_size = mSlipBoundaryList.size(); -#pragma omp parallel for firstprivate(slip_size) + + #pragma omp parallel for firstprivate(slip_size) for (int i_slip = 0; i_slip < slip_size; i_slip++) { unsigned int i_node = mSlipBoundaryList[i_slip]; @@ -1932,7 +1873,8 @@ namespace Kratos std::cout << "max angle between normals found in the model = " << max_angle_overall << std::endl; int edge_size = medge_nodes.size(); -#pragma omp parallel for firstprivate(edge_size) + + #pragma omp parallel for firstprivate(edge_size) for (int i = 0; i < edge_size; i++) { int i_node = medge_nodes[i]; @@ -2168,7 +2110,6 @@ namespace Kratos } private: - double mMolecularViscosity; double mcorner_coefficient; double medge_coefficient; double mmax_dt; @@ -2343,7 +2284,8 @@ namespace Kratos double stab_high; array_1d a_i; array_1d a_j; -#pragma omp parallel for private(stab_low, stab_high, a_i, a_j) + + #pragma omp parallel for private(stab_low, stab_high, a_i, a_j) for (int i_node = 0; i_node < n_nodes; i_node++) { double &rhs_i = rhs[i_node]; @@ -2588,18 +2530,22 @@ namespace Kratos { double ym = mY_wall; - KRATOS_ERROR_IF(mViscosity[0] == 0) << "it is not possible to use the wall law with 0 viscosity" << std::endl; + // slip condition int slip_size = mSlipBoundaryList.size(); -#pragma omp parallel for firstprivate(slip_size, ym) + + #pragma omp parallel for firstprivate(slip_size, ym) for (int i_slip = 0; i_slip < slip_size; i_slip++) { unsigned int i_node = mSlipBoundaryList[i_slip]; double dist = mdistances[i_node]; if (dist <= 0.0) { - double nu = mMolecularViscosity; + KRATOS_ERROR_IF(mViscosity[i_node] == 0) << "it is not possible to use the wall law with 0 viscosity" << std::endl; + + double nu = mViscosity[i_node]; + const array_1d &U_i = vel[i_node]; const array_1d &an_i = mSlipNormal[i_node]; @@ -2635,7 +2581,8 @@ namespace Kratos int n_nodes = rNodes.size(); mr_matrix_container.FillVectorFromDatabase(VELOCITY, mvel_n1, rNodes); array_1d stab_high; -#pragma omp parallel for private(grad_vx, grad_vy, grad_vz) + + #pragma omp parallel for private(grad_vx, grad_vy, grad_vz) for (int i_node = 0; i_node < n_nodes; i_node++) { // set to zero the gradients @@ -2734,7 +2681,8 @@ namespace Kratos // calculating the convective projection array_1d a_i; array_1d a_j; -#pragma omp parallel for private(a_i, a_j) + + #pragma omp parallel for private(a_i, a_j) for (int i_node = 0; i_node < n_nodes; i_node++) { array_1d &pi_i = mPiConvection[i_node]; diff --git a/applications/FreeSurfaceApplication/python_scripts/edgebased_levelset_solver.py b/applications/FreeSurfaceApplication/python_scripts/edgebased_levelset_solver.py index 7e1f92226b91..8d07ef5e7a6a 100644 --- a/applications/FreeSurfaceApplication/python_scripts/edgebased_levelset_solver.py +++ b/applications/FreeSurfaceApplication/python_scripts/edgebased_levelset_solver.py @@ -143,12 +143,12 @@ def AddVariables(self) -> None: self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DISTANCE) self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.BODY_FORCE) self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.PRESS_PROJ) - self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.POROSITY) - self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.VISCOSITY) self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DENSITY) - self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DIAMETER) + self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.VISCOSITY) + self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.POROSITY) self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.LIN_DARCY_COEF) self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.NONLIN_DARCY_COEF) + self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.DIAMETER) self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.NODAL_AREA) self.model_part.AddNodalSolutionStepVariable(KratosMultiphysics.STRUCTURE_VELOCITY) @@ -165,6 +165,10 @@ def ImportModelPart(self) -> None: def PrepareModelPart(self) -> None: self.model_part.SetBufferSize(self.GetMinimumBufferSize()) + # Set ProcessInfo variables + self.model_part.ProcessInfo.SetValue(KratosMultiphysics.TIME, 0) + self.model_part.ProcessInfo.SetValue(KratosMultiphysics.STEP, 0) + def Check(self) -> None: super().Check() @@ -196,14 +200,12 @@ def Initialize(self) -> None: self.fluid_solver.ActivateWallResistance(self.wall_law_y) def InitializeSolutionStep(self) -> None: - # Transfer density and viscosity to the nodes + # Transfer density to the nodes for node in self.model_part.Nodes: if node.GetSolutionStepValue(KratosMultiphysics.DISTANCE) <= 0.0: node.SetSolutionStepValue(KratosMultiphysics.DENSITY, self.density) - node.SetSolutionStepValue(KratosMultiphysics.VISCOSITY, self.viscosity) else: node.SetSolutionStepValue(KratosMultiphysics.DENSITY, 0.0) - node.SetSolutionStepValue(KratosMultiphysics.VISCOSITY, 0.0) def AdvanceInTime(self, current_time: float) -> float: """ diff --git a/applications/FreeSurfaceApplication/python_scripts/free_surface_analysis.py b/applications/FreeSurfaceApplication/python_scripts/free_surface_analysis.py index 6c98ea2f7f45..59846302dae3 100644 --- a/applications/FreeSurfaceApplication/python_scripts/free_surface_analysis.py +++ b/applications/FreeSurfaceApplication/python_scripts/free_surface_analysis.py @@ -37,15 +37,14 @@ def Initialize(self): break for node in self.model.GetModelPart(model_part_name).Nodes: - # Initialize DENSITY - node.SetSolutionStepValue(KratosMultiphysics.DENSITY, density) - - # Initialize DISTANCE + # Initialize DISTANCE and DENSITY if node.GetSolutionStepValue(KratosMultiphysics.DISTANCE) < 0.0: active_node_count += 1 node.SetSolutionStepValue(KratosMultiphysics.DISTANCE, -small_value) + node.SetSolutionStepValue(KratosMultiphysics.DENSITY, density) else: node.SetSolutionStepValue(KratosMultiphysics.DISTANCE, small_value) + node.SetSolutionStepValue(KratosMultiphysics.DENSITY, 0.0) # Make sure no node has null porosity and diameter if node.GetSolutionStepValue(KratosMultiphysics.POROSITY) == 0.0: