diff --git a/SPlisHSPlasH/CMakeLists.txt b/SPlisHSPlasH/CMakeLists.txt index b0238866..96bf8726 100644 --- a/SPlisHSPlasH/CMakeLists.txt +++ b/SPlisHSPlasH/CMakeLists.txt @@ -75,6 +75,7 @@ set(SURFACETENSION_HEADER_FILES SurfaceTension/SurfaceTension_Becker2007.h SurfaceTension/SurfaceTension_Akinci2013.h SurfaceTension/SurfaceTension_He2014.h + SurfaceTension/SurfaceTension_Jeske2023.h ) set(SURFACETENSION_SOURCE_FILES @@ -82,6 +83,7 @@ set(SURFACETENSION_SOURCE_FILES SurfaceTension/SurfaceTension_Becker2007.cpp SurfaceTension/SurfaceTension_Akinci2013.cpp SurfaceTension/SurfaceTension_He2014.cpp + SurfaceTension/SurfaceTension_Jeske2023.cpp ) set(VISCOSITY_HEADER_FILES diff --git a/SPlisHSPlasH/NonPressureForceRegistration.cpp b/SPlisHSPlasH/NonPressureForceRegistration.cpp index 6ccccf6a..3a90bc75 100644 --- a/SPlisHSPlasH/NonPressureForceRegistration.cpp +++ b/SPlisHSPlasH/NonPressureForceRegistration.cpp @@ -20,6 +20,7 @@ #include "SurfaceTension/SurfaceTension_Becker2007.h" #include "SurfaceTension/SurfaceTension_Akinci2013.h" #include "SurfaceTension/SurfaceTension_He2014.h" +#include "SurfaceTension/SurfaceTension_Jeske2023.h" #ifdef USE_THIRD_PARTY_METHODS #include "SurfaceTension/SurfaceTension_ZorillaRitter2020.h" #endif @@ -42,6 +43,7 @@ void Simulation::registerNonpressureForces() addSurfaceTensionMethod("Becker & Teschner 2007", SurfaceTension_Becker2007::creator); addSurfaceTensionMethod("Akinci et al. 2013", SurfaceTension_Akinci2013::creator); addSurfaceTensionMethod("He et al. 2014", SurfaceTension_He2014::creator); + addSurfaceTensionMethod("Jeske et al. 2023", SurfaceTension_Jeske2023::creator); #ifdef USE_THIRD_PARTY_METHODS addSurfaceTensionMethod("Zorilla, Ritter, et al. 2020", SurfaceTension_ZorillaRitter2020::creator); #endif diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Jeske2023.cpp b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Jeske2023.cpp new file mode 100644 index 00000000..53432d88 --- /dev/null +++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Jeske2023.cpp @@ -0,0 +1,1329 @@ +#include "SurfaceTension_Jeske2023.h" +#include "SPlisHSPlasH/TimeManager.h" +#include "Utilities/Timing.h" +#include "Utilities/Counting.h" +#include "SPlisHSPlasH/Simulation.h" +#include "SPlisHSPlasH/BoundaryModel_Akinci2012.h" +#include "SPlisHSPlasH/BoundaryModel_Koschier2017.h" +#include "SPlisHSPlasH/BoundaryModel_Bender2019.h" +#include "SPlisHSPlasH/Utilities/MathFunctions.h" + +using namespace SPH; +using namespace GenParam; + +int SurfaceTension_Jeske2023::ITERATIONS = -1; +int SurfaceTension_Jeske2023::MAX_ITERATIONS = -1; +int SurfaceTension_Jeske2023::MAX_ERROR = -1; +int SurfaceTension_Jeske2023::XSPH = -1; + +int SurfaceTension_Jeske2023::VISCOSITY_COEFFICIENT = -1; +int SurfaceTension_Jeske2023::VISCOSITY_COEFFICIENT_BOUNDARY = -1; + +SurfaceTension_Jeske2023::SurfaceTension_Jeske2023(FluidModel *model) : + SurfaceTensionBase(model), m_vDiff(), m_gradRho() { + m_viscosity = 0; + m_boundaryViscosity = 0; + + m_maxIter = 100; + m_maxError = static_cast(0.001); // This seems to be a better default than 1e-2 avoiding spinning problems etc + m_iterations = 0; + m_xsph = 0; + + m_vDiff.resize(model->numParticles(), Vector3r::Zero()); + m_gradRho.resize(model->numParticles(), static_cast(0.0)); + m_surfaceEnergy.resize(model->numParticles(), static_cast(0.0)); + m_color.resize(model->numParticles(), static_cast(0.0)); + m_colorGrad.resize(model->numParticles(), Vector3r::Zero()); + m_nonlinearAcc.resize(model->numParticles(), Vector3r::Zero()); + m_nonlinearRes.resize(model->numParticles(), Vector3r::Zero()); + m_nonlinearGrad.resize(model->numParticles(), Vector3r::Zero()); + + + model->addField( {"velocity difference", FieldType::Vector3, [&](const unsigned int i) -> Real * { return &m_vDiff[i][0]; }, true}); + model->addField( {"density gradient", FieldType::Scalar, [&](const unsigned int i) -> Real * { return &m_gradRho[i]; }, true}); + model->addField( {"surface energy", FieldType::Scalar, [&](const unsigned int i) -> Real * { return &m_surfaceEnergy[i]; }, true}); + model->addField( {"surface color", FieldType::Scalar, [&](const unsigned int i) -> Real * { return &m_color[i]; }, true}); + model->addField( {"surface color grad", FieldType::Vector3, [&](const unsigned int i) -> Real * { return &m_colorGrad[i][0]; }, true}); + model->addField( {"surface nonlinear acc", FieldType::Vector3, [&](const unsigned int i) -> Real * { return &m_nonlinearAcc[i][0]; }, true}); + model->addField( {"surface nonlinear res", FieldType::Vector3, [&](const unsigned int i) -> Real * { return &m_nonlinearRes[i][0]; }, true}); + model->addField( {"surface nonlinear grad", FieldType::Vector3, [&](const unsigned int i) -> Real * { return &m_nonlinearGrad[i][0]; }, true}); +} + +SurfaceTension_Jeske2023::~SurfaceTension_Jeske2023(void) { + m_model->removeFieldByName("velocity difference"); + m_model->removeFieldByName("density gradient"); + m_model->removeFieldByName("surface energy"); + m_model->removeFieldByName("surface color"); + m_model->removeFieldByName("surface color grad"); + m_model->removeFieldByName("surface nonlinear acc"); + m_model->removeFieldByName("surface nonlinear res"); + m_model->removeFieldByName("surface nonlinear grad"); + + m_vDiff.clear(); +} + +void SurfaceTension_Jeske2023::initParameters() { + SurfaceTensionBase::initParameters(); + + VISCOSITY_COEFFICIENT = createNumericParameter("surfaceTensionViscosity", "Viscosity coefficient ", &m_viscosity); + setGroup(VISCOSITY_COEFFICIENT, "Fluid Model|Surface tension"); + setDescription(VISCOSITY_COEFFICIENT, "Coefficient for the viscosity force computation."); + RealParameter* rparam = static_cast(getParameter(VISCOSITY_COEFFICIENT)); + rparam->setMinValue(0.0); + + VISCOSITY_COEFFICIENT_BOUNDARY = createNumericParameter("surfaceTensionViscosityBoundary", "Viscosity coefficient (Boundary)", &m_boundaryViscosity); + setGroup(VISCOSITY_COEFFICIENT_BOUNDARY, "Fluid Model|Surface tension"); + setDescription(VISCOSITY_COEFFICIENT_BOUNDARY, "Coefficient for the viscosity force computation at the boundary."); + rparam = static_cast(getParameter(VISCOSITY_COEFFICIENT_BOUNDARY)); + rparam->setMinValue(0.0); + + ITERATIONS = createNumericParameter("surfaceTensionIterations", "Iterations", &m_iterations); + setGroup(ITERATIONS, "Fluid Model|Surface tension"); + setDescription(ITERATIONS, "Iterations required by the viscosity solver."); + getParameter(ITERATIONS)->setReadOnly(true); + + MAX_ITERATIONS = createNumericParameter("surfaceTensionMaxIter", "Max. iterations (visco)", &m_maxIter); + setGroup(MAX_ITERATIONS, "Fluid Model|Surface tension"); + setDescription(MAX_ITERATIONS, "Max. iterations of the viscosity solver."); + static_cast *>(getParameter(MAX_ITERATIONS))->setMinValue(1); + + MAX_ERROR = createNumericParameter("surfaceTensionMaxError", "Max. surface tension error", &m_maxError); + setGroup(MAX_ERROR, "Fluid Model|Surface tension"); + setDescription(MAX_ERROR, "Max. error of the surface tension solver."); + // rparam = static_cast(getParameter(MAX_ERROR)); + // rparam->setMinValue(static_cast(1e-6)); + + XSPH = createNumericParameter("surfaceTensionXSPH", "XSPH Smoothing Factor", &m_xsph); + setGroup(XSPH, "Fluid Model|Surface tension"); + setDescription(XSPH, "Factor for xsph smoothing"); +} + +#ifdef USE_AVX + +void SurfaceTension_Jeske2023::step() { + const int numParticles = (int)m_model->numActiveParticles(); + // prevent solver from running with a zero-length vector + if (numParticles == 0) + return; + const Real density0 = m_model->getDensity0(); + const Real h = TimeManager::getCurrent()->getTimeStepSize(); + + // Compute density gradient first + computeDensityGradient(); + + ////////////////////////////////////////////////////////////////////////// + // Init linear system solver and preconditioner + ////////////////////////////////////////////////////////////////////////// + MatrixReplacement A(3 * m_model->numActiveParticles(), matrixVecProd, (void*)this); + + m_solver.setTolerance(m_maxError); + m_solver.setMaxIterations(m_maxIter); + m_solver.compute(A); + + VectorXr b(3 * numParticles); + VectorXr x(3 * numParticles); + VectorXr g(3 * numParticles); + + computeRHS(b, g); + + ////////////////////////////////////////////////////////////////////////// + // Solve linear system + ////////////////////////////////////////////////////////////////////////// + START_TIMING("CG solve"); + x = m_solver.solveWithGuess(b, g); + m_iterations = (int)m_solver.iterations(); + STOP_TIMING_AVG; + INCREASE_COUNTER("Surface tension iterations", static_cast(m_iterations)); + + applyForces(x); +} + +void SurfaceTension_Jeske2023::matrixVecProd(const Real *vec, Real *result, void *userData) { + Simulation *sim = Simulation::getCurrent(); + SurfaceTension_Jeske2023 *surf_tens = (SurfaceTension_Jeske2023 *) userData; + FluidModel *model = surf_tens->getModel(); + const unsigned int numParticles = model->numActiveParticles(); + const unsigned int fluidModelIndex = model->getPointSetIndex(); + const unsigned int nFluids = sim->numberOfFluidModels(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + + const Real h = sim->getSupportRadius(); + const Real h2 = h * h; + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + const Real density0 = model->getDensity0(); + const Real mu = surf_tens->m_viscosity * density0; + const Real mub = surf_tens->m_boundaryViscosity * density0; + const Real mu_xsph = surf_tens->m_xsph; + const Real sphereVolume = static_cast(4.0 / 3.0 * M_PI) * h2 * h; + + const Real cohesion = surf_tens->m_surfaceTension; + const Real adhesion = surf_tens->m_surfaceTensionBoundary; + const Real radius = sim->getValue(Simulation::PARTICLE_RADIUS); + const Real diameter = static_cast(2.0) * radius; + const Real diameter2 = diameter * diameter; + const Real W_min = sim->W(Vector3r(diameter, 0, 0)); + + Real d = 10.0; + if (sim->is2DSimulation()) + d = 8.0; + + const Scalarf8 d_mu_rho0(d * mu * density0); + const Scalarf8 h2_001(0.01f * h2); + const Scalarf8 W_min_avx(W_min); + const Scalarf8 dt_avx(dt); + const Scalarf8 diameter2_avx(diameter * diameter); + const Scalarf8 density0_avx(density0); + const Scalarf8 dt_cohesion_avx(dt * cohesion); + const Scalarf8 mu_xsph_avx(mu_xsph); + Vector3f8 zero; + zero.setZero(); + +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) + for (int i = 0; i < (int) numParticles; i++) { + const Vector3r &xi = model->getPosition(i); + Vector3r ai; + ai.setZero(); + const Real density_i = model->getDensity(i); + const Vector3r &vi = Eigen::Map(&vec[3 * i]); + const Vector3r &vi_old = model->getVelocity(i); + // const Vector3r vi_old = Vector3r::Zero(); + const Real m_i = model->getMass(i); + + const Vector3f8 xi_avx(xi); + const Vector3f8 vi_avx(vi); + const Vector3f8 vi_old_avx(vi_old); + const Scalarf8 densityGrad_i_avx(surf_tens->getDensityGrad(i)); + Vector3f8 delta_ai_avx; + const Scalarf8 density_i_avx(density_i); + delta_ai_avx.setZero(); + + ////////////////////////////////////////////////////////////////////////// + // Fluid + ////////////////////////////////////////////////////////////////////////// + forall_fluid_neighbors_in_same_phase_avx( + compute_Vj(model); + compute_Vj_gradW_samephase(); + const Scalarf8 density_j_avx = convert_one(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &model->getDensity(0), count); + const Scalarf8 densityGrad_j_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &surf_tens->getDensityGrad(0), count); + const Scalarf8 Vj_avx = convert_zero(model->getVolume(0), count); + const Scalarf8 mass_j_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &model->getMass(0), count); + const Vector3f8 xixj = xi_avx - xj_avx; + const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &vec[0], count); + const Vector3f8 vj_old_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &model->getVelocity(0), count); + + delta_ai_avx += V_gradW * ((d_mu_rho0 / density_j_avx) * (vi_avx - vj_avx).dot(xixj) / (xixj.squaredNorm() + h2_001)); + delta_ai_avx -= (vi_avx - vj_avx) * mu_xsph_avx * mass_j_avx/density_j_avx * CubicKernel_AVX::W(xixj) * density_i_avx / dt_avx; + + const Scalarf8 W_cohesion = Scalarf8(10.f/7) * min(CubicKernel_AVX::W(xixj), W_min_avx); + const Vector3f8 gradW = CubicKernel_AVX::gradW(xixj); + const Vector3f8 gradW_cohesion = Vector3f8::blend((xixj).squaredNorm() > diameter2_avx, gradW, zero) * Scalarf8(10.f/7) ; + const Scalarf8 rhoij_avx = density_i_avx + density_j_avx; + + const Scalarf8 W_firstOrder = convert_zero(2.0, count) / rhoij_avx * (W_cohesion + dt_avx * (gradW_cohesion.dot(vi_old_avx - vj_old_avx) + - W_cohesion * (densityGrad_i_avx + densityGrad_j_avx) / rhoij_avx) + ); + + delta_ai_avx -= (vi_avx - vj_avx) * dt_cohesion_avx * density_i_avx * W_firstOrder; + + ); + + ai[0] += delta_ai_avx.x().reduce(); + ai[1] += delta_ai_avx.y().reduce(); + ai[2] += delta_ai_avx.z().reduce(); + + ////////////////////////////////////////////////////////////////////////// + // Boundary + ////////////////////////////////////////////////////////////////////////// + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + // bm_neighbor->addForce(xj, -model->getMass(i) / density_i * a); + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + // ToDo: implement adhesion for moving boundaries, i.e. vj not equal to zero + forall_volume_maps( + Vector3r vj_old; + bm_neighbor->getPointVelocity(xj, vj_old); + const Vector3r xixj = xi - xj; + //const Real rho_iT1 = surf_tens->m_rhoT1[i]; + const Real m_j = Vj * density0; + const Real density_j = m_j / Vj; + //const Real Vij = (m_i/rho_iT1 + Vj)*0.5; + const Real W_cohesion = std::min(sim->W(xixj), W_min); + const Vector3r gradW = sim->gradW(xixj); + const Vector3r gradW_cohesion = (xixj).norm() > diameter ? gradW : Vector3r::Zero(); + + // First order + // ai -= (vi * dt - vj * dt) * 2/(density_i + density_j) * cohesion * density_i * (std::min(sim->W(xixj), W_min) + dt * gradW_cohesion.dot(vi_old) - dt*gradW_cohesion.dot(vj_old)); + + const Real rhoij = density_i + density_j; + const Real m_ij = (m_i + m_j)/m_i; + const Real W_firstOrder = m_ij * (W_cohesion / rhoij + dt * gradW_cohesion.dot(vi_old) / rhoij + - dt * gradW_cohesion.dot(vj_old) / rhoij + - dt * W_cohesion * (surf_tens->getDensityGrad(i) + // + surf_tens->getDensityGrad(neighborIndex) + ) / (rhoij * rhoij) + ); + ai -= (vi * dt + // - vj * dt + ) * adhesion * density_i * W_firstOrder; + //bm_neighbor->addForce(xi, -model->getMass(i) * ai); + ); + } + + if (mub != 0.0) + { + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + const Vector3r a = d * mub * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vi).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW; + ai += a; + // bm_neighbor->addForce(xj, -model->getMass(i) / density_i * a); + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = surf_tens->m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1 * dist; + const Vector3r x2 = xj + t1 * dist; + const Vector3r x3 = xj - t2 * dist; + const Vector3r x4 = xj + t2 * dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * rho * sphereVolume; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (vi).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (vi).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (vi).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (vi).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + ai += a1 + a2 + a3 + a4; + + //bm_neighbor->addForce(x1, -model->getMass(i)/density_i * a1); + //bm_neighbor->addForce(x2, -model->getMass(i)/density_i * a2); + //bm_neighbor->addForce(x3, -model->getMass(i)/density_i * a3); + //bm_neighbor->addForce(x4, -model->getMass(i)/density_i * a4); + } + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + forall_volume_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = surf_tens->m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1*dist; + const Vector3r x2 = xj + t1*dist; + const Vector3r x3 = xj - t2*dist; + const Vector3r x4 = xj + t2*dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * Vj; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (vi).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (vi).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (vi).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (vi).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + ai += a1 + a2 + a3 + a4; + + //bm_neighbor->addForce(x1, -model->getMass(i)/density_i * a1); + //bm_neighbor->addForce(x2, -model->getMass(i)/density_i * a2); + //bm_neighbor->addForce(x3, -model->getMass(i)/density_i * a3); + //bm_neighbor->addForce(x4, -model->getMass(i)/density_i * a4); + } + ); + } + } + + result[3 * i] = vec[3 * i] - dt / density_i * ai[0]; + result[3 * i + 1] = vec[3 * i + 1] - dt / density_i * ai[1]; + result[3 * i + 2] = vec[3 * i + 2] - dt / density_i * ai[2]; + } + } +} + +void SurfaceTension_Jeske2023::computeDensityGradient() { + const int numParticles = (int) m_model->numActiveParticles(); + Simulation *sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + const unsigned int fluidModelIndex = m_model->getPointSetIndex(); + FluidModel *model = m_model; + const Real density0 = model->getDensity0(); + const Scalarf8 density0_avx(density0); + Vector3f8 zero; + zero.setZero(); + + ////////////////////////////////////////////////////////////////////////// + // Compute RHS + ////////////////////////////////////////////////////////////////////////// +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) nowait + for (int i = 0; i < (int) numParticles; i++) { + const Vector3r &xi = m_model->getPosition(i); + const Vector3r &vi = m_model->getVelocity(i); + Scalarf8 gradRho_avx(0.0); + const Vector3f8 xi_avx(xi); + const Vector3f8 vi_avx(vi); + + // Compute density gradient + forall_fluid_neighbors_in_same_phase_avx( + compute_Vj(model); + compute_Vj_gradW_samephase(); + + const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &model->getVelocity(0), count); + const Vector3f8 xixj = xi_avx - xj_avx; + + gradRho_avx += density0_avx * (vi_avx - vj_avx).dot(V_gradW); + ); + m_gradRho[i] = gradRho_avx.reduce(); + } + } +} + +void SurfaceTension_Jeske2023::computeRHS(VectorXr &b, VectorXr &g) { + const int numParticles = (int) m_model->numActiveParticles(); + Simulation *sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + const unsigned int fluidModelIndex = m_model->getPointSetIndex(); + FluidModel *model = m_model; + const Real h = sim->getSupportRadius(); + const Real h2 = h * h; + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + const Real density0 = m_model->getDensity0(); + const Real mu = m_viscosity * density0; + const Real mub = m_boundaryViscosity * density0; + const Real cohesion = m_surfaceTension; + const Real adhesion = m_surfaceTensionBoundary; + const Real sphereVolume = static_cast(4.0 / 3.0 * M_PI) * h2 * h; + const Real radius = sim->getValue(Simulation::PARTICLE_RADIUS); + const Real diameter = static_cast(2.0) * radius; + const Real diameter2 = diameter * diameter; + const Real W_min = sim->W(Vector3r(diameter, 0, 0)); + Real d = 10.0; + if (sim->is2DSimulation()) + d = 8.0; + + const Scalarf8 W_min_avx(W_min); + const Scalarf8 dt_avx(dt); + const Scalarf8 diameter2_avx(diameter * diameter); + const Scalarf8 density0_avx(density0); + const Scalarf8 cohesion_avx(cohesion); + Vector3f8 zero; + zero.setZero(); + + ////////////////////////////////////////////////////////////////////////// + // Compute RHS + ////////////////////////////////////////////////////////////////////////// +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) nowait + for (int i = 0; i < (int) numParticles; i++) { + const Vector3r &vi = m_model->getVelocity(i); + const Vector3r vi_old = vi; + // const Vector3r vi_old = Vector3r::Zero(); + const Vector3r &xi = m_model->getPosition(i); + const Real density_i = m_model->getDensity(i); + const Real m_i = m_model->getMass(i); + + const Vector3f8 xi_avx(xi); + const Vector3f8 vi_avx(vi); + Vector3f8 delta_bi_avx; + const Scalarf8 density_i_avx(density_i); + delta_bi_avx.setZero(); + const Scalarf8 densityGrad_i_avx(getDensityGrad(i)); + + // Implicit cohesion + forall_fluid_neighbors_in_same_phase_avx( + compute_Vj(model); + compute_Vj_gradW_samephase(); + const Vector3f8 xixj = xi_avx - xj_avx; + const Scalarf8 Vj_avx = convert_zero(model->getVolume(0), count); + const Scalarf8 density_j_avx = convert_one(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &model->getDensity(0), count); + const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &model->getVelocity(0), count); + const Scalarf8 densityGrad_j_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, fluidModelIndex, i)[j], &getDensityGrad(0), count); + + // Cohesion + const Scalarf8 W_cohesion = Scalarf8(10.f/7) * min(CubicKernel_AVX::W(xixj), W_min_avx); + const Vector3f8 gradW = CubicKernel_AVX::gradW(xixj); + const Vector3f8 gradW_cohesion = Vector3f8::blend((xixj).squaredNorm() > diameter2_avx, gradW, zero) * Scalarf8(10.f/7); + + // First order + const Scalarf8 rhoij_avx = density_i_avx + density_j_avx; + const Scalarf8 W_firstOrder = convert_zero(2.0, count) / rhoij_avx * (W_cohesion + dt_avx * (gradW_cohesion.dot(vi_avx) + - gradW_cohesion.dot(vj_avx) + - W_cohesion * (densityGrad_i_avx + densityGrad_j_avx) / rhoij_avx) + ); + + delta_bi_avx += xixj * cohesion_avx * density_i_avx * W_firstOrder; + ); + Vector3r bi = delta_bi_avx.reduce(); + + ////////////////////////////////////////////////////////////////////////// + // Boundary + ////////////////////////////////////////////////////////////////////////// + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + // bm_neighbor->addForce(xj, -model->getMass(i) / density_i * a); + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + // ToDo: implement adhesion for moving boundaries, i.e. vj not equal to zero + forall_volume_maps( + Vector3r vj; + bm_neighbor->getPointVelocity(xj, vj); + const Vector3r xixj = xi - xj; + //const Real rho_iT1 = m_rhoT1[i]; + const Real m_j = Vj * density0; + const Real density_j = m_j / Vj; + //const Real Vij = (m_i/rho_iT1 + Vj)*0.5; + const Real W_cohesion = std::min(sim->W(xixj), W_min); + const Vector3r gradW = sim->gradW(xixj); + const Vector3r gradW_cohesion = (xixj).norm() > diameter ? gradW : Vector3r::Zero(); + + const Real m_ij = (m_i + m_j)/m_i; + const Real rhoij = density_i + density_j; + // TODO: some issue with adhesion here + const Real W_firstOrder = m_ij * ( W_cohesion / rhoij + + dt * gradW_cohesion.dot(vi) / rhoij + - dt * gradW_cohesion.dot(vj) / rhoij + - dt * W_cohesion * (getDensityGrad(i) + // + getDensityGrad(neighborIndex) + ) / (rhoij * rhoij) + ); + const Vector3r a = adhesion * xixj * density_i * W_firstOrder; + bi += a; + bm_neighbor->addForce(xi, -m_model->getMass(i) / density_i * a); + ); + } + if (mub != 0.0) + { + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + const Vector3r a = d * mub * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW; + bi += a; + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1 * dist; + const Vector3r x2 = xj + t1 * dist; + const Vector3r x3 = xj - t2 * dist; + const Vector3r x4 = xj + t2 * dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * rho * sphereVolume; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (v1).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (v2).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (v3).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (v4).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + bi += a1 + a2 + a3 + a4; + } + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + forall_volume_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1*dist; + const Vector3r x2 = xj + t1*dist; + const Vector3r x3 = xj - t2*dist; + const Vector3r x4 = xj + t2*dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * Vj; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (v1).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (v2).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (v3).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (v4).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + bi += a1 + a2 + a3 + a4; + } + ); + } + } + + b[3 * i] = vi[0] - dt / density_i * bi[0]; + b[3 * i + 1] = vi[1] - dt / density_i * bi[1]; + b[3 * i + 2] = vi[2] - dt / density_i * bi[2]; + + // Warmstart + g[3 * i] = vi[0] + m_vDiff[i][0]; + g[3 * i + 1] = vi[1] + m_vDiff[i][1]; + g[3 * i + 2] = vi[2] + m_vDiff[i][2]; + } + } +} + + +#else + +void SurfaceTension_Jeske2023::matrixVecProd(const Real *vec, Real *result, void *userData) { + Simulation *sim = Simulation::getCurrent(); + SurfaceTension_Jeske2023 *surf_tens = (SurfaceTension_Jeske2023 *) userData; + FluidModel *model = surf_tens->getModel(); + const unsigned int numParticles = model->numActiveParticles(); + const unsigned int fluidModelIndex = model->getPointSetIndex(); + const unsigned int nFluids = sim->numberOfFluidModels(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + + const Real h = sim->getSupportRadius(); + const Real h2 = h * h; + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + const Real density0 = model->getDensity0(); + const Real mu = surf_tens->m_viscosity * density0; + const Real mub = surf_tens->m_boundaryViscosity * density0; + const Real sphereVolume = static_cast(4.0 / 3.0 * M_PI) * h2 * h; + + const Real cohesion = surf_tens->m_surfaceTension; + const Real adhesion = surf_tens->m_surfaceTensionBoundary; + const Real mu_xsph = surf_tens->m_xsph; + const Real radius = sim->getValue(Simulation::PARTICLE_RADIUS); + const Real diameter = static_cast(2.0) * radius; + const Real diameter2 = diameter * diameter; + const Real W_min = sim->W(Vector3r(diameter, 0, 0)); + + Real d = 10.0; + if (sim->is2DSimulation()) + d = 8.0; + +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) + for (int i = 0; i < (int) numParticles; i++) { + const Vector3r &xi = model->getPosition(i); + Vector3r ai; + ai.setZero(); + const Real density_i = model->getDensity(i); + const Vector3r &vi = Eigen::Map(&vec[3 * i]); + const Vector3r &vi_old = model->getVelocity(i); + // const Vector3r vi_old = Vector3r::Zero(); + const Real m_i = model->getMass(i); + + ////////////////////////////////////////////////////////////////////////// + // Fluid + ////////////////////////////////////////////////////////////////////////// + forall_fluid_neighbors_in_same_phase( + const Real density_j = model->getDensity(neighborIndex); + const Real mass_j = model->getMass(neighborIndex); + const Vector3r gradW = sim->gradW(xi - xj); + + const Vector3r &vj = Eigen::Map(&vec[3 * neighborIndex]); + const Vector3r &vj_old = model->getVelocity(neighborIndex); + // const Vector3r vj_old = Vector3r::Zero(); + const Vector3r xixj = xi - xj; + const Real m_ij = m_i + mass_j; + + ai += d * mu * (model->getMass(neighborIndex) / density_j) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW; + ai -= (vi - vj) * mu_xsph * mass_j / density_j * CubicKernel::W(xixj) * density_i / dt; + + // Cohesion + const Real W_cohesion = (10./7.) * std::min(sim->W(xixj), W_min); + Vector3r gradW_cohesion = (xi - xj).norm() > diameter ? gradW : Vector3r::Zero(); + gradW_cohesion *= (10./7.); + + // First order + // ai -= (vi * dt - vj * dt) * 2/(density_i + density_j) * cohesion * density_i * (std::min(sim->W(xixj), W_min) + dt * gradW_cohesion.dot(vi_old) - dt*gradW_cohesion.dot(vj_old)); + + const Real rhoij = density_i + density_j; + const Real W_firstOrder = 2.0 * (W_cohesion / rhoij + dt * gradW_cohesion.dot(vi_old) / rhoij // - dt * W_cohesion * (surf_tens->getDensityGrad(i) - m_i * gradW ).dot(vi_old) / (rhoij * rhoij) + - dt * gradW_cohesion.dot(vj_old) / rhoij // - dt * W_cohesion * (surf_tens->getDensityGrad(neighborIndex) - mass_j * gradW ).dot(vj_old) / (rhoij * rhoij) + - dt * W_cohesion * (surf_tens->getDensityGrad(i) + surf_tens->getDensityGrad(neighborIndex)) / (rhoij * rhoij) + ); + ai -= (vi * dt - vj * dt) * cohesion * density_i * W_firstOrder; + + ); + + ////////////////////////////////////////////////////////////////////////// + // Boundary + ////////////////////////////////////////////////////////////////////////// + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + // bm_neighbor->addForce(xj, -model->getMass(i) / density_i * a); + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + // ToDo: implement adhesion for moving boundaries, i.e. vj not equal to zero + forall_volume_maps( + Vector3r vj_old; + bm_neighbor->getPointVelocity(xj, vj_old); + const Vector3r xixj = xi - xj; + const Real m_j = Vj * density0; + const Real density_j = m_j / Vj; + const Real W_cohesion = std::min(sim->W(xixj), W_min); + const Vector3r gradW = sim->gradW(xixj); + const Vector3r gradW_cohesion = (xixj).norm() > diameter ? gradW : Vector3r::Zero(); + + const Real rhoij = density_i + density_j; + const Real m_ij = (m_i + m_j)/m_i; + const Real W_firstOrder = m_ij * (W_cohesion / rhoij + dt * gradW_cohesion.dot(vi_old) / rhoij + - dt * gradW_cohesion.dot(vj_old) / rhoij + - dt * W_cohesion * (surf_tens->getDensityGrad(i) + // + surf_tens->getDensityGrad(neighborIndex) + ) / (rhoij * rhoij) + ); + ai -= (vi * dt + ) * adhesion * density_i * W_firstOrder; + ); + } + + if (mub != 0.0) + { + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + const Vector3r a = d * mub * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vi).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW; + ai += a; + // bm_neighbor->addForce(xj, -model->getMass(i) / density_i * a); + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = surf_tens->m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1 * dist; + const Vector3r x2 = xj + t1 * dist; + const Vector3r x3 = xj - t2 * dist; + const Vector3r x4 = xj + t2 * dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * rho * sphereVolume; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (vi).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (vi).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (vi).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (vi).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + ai += a1 + a2 + a3 + a4; + + // bm_neighbor->addForce(x1, -model->getMass(i)/density_i * a1); + // bm_neighbor->addForce(x2, -model->getMass(i)/density_i * a2); + // bm_neighbor->addForce(x3, -model->getMass(i)/density_i * a3); + // bm_neighbor->addForce(x4, -model->getMass(i)/density_i * a4); + } + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + forall_volume_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = surf_tens->m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1*dist; + const Vector3r x2 = xj + t1*dist; + const Vector3r x3 = xj - t2*dist; + const Vector3r x4 = xj + t2*dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * Vj; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (vi).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (vi).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (vi).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (vi).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + ai += a1 + a2 + a3 + a4; + + // bm_neighbor->addForce(x1, -model->getMass(i)/density_i * a1); + // bm_neighbor->addForce(x2, -model->getMass(i)/density_i * a2); + // bm_neighbor->addForce(x3, -model->getMass(i)/density_i * a3); + // bm_neighbor->addForce(x4, -model->getMass(i)/density_i * a4); + } + ); + } + } + + result[3 * i] = vec[3 * i] - dt / density_i * ai[0]; + result[3 * i + 1] = vec[3 * i + 1] - dt / density_i * ai[1]; + result[3 * i + 2] = vec[3 * i + 2] - dt / density_i * ai[2]; + } + } +} + +void SurfaceTension_Jeske2023::step() { + const int numParticles = (int)m_model->numActiveParticles(); + // prevent solver from running with a zero-length vector + if (numParticles == 0) + return; + const Real density0 = m_model->getDensity0(); + const Real h = TimeManager::getCurrent()->getTimeStepSize(); + + // Compute xsph first + // step_xsph(); +// solve_nonlinear(); + + // Compute density gradient first + computeDensityGradient(); + + ////////////////////////////////////////////////////////////////////////// + // Init linear system solver and preconditioner + ////////////////////////////////////////////////////////////////////////// + MatrixReplacement A(3 * m_model->numActiveParticles(), matrixVecProd, (void*)this); + // m_solver.preconditioner().init(m_model->numActiveParticles(), diagonalMatrixElement, (void*)this); + + m_solver.setTolerance(m_maxError); + m_solver.setMaxIterations(m_maxIter); + m_solver.compute(A); + + VectorXr b(3 * numParticles); + VectorXr x(3 * numParticles); + VectorXr g(3 * numParticles); + + computeRHS(b, g); + + ////////////////////////////////////////////////////////////////////////// + // Solve linear system + ////////////////////////////////////////////////////////////////////////// + START_TIMING("CG solve"); + x = m_solver.solveWithGuess(b, g); + m_iterations = (int)m_solver.iterations(); + STOP_TIMING_AVG; + INCREASE_COUNTER("Surface tension iterations", static_cast(m_iterations)); + + applyForces(x); + +} + +void SurfaceTension_Jeske2023::computeDensityGradient() { + const int numParticles = (int) m_model->numActiveParticles(); + Simulation *sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + const unsigned int fluidModelIndex = m_model->getPointSetIndex(); + FluidModel *model = m_model; + const Real radius = sim->getValue(Simulation::PARTICLE_RADIUS); + const Real diameter = static_cast(2.0) * radius; + const Real diameter2 = diameter * diameter; + // Reset density gradients + std::fill(m_gradRho.begin(), m_gradRho.end(), static_cast(0.0)); + + ////////////////////////////////////////////////////////////////////////// + // Compute RHS + ////////////////////////////////////////////////////////////////////////// +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) nowait + for (int i = 0; i < (int) numParticles; i++) { + const Vector3r &xi = m_model->getPosition(i); + const Vector3r &vi = m_model->getVelocity(i); + + // Compute density gradient + forall_fluid_neighbors_in_same_phase( + const Vector3r xixj = xi - xj; + const Vector3r &vj = m_model->getVelocity(neighborIndex); + const Vector3r gradW = sim->gradW(xixj); + const Real m_j = model->getMass(neighborIndex); + m_gradRho[i] += m_j * (vi - vj).dot(gradW); + ) + } + } +} + +void SurfaceTension_Jeske2023::computeRHS(VectorXr &b, VectorXr &g) { + const int numParticles = (int) m_model->numActiveParticles(); + Simulation *sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + const unsigned int fluidModelIndex = m_model->getPointSetIndex(); + FluidModel *model = m_model; + const Real h = sim->getSupportRadius(); + const Real h2 = h * h; + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + const Real density0 = m_model->getDensity0(); + const Real mu = m_viscosity * density0; + const Real mub = m_boundaryViscosity * density0; + const Real cohesion = m_surfaceTension; + const Real adhesion = m_surfaceTensionBoundary; + const Real sphereVolume = static_cast(4.0 / 3.0 * M_PI) * h2 * h; + const Real radius = sim->getValue(Simulation::PARTICLE_RADIUS); + const Real diameter = static_cast(2.0) * radius; + const Real diameter2 = diameter * diameter; + const Real W_min = sim->W(Vector3r(diameter, 0, 0)); + Real d = 10.0; + if (sim->is2DSimulation()) + d = 8.0; + + ////////////////////////////////////////////////////////////////////////// + // Compute RHS + ////////////////////////////////////////////////////////////////////////// +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) nowait + for (int i = 0; i < (int) numParticles; i++) { + const Vector3r &vi = m_model->getVelocity(i); + const Vector3r vi_old = vi; + const Vector3r &xi = m_model->getPosition(i); + const Real density_i = m_model->getDensity(i); + const Real m_i = m_model->getMass(i); + Vector3r bi = Vector3r::Zero(); + + // Implicit cohesion + forall_fluid_neighbors_in_same_phase( + const Vector3r xixj = xi - xj; + const Real m_j = model->getMass(neighborIndex); + const Real m_ij = m_i + m_j; + const Real density_j = m_model->getDensity(neighborIndex); + const Vector3r &vj = model->getVelocity(neighborIndex); + const Vector3r vj_old = vj; + const Vector3r gradW = sim->gradW(xi - xj); + + // Cohesion + const Real W_cohesion = (10./7.) * std::min(sim->W(xixj), W_min); + Vector3r gradW_cohesion = (xi - xj).norm() > diameter ? gradW : Vector3r::Zero(); + gradW_cohesion *= (10./7.); + + const Real rhoij = density_i + density_j; + const Real W_firstOrder = 2.0 * ( W_cohesion / rhoij + + dt * gradW_cohesion.dot(vi) / rhoij + - dt * gradW_cohesion.dot(vj) / rhoij + - dt * W_cohesion * (getDensityGrad(i) + getDensityGrad(neighborIndex)) / (rhoij * rhoij) + ); + bi += cohesion * xixj * density_i * W_firstOrder; + + ) + + ////////////////////////////////////////////////////////////////////////// + // Boundary + ////////////////////////////////////////////////////////////////////////// + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + // bm_neighbor->addForce(xj, -model->getMass(i) / density_i * a); + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + // ToDo: implement adhesion for moving boundaries, i.e. vj not equal to zero + forall_volume_maps( + Vector3r vj; + bm_neighbor->getPointVelocity(xj, vj); + const Vector3r xixj = xi - xj; + //const Real rho_iT1 = m_rhoT1[i]; + const Real m_j = Vj * density0; + const Real density_j = m_j / Vj; + //const Real Vij = (m_i/rho_iT1 + Vj)*0.5; + const Real W_cohesion = std::min(sim->W(xixj), W_min); + const Vector3r gradW = sim->gradW(xixj); + const Vector3r gradW_cohesion = (xixj).norm() > diameter ? gradW : Vector3r::Zero(); + + const Real m_ij = (m_i + m_j)/m_i; + const Real rhoij = density_i + density_j; + // TODO: some issue with adhesion here + const Real W_firstOrder = m_ij * ( W_cohesion / rhoij + + dt * gradW_cohesion.dot(vi) / rhoij + - dt * gradW_cohesion.dot(vj) / rhoij + - dt * W_cohesion * (getDensityGrad(i) + // + getDensityGrad(neighborIndex) + ) / (rhoij * rhoij) + ); + bi += adhesion * xixj * density_i * W_firstOrder; + + ); + } + if (mub != 0.0) + { + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + { + forall_boundary_neighbors( + const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); + const Vector3r xixj = xi - xj; + const Vector3r gradW = sim->gradW(xixj); + const Vector3r a = d * mub * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW; + bi += a; + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + { + forall_density_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1 * dist; + const Vector3r x2 = xj + t1 * dist; + const Vector3r x3 = xj - t2 * dist; + const Vector3r x4 = xj + t2 * dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * rho * sphereVolume; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (v1).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (v2).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (v3).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (v4).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + bi += a1 + a2 + a3 + a4; + } + ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + { + forall_volume_maps( + const Vector3r xixj = xi - xj; + Vector3r normal = -xixj; + const Real nl = normal.norm(); + if (nl > static_cast(0.0001)) + { + normal /= nl; + + Vector3r t1; + Vector3r t2; + MathFunctions::getOrthogonalVectors(normal, t1, t2); + + const Real dist = m_tangentialDistanceFactor * h; + const Vector3r x1 = xj - t1*dist; + const Vector3r x2 = xj + t1*dist; + const Vector3r x3 = xj - t2*dist; + const Vector3r x4 = xj + t2*dist; + + const Vector3r xix1 = xi - x1; + const Vector3r xix2 = xi - x2; + const Vector3r xix3 = xi - x3; + const Vector3r xix4 = xi - x4; + + const Vector3r gradW1 = sim->gradW(xix1); + const Vector3r gradW2 = sim->gradW(xix2); + const Vector3r gradW3 = sim->gradW(xix3); + const Vector3r gradW4 = sim->gradW(xix4); + + // each sample point represents the quarter of the volume inside of the boundary + const Real vol = static_cast(0.25) * Vj; + + Vector3r v1; + Vector3r v2; + Vector3r v3; + Vector3r v4; + bm_neighbor->getPointVelocity(x1, v1); + bm_neighbor->getPointVelocity(x2, v2); + bm_neighbor->getPointVelocity(x3, v3); + bm_neighbor->getPointVelocity(x4, v4); + + // compute forces for both sample point + const Vector3r a1 = d * mub * vol * (v1).dot(xix1) / (xix1.squaredNorm() + 0.01*h2) * gradW1; + const Vector3r a2 = d * mub * vol * (v2).dot(xix2) / (xix2.squaredNorm() + 0.01*h2) * gradW2; + const Vector3r a3 = d * mub * vol * (v3).dot(xix3) / (xix3.squaredNorm() + 0.01*h2) * gradW3; + const Vector3r a4 = d * mub * vol * (v4).dot(xix4) / (xix4.squaredNorm() + 0.01*h2) * gradW4; + bi += a1 + a2 + a3 + a4; + } + ); + } + } + + b[3 * i] = vi[0] - dt / density_i * bi[0]; + b[3 * i + 1] = vi[1] - dt / density_i * bi[1]; + b[3 * i + 2] = vi[2] - dt / density_i * bi[2]; + + // Warmstart + g[3 * i] = vi[0] + m_vDiff[i][0]; + g[3 * i + 1] = vi[1] + m_vDiff[i][1]; + g[3 * i + 2] = vi[2] + m_vDiff[i][2]; + } + } +} + + +#endif + + +void SurfaceTension_Jeske2023::applyForces(const VectorXr &x) { + const int numParticles = (int) m_model->numActiveParticles(); + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + +#pragma omp parallel default(shared) + { +#pragma omp for schedule(static) + for (int i = 0; i < (int) numParticles; i++) { + // Compute the acceleration from the velocity change + Vector3r &ai = m_model->getAcceleration(i); + const Vector3r newVi(x[3 * i], x[3 * i + 1], x[3 * i + 2]); + ai += (1.0 / dt) * (newVi - m_model->getVelocity(i)); + m_vDiff[i] = (newVi - m_model->getVelocity(i)); + } + } +} + + + +void SurfaceTension_Jeske2023::reset() { + std::fill(m_vDiff.begin(), m_vDiff.end(), Vector3r::Zero()); +} + +void SurfaceTension_Jeske2023::performNeighborhoodSearchSort() { + const unsigned int numPart = m_model->numActiveParticles(); + if (numPart == 0) + return; + + Simulation *sim = Simulation::getCurrent(); + auto const& d = sim->getNeighborhoodSearch()->point_set(m_model->getPointSetIndex()); + d.sort_field(&m_vDiff[0]); +} diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Jeske2023.h b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Jeske2023.h new file mode 100644 index 00000000..7ce84b98 --- /dev/null +++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Jeske2023.h @@ -0,0 +1,115 @@ +#ifndef __SurfaceTension_Jeske2023_Surface_Tension_h__ +#define __SurfaceTension_Jeske2023_Surface_Tension_h__ + +#include "SPlisHSPlasH/Common.h" +#include "SPlisHSPlasH/FluidModel.h" +#include "SurfaceTensionBase.h" +#include "SPlisHSPlasH/Utilities/MatrixFreeSolver.h" +#include "Utilities/Logger.h" + +namespace SPH +{ + /** \brief This class implements the implicit surface tension method + * by Jeske et al. 2023 [JWL+23]. + * + * References: + * - [JWL+23] Jeske, Stefan Rhys, Lukas Westhofen, Fabian Löschner, José Antonio Fernández-Fernández, and Jan Bender. “Implicit Surface Tension for SPH Fluid Simulation.” ACM Transactions on Graphics, November 7, 2023. URL: https://doi.org/10.1145/3631936. + + */ + class SurfaceTension_Jeske2023 : public SurfaceTensionBase + { + protected: + Real m_viscosity; + Real m_boundaryViscosity; + + unsigned int m_maxIter; + Real m_maxError; + unsigned int m_iterations; + std::vector m_vDiff; + std::vector m_gradRho; + std::vector m_surfaceEnergy; + std::vector m_color; + std::vector m_colorGrad; + std::vector m_nonlinearAcc; + std::vector m_nonlinearRes; + std::vector m_nonlinearGrad; + + Real m_tangentialDistanceFactor; + bool m_weakPhaseCoupling; + Real m_xsph; + + typedef Eigen::ConjugateGradient Solver; + + Solver m_solver; + + + virtual void initParameters(); + + public: + static int ITERATIONS; + static int MAX_ITERATIONS; + static int MAX_ERROR; + static int VISCOSITY_COEFFICIENT; + static int VISCOSITY_COEFFICIENT_BOUNDARY; + static int XSPH; + + SurfaceTension_Jeske2023(FluidModel *model); + virtual ~SurfaceTension_Jeske2023(void); + + static NonPressureForceBase* creator(FluidModel* model) { return new SurfaceTension_Jeske2023(model); } + + virtual void step(); + virtual void reset(); + + virtual void performNeighborhoodSearchSort(); + + static void matrixVecProd(const Real* vec, Real *result, void *userData); + + FORCE_INLINE const Vector3r& getVDiff(const unsigned int i) const + { + return m_vDiff[i]; + } + + FORCE_INLINE Vector3r& getVDiff(const unsigned int i) + { + return m_vDiff[i]; + } + + FORCE_INLINE void setVDiff(const unsigned int i, const Vector3r& val) + { + m_vDiff[i] = val; + } + + FORCE_INLINE const Real& getDensityGrad(const unsigned int i) const + { + return m_gradRho[i]; + } + + FORCE_INLINE Real& getDensityGrad(const unsigned int i) + { + return m_gradRho[i]; + } + + FORCE_INLINE void setDensityGrad(const unsigned int i, const Real& val) + { + m_gradRho[i] = val; + } + + void computeRHS(VectorXr &b, VectorXr &g); + + void applyForces(const VectorXr &x); + + Real getMaxSolverError(){ return m_maxError; } + void setMaxSolverError(Real error){ m_maxError = error; } + + bool getWeakCoupling(){ return m_weakPhaseCoupling; } + void setWeakCoupling(bool val){ m_weakPhaseCoupling = val; } + + bool getViscosity(){ return m_viscosity; } + void setViscosity(Real val){ m_viscosity = val; } + + void computeDensityGradient(); + }; +} + +#endif diff --git a/SPlisHSPlasH/Utilities/AVX_math.h b/SPlisHSPlasH/Utilities/AVX_math.h index 376fa7c9..674fc4ca 100644 --- a/SPlisHSPlasH/Utilities/AVX_math.h +++ b/SPlisHSPlasH/Utilities/AVX_math.h @@ -153,6 +153,10 @@ static inline Scalarf8 operator >= (Scalarf8 const & a, Scalarf8 const & b) { return _mm256_cmp_ps(b.v, a.v, 2); } +static inline Scalarf8 min(Scalarf8 const & a, Scalarf8 const & b) { + return _mm256_min_ps(a.v, b.v); +} + static inline Scalarf8 max(Scalarf8 const & a, Scalarf8 const & b) { return _mm256_max_ps(a.v, b.v); } diff --git a/data/Scenes/SurfaceTension_BoxHole_ZR2020.json b/data/Scenes/SurfaceTension_BoxHole_ZR2020.json index 32fb3493..b6aa01fc 100644 --- a/data/Scenes/SurfaceTension_BoxHole_ZR2020.json +++ b/data/Scenes/SurfaceTension_BoxHole_ZR2020.json @@ -39,7 +39,7 @@ { "id": "Fluid", "surfaceTension": 0.15, - "surfaceTensionMethod": 4, + "surfaceTensionMethod": 5, "viscosity": 0.01, "viscosityMethod": 1, "surfTZRversion": 1, diff --git a/data/Scenes/SurfaceTension_BreakDamZR2020.json b/data/Scenes/SurfaceTension_BreakDamZR2020.json index 2f34301a..e7dec799 100644 --- a/data/Scenes/SurfaceTension_BreakDamZR2020.json +++ b/data/Scenes/SurfaceTension_BreakDamZR2020.json @@ -24,7 +24,7 @@ { "id": "Fluid", "surfaceTension": 0.1, - "surfaceTensionMethod": 4, + "surfaceTensionMethod": 5, "viscosity": 0.01, "viscosityMethod": 1, "surfTZRversion": 0, diff --git a/data/Scenes/SurfaceTension_CoveredSphere_JWL+23.json b/data/Scenes/SurfaceTension_CoveredSphere_JWL+23.json new file mode 100644 index 00000000..1aed3a8e --- /dev/null +++ b/data/Scenes/SurfaceTension_CoveredSphere_JWL+23.json @@ -0,0 +1,234 @@ +{ + "Configuration": { + "pause": true, + "pauseAt": 20.0, + "stopAt": -1.0, + "sceneScale": 1, + "numberOfStepsPerRenderUpdate": 4, + "renderWalls": 4, + "enablePartioExport": false, + "enableVTKExport": false, + "enableRigidBodyExport": false, + "enableRigidBodyVTKExport": false, + "dataExportFPS": 25.0, + "particleAttributes": "velocity;density", + "enableStateExport": false, + "stateExportFPS": 1.0, + "enableObjectSplitting": true, + "timeStepSize": 0.0010000000474974513, + "particleRadius": 0.02500000037252903, + "sim2D": false, + "enableZSort": true, + "gravitation": [ + 0.0, + -9.8100004196167, + 4.288087325221568e-07 + ], + "maxIterations": 100, + "maxError": 0.009999999776482582, + "boundaryHandlingMethod": 2, + "simulationMethod": 4, + "stiffness": 50000.0, + "exponent": 7.0, + "velocityUpdateMethod": 0, + "enableDivergenceSolver": true, + "maxIterationsV": 100, + "maxErrorV": 0.10000000149011612, + "kernel": 0, + "gradKernel": 0, + "cflMethod": 0, + "cflFactor": 0.8999999761581421, + "cflMinTimeStepSize": 9.999999747378752e-06, + "cflMaxTimeStepSize": 0.009999999776482582, + "cameraPosition": [ + 10.11301040649414, + 5.172463417053223, + 2.1483840942382812 + ], + "cameraLookat": [ + -2.5401246602996252e-05, + -0.20855534076690674, + 0.9780105352401733 + ] + }, + "RigidBodies": [ + { + "tag": "", + "isDynamic": false, + "isWall": false, + "color": [ + 1.0, + 0.0, + 0.0, + 1.0 + ], + "mapInvert": false, + "mapThickness": 0.0, + "mapResolution": [ + 20, + 20, + 20 + ], + "samplingMode": 0, + "translation": [ + 0.0, + 1.9080313444137573, + -0.0 + ], + "scale": [ + 0.45456400513648987, + 0.45456400513648987, + 0.45456400513648987 + ], + "rotationAxis": [ + 0.0, + 0.0, + -1.0 + ], + "rotationAngle": 0.0, + "geometryFile": "../models/sphere.obj", + "boundingBoxMin": [ + -0.45456355810165405, + 1.4534673690795898, + -0.4545634388923645 + ], + "boundingBoxMax": [ + 0.45456385612487793, + 2.362595319747925, + 0.4545637369155884 + ] + }, + { + "tag": "", + "isDynamic": false, + "isWall": true, + "color": [ + 1.0, + 0.0, + 0.0, + 1.0 + ], + "mapInvert": true, + "mapThickness": 0.0, + "mapResolution": [ + 20, + 20, + 20 + ], + "samplingMode": 0, + "translation": [ + 0.0, + 2.0, + -0.0 + ], + "scale": [ + 5.389601707458496, + 5.389601707458496, + 5.389601707458496 + ], + "rotationAxis": [ + 0.0, + 0.0, + -1.0 + ], + "rotationAngle": 0.0, + "geometryFile": "../models/UnitBox.obj", + "boundingBoxMin": [ + -5.389601707458496, + -0.5592994689941406, + -5.389601707458496 + ], + "boundingBoxMax": [ + 5.389601707458496, + 10.219903945922852, + 5.389601707458496 + ] + } + ], + "FluidModels": [], + "FluidBlocks": [], + "Emitters": [ + { + "id": "Fluid_0", + "type": 1, + "velocity": 1.0, + "emitStartTime": 0.0, + "emitEndTime": 100.0, + "tag": "", + "translation": [ + 0.0, + 2.9187259674072266, + -0.0 + ], + "rotationAxis": [ + 0.0, + 0.0, + -1.0 + ], + "rotationAngle": 1.570796012878418, + "width": 6.257507707811778, + "boundingBoxMin": [ + -0.3128753900527954, + 2.918724775314331, + -0.3128753900527954 + ], + "boundingBoxMax": [ + 0.3128753900527954, + 2.918727159500122, + 0.3128753900527954 + ] + } + ], + "AnimationFields": [], + "Materials": [ + { + "id": "Fluid_0", + "density0": 1000.0, + "colorField": "velocity", + "colorMapType": 1, + "renderMinValue": 0.0, + "renderMaxValue": 1.0, + "viscosityMethod": 0, + "viscosity": 0.20000000298023224, + "viscoMaxIter": 100, + "viscoMaxError": 0.009999999776482582, + "viscoMaxIterOmega": 100, + "viscoMaxErrorOmega": 0.009999999776482582, + "viscosityBoundary": 0.0, + "vorticityMethod": 0, + "vorticity": 0.0, + "viscosityOmega": 0.009999999776482582, + "inertiaInverse": 0.0, + "dragMethod": 0, + "drag": 0.0, + "surfaceTensionMethod": 4, + "surfaceTension": 125.00000558793553, + "surfaceTensionBoundary": 1500.0000670552263, + "surfaceTensionViscosity": 0.20000000298023224, + "surfaceTensionViscosityBoundary": 0.0, + "surfaceTensionMaxIter": 100, + "surfaceTensionMaxError": 0.0010000000474974513, + "surfaceTensionMode": 1, + "elasticityMethod": 0, + "youngsModulus": 100000.0, + "poissonsRatio": 0.30000001192092896, + "alpha": 0.0, + "elasticityMaxIter": 100, + "elasticityMaxError": 0.00039999998989515007, + "maxEmitterParticles": 1000000, + "emitterReuseParticles": false, + "emitterBoxMin": [ + 0.0, + 0.0, + 0.0 + ], + "emitterBoxMax": [ + 0.0, + 0.0, + 0.0 + ], + "surfTZRr-ratio": 0.800000011920929, + "surfTZRnormal-mode": 2 + } + ] +} \ No newline at end of file diff --git a/data/Scenes/SurfaceTension_CoveredSphere_ZR2020.json b/data/Scenes/SurfaceTension_CoveredSphere_ZR2020.json index 55536cc6..b74f34b3 100644 --- a/data/Scenes/SurfaceTension_CoveredSphere_ZR2020.json +++ b/data/Scenes/SurfaceTension_CoveredSphere_ZR2020.json @@ -24,7 +24,7 @@ { "id": "Fluid", "surfaceTension": 0.1, - "surfaceTensionMethod": 4, + "surfaceTensionMethod": 5, "renderMaxValue": 13, "viscosity": 0.01, "viscosityMethod": 1, diff --git a/data/Scenes/SurfaceTension_Crown_ZR2020.json b/data/Scenes/SurfaceTension_Crown_ZR2020.json index 8aad029a..b473df2e 100644 --- a/data/Scenes/SurfaceTension_Crown_ZR2020.json +++ b/data/Scenes/SurfaceTension_Crown_ZR2020.json @@ -24,7 +24,7 @@ { "id": "Fluid", "surfaceTension": 0.1, - "surfaceTensionMethod": 4, + "surfaceTensionMethod": 5, "renderMaxValue": 15, "renderMinValue": 0.0, "colorMapType": 1, diff --git a/data/Scenes/SurfaceTension_DoubleDroplet_ZR2020.json b/data/Scenes/SurfaceTension_DoubleDroplet_ZR2020.json index b182fd4f..b096a456 100644 --- a/data/Scenes/SurfaceTension_DoubleDroplet_ZR2020.json +++ b/data/Scenes/SurfaceTension_DoubleDroplet_ZR2020.json @@ -28,7 +28,7 @@ "renderMinValue": 0.0, "colorMapType": 1, "surfaceTension": 0.1 , - "surfaceTensionMethod": 4, + "surfaceTensionMethod": 5, "viscosity": 0.01, "viscosityMethod": 1, "surfTZRversion": 0, diff --git a/data/Scenes/SurfaceTension_FluidChain_JWL+23.json b/data/Scenes/SurfaceTension_FluidChain_JWL+23.json new file mode 100644 index 00000000..10a0f071 --- /dev/null +++ b/data/Scenes/SurfaceTension_FluidChain_JWL+23.json @@ -0,0 +1,142 @@ +{ + "AnimationFields": [], + "Configuration": { + "boundaryHandlingMethod": 2, + "cameraLookat": [ + 0.0, + 0.0, + 0.0 + ], + "cameraPosition": [ + 0.0, + 3.0, + 8.0 + ], + "cflFactor": 0.5, + "cflMaxTimeStepSize": 0.004999999888241291, + "cflMethod": 0, + "cflMinTimeStepSize": 9.999999747378752e-05, + "dataExportFPS": 25.0, + "enableDivergenceSolver": true, + "enablePartioExport": true, + "enableRigidBodyExport": false, + "enableRigidBodyVTKExport": false, + "enableStateExport": false, + "enableVTKExport": false, + "enableZSort": true, + "exponent": 7.0, + "gradKernel": 4, + "gravitation": [ + 0.0, + -9.8100004196167, + 0.0 + ], + "kernel": 4, + "maxError": 0.009999999776482582, + "maxErrorV": 0.009999999776482582, + "maxIterations": 100, + "maxIterationsV": 100, + "numberOfStepsPerRenderUpdate": 4, + "particleAttributes": "velocity", + "particleRadius": 0.01, + "pause": true, + "pauseAt": -1.0, + "renderWalls": 4, + "sim2D": false, + "simulationMethod": 4, + "stateExportFPS": 1.0, + "stiffness": 10000.0, + "stiffnessPF": 50000.0, + "stopAt": -1.0, + "timeStepSize": 0.001, + "velocityUpdateMethod": 0, + "scriptFile_": "two_emitters.py" + }, + "Emitters": [ + { + "emitEndTime": 100000.0, + "emitStartTime": 0.0, + "id": "Fluid", + "rotationAngle": 0.7853981633974484, + "rotationAxis": [ + 0.0, + 0.0, + -1.0 + ], + "translation": [ + -0.5, + 0.0, + 0.0 + ], + "type": 1, + "velocity": 5, + "width": 17 + }, + { + "emitEndTime": 100000.0, + "emitStartTime": 0.0, + "id": "Fluid", + "rotationAngle": -2.356194490192345, + "rotationAxis": [ + 0.0, + 0.0, + 1.0 + ], + "translation": [ + 0.5, + 0.0, + 0.0 + ], + "type": 1, + "velocity": 5, + "width": 17 + } + ], + "FluidModels": [], + "Materials": [ + { + "alpha": 0.0, + "colorField": "velocity", + "colorMapType": 1, + "density0": 1000.0, + "viscosity": 0.0125, + "viscosityMethod": 0, + "drag": 0.009999999776482582, + "dragMethod": 0, + "elasticityMaxError": 9.999999747378752e-05, + "elasticityMaxIter": 100, + "elasticityMethod": 0, + "emitterBoxMax": [ + 10.0, + 10.0, + 10.0 + ], + "emitterBoxMax0": 10.0, + "emitterBoxMax1": 10.0, + "emitterBoxMax2": 10.0, + "emitterBoxMin": [ + -10.0, + -10.0, + -10.0 + ], + "emitterBoxMin0": -10.0, + "emitterBoxMin1": -3.0, + "emitterBoxMin2": -10.0, + "emitterReuseParticles": true, + "id": "Fluid", + "inertiaInverse": 0.5, + "maxEmitterParticles": 10000000, + "poissonsRatio": 0.30000001192092896, + "renderMaxValue": 5.0, + "renderMinValue": 0.0, + "surfaceTensionMethod": 4, + "surfaceTension": 12000.000000000004, + "surfaceTensionBoundary": 0.0, + "surfaceTensionViscosity": 0.0075, + "surfaceTensionViscosityBoundary": 0.0, + "surfaceTensionMaxIter": 100, + "surfaceTensionMaxError": 0.0001 + } + ], + "RigidBodies": [] +} diff --git a/data/Scenes/SurfaceTension_NoGravCube_ZR2020.json b/data/Scenes/SurfaceTension_NoGravCube_ZR2020.json index ec0647cd..5bbdf52e 100644 --- a/data/Scenes/SurfaceTension_NoGravCube_ZR2020.json +++ b/data/Scenes/SurfaceTension_NoGravCube_ZR2020.json @@ -26,7 +26,7 @@ "id": "Fluid", "renderMaxValue": 1.0, "surfaceTension": 0.1, - "surfaceTensionMethod": 4, + "surfaceTensionMethod": 5, "surfTZRversion": 0, "surfTZRCsd": 60000, "surfTZRr-ratio": 0.8, diff --git a/data/Scenes/SurfaceTension_SnakeTestZR2020.json b/data/Scenes/SurfaceTension_SnakeTestZR2020.json index d74f28b4..7f00c57f 100644 --- a/data/Scenes/SurfaceTension_SnakeTestZR2020.json +++ b/data/Scenes/SurfaceTension_SnakeTestZR2020.json @@ -29,7 +29,7 @@ "colorMapType": 1, "colorField": "velocity", "surfaceTension": 0.05 , - "surfaceTensionMethod": 4, + "surfaceTensionMethod": 5, "viscosity": 0.01, "viscosityMethod": 1, " ": 1, diff --git a/data/Scenes/SurfaceTension_WaterBell_JWL+23.json b/data/Scenes/SurfaceTension_WaterBell_JWL+23.json new file mode 100644 index 00000000..9669fb93 --- /dev/null +++ b/data/Scenes/SurfaceTension_WaterBell_JWL+23.json @@ -0,0 +1,171 @@ +{ + "Configuration": { + "pause": true, + "pauseAt": 20.0, + "stopAt": -1.0, + "sceneScale": 1, + "numberOfStepsPerRenderUpdate": 4, + "renderWalls": 4, + "enablePartioExport": false, + "enableVTKExport": false, + "enableRigidBodyExport": false, + "enableRigidBodyVTKExport": false, + "dataExportFPS": 100.0, + "particleAttributes": "velocity;density", + "enableStateExport": false, + "stateExportFPS": 1.0, + "enableObjectSplitting": false, + "timeStepSize": 0.0010000000474974513, + "particleRadius": 0.014999999664723873, + "sim2D": false, + "enableZSort": true, + "gravitation": [ + 0.0, + -9.8100004196167, + 4.288087325221568e-07 + ], + "maxIterations": 100, + "maxError": 0.009999999776482582, + "boundaryHandlingMethod": 2, + "simulationMethod": 4, + "stiffness": 50000.0, + "exponent": 7.0, + "velocityUpdateMethod": 0, + "enableDivergenceSolver": true, + "maxIterationsV": 100, + "maxErrorV": 0.10000000149011612, + "kernel": 0, + "gradKernel": 0, + "cflMethod": 0, + "cflFactor": 0.8999999761581421, + "cflMinTimeStepSize": 9.999999747378752e-06, + "cflMaxTimeStepSize": 0.009999999776482582, + "cameraPosition": [ + 10.11301040649414, + 5.172463417053223, + 2.1483840942382812 + ], + "cameraLookat": [ + -2.5401246602996252e-05, + -0.20855534076690674, + 0.9780105352401733 + ] + }, + "RigidBodies": [], + "FluidModels": [], + "FluidBlocks": [], + "Emitters": [ + { + "id": "Fluid_0", + "type": 1, + "velocity": 8.0, + "emitStartTime": 0.0, + "emitEndTime": 100.0, + "tag": "", + "translation": [ + 0.0, + 2.9187259674072266, + -0.0 + ], + "rotationAxis": [ + 0.0, + 0.0, + -1.0 + ], + "rotationAngle": 1.570796012878418, + "width": 16.666667039195705, + "boundingBoxMin": [ + -0.5, + 2.9187240600585938, + -0.5 + ], + "boundingBoxMax": [ + 0.5, + 2.9187278747558594, + 0.5 + ] + }, + { + "id": "Fluid_0", + "type": 1, + "velocity": 8.0, + "emitStartTime": 0.0, + "emitEndTime": 100.0, + "tag": "", + "translation": [ + 0.0, + 2.017160654067993, + -0.0 + ], + "rotationAxis": [ + 0.0, + 0.0, + 0.9999999403953552 + ], + "rotationAngle": 1.570796251296997, + "width": 16.666667039195705, + "boundingBoxMin": [ + -0.5, + 2.0171587467193604, + -0.5 + ], + "boundingBoxMax": [ + 0.5, + 2.017162561416626, + 0.5 + ] + } + ], + "AnimationFields": [], + "Materials": [ + { + "id": "Fluid_0", + "density0": 1000.0, + "colorField": "velocity", + "colorMapType": 1, + "renderMinValue": 0.0, + "renderMaxValue": 1.0, + "viscosityMethod": 0, + "viscosity": 0.20000000298023224, + "viscoMaxIter": 100, + "viscoMaxError": 0.009999999776482582, + "viscoMaxIterOmega": 100, + "viscoMaxErrorOmega": 0.009999999776482582, + "viscosityBoundary": 0.0, + "vorticityMethod": 0, + "vorticity": 0.0, + "viscosityOmega": 0.009999999776482582, + "inertiaInverse": 0.0, + "dragMethod": 0, + "drag": 0.0, + "surfaceTensionMethod": 4, + "surfaceTension": 10799.999275803582, + "surfaceTensionBoundary": 0.0, + "surfaceTensionViscosity": 0.05, + "surfaceTensionViscosityBoundary": 0.0, + "surfaceTensionMaxIter": 100, + "surfaceTensionMaxError": 0.0010000000474974513, + "surfaceTensionMode": 1, + "elasticityMethod": 0, + "youngsModulus": 100000.0, + "poissonsRatio": 0.30000001192092896, + "alpha": 0.0, + "elasticityMaxIter": 100, + "elasticityMaxError": 0.00039999998989515007, + "maxEmitterParticles": 10000000, + "emitterReuseParticles": false, + "emitterBoxMin": [ + 0.0, + 0.0, + 0.0 + ], + "emitterBoxMax": [ + 0.0, + 0.0, + 0.0 + ], + "surfTZRr-ratio": 0.800000011920929, + "surfTZRnormal-mode": 2 + } + ] +} diff --git a/pySPlisHSPlasH/SurfaceTensionModule.cpp b/pySPlisHSPlasH/SurfaceTensionModule.cpp index cd2b6233..fb9a16b7 100644 --- a/pySPlisHSPlasH/SurfaceTensionModule.cpp +++ b/pySPlisHSPlasH/SurfaceTensionModule.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include @@ -48,4 +49,25 @@ void SurfaceTensionModule(py::module m_sub) { .def("getGradC2", (const Real (SPH::SurfaceTension_He2014::*)(const unsigned int)const)(&SPH::SurfaceTension_He2014::getGradC2)) // .def("getGradC2", (Real& (SPH::SurfaceTension_He2014::*)(const unsigned int))(&SPH::SurfaceTension_He2014::getGradC2)) // TODO: wont work by reference .def("setGradC2", &SPH::SurfaceTension_He2014::setGradC2); + + // --------------------------------------- + // Class Surface Tension Jeske 2023 + // --------------------------------------- + py::class_(m_sub, "SurfaceTension_Jeske2023") + .def(py::init()) + .def_readonly_static("ITERATIONS", &SPH::SurfaceTension_Jeske2023::ITERATIONS) + .def_readonly_static("MAX_ITERATIONS", &SPH::SurfaceTension_Jeske2023::MAX_ITERATIONS) + .def_readonly_static("MAX_ERROR", &SPH::SurfaceTension_Jeske2023::MAX_ERROR) + .def_readonly_static("VISCOSITY_COEFFICIENT", &SPH::SurfaceTension_Jeske2023::VISCOSITY_COEFFICIENT) + .def_readonly_static("VISCOSITY_COEFFICIENT_BOUNDARY", &SPH::SurfaceTension_Jeske2023::VISCOSITY_COEFFICIENT_BOUNDARY) + .def_readonly_static("XSPH", &SPH::SurfaceTension_Jeske2023::XSPH) + + .def("getMaxSolverError", &SPH::SurfaceTension_Jeske2023::getMaxSolverError) + .def("setMaxSolverError", &SPH::SurfaceTension_Jeske2023::setMaxSolverError) + + .def("getWeakCoupling", &SPH::SurfaceTension_Jeske2023::getWeakCoupling) + .def("setWeakCoupling", &SPH::SurfaceTension_Jeske2023::setWeakCoupling) + + .def("getViscosity", &SPH::SurfaceTension_Jeske2023::getViscosity) + .def("setViscosity", &SPH::SurfaceTension_Jeske2023::setViscosity); } \ No newline at end of file diff --git a/version.txt b/version.txt index 064eb29c..a3ebb9f5 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -2.12.5 \ No newline at end of file +2.13.0 \ No newline at end of file