Skip to content

Speedup improvements: faster svd and better stencils #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
larsblatny wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Speedup improvements: faster svd and better stencils #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1aad69a
replaced eigen svd with faster svd
larsblatny Sep 9, 2026
3867168
reducing stencil for quadratic spline
larsblatny Sep 9, 2026
bba5a31
added test for sparse and musl
larsblatny Sep 9, 2026
9407b3f
removed unnecessary init of arrays in g2p
larsblatny Sep 22, 2026
584c781
improved N-R scheme for VMVisc and DPVisc
larsblatny Oct 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions src/simulation/deformation_update.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,22 +17,22 @@ void Simulation::deformationUpdate(){

TM v_grad = TM::Zero();
TV xp = particles.x[p];
unsigned int i_base = std::floor((xp(0)-grid.xc)*one_over_dx) - 1; // the subtraction of one is valid for both quadratic and cubic splines
unsigned int j_base = std::floor((xp(1)-grid.yc)*one_over_dx) - 1;
int i_base = stencilbase((xp(0)-grid.xc)*one_over_dx);
int j_base = stencilbase((xp(1)-grid.yc)*one_over_dx);
#ifdef THREEDIM
unsigned int k_base = std::floor((xp(2)-grid.zc)*one_over_dx) - 1;
int k_base = stencilbase((xp(2)-grid.zc)*one_over_dx);
#endif

for(int i = i_base; i < i_base+4; i++){
for(int i = i_base; i < i_base+STENCILWIDTH; i++){
T xi = grid.x[i];
T wi = N((xp(0)-xi)*one_over_dx);
T wi_grad = dNdu((xp(0) - xi) * one_over_dx) * one_over_dx;
for(int j = j_base; j < j_base+4; j++){
for(int j = j_base; j < j_base+STENCILWIDTH; j++){
T yi = grid.y[j];
T wj = N((xp(1) - yi)*one_over_dx);
T wj_grad = dNdu((xp(1) - yi) * one_over_dx) * one_over_dx;
#ifdef THREEDIM
for(int k = k_base; k < k_base+4; k++){
for(int k = k_base; k < k_base+STENCILWIDTH; k++){
T zi = grid.z[k];
T wk = N((xp(2) - zi)*one_over_dx);
T wk_grad = dNdu((xp(2) - zi) * one_over_dx) * one_over_dx;
Expand Down
30 changes: 13 additions & 17 deletions src/simulation/g2p.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,10 +12,6 @@ void Simulation::G2P(){
debug("G2P");
#endif

std::fill( particles.pic.begin(), particles.pic.end(), TV::Zero() );
std::fill( particles.flip.begin(), particles.flip.end(), TV::Zero() );
std::fill( particles.Bmat.begin(), particles.Bmat.end(), TM::Zero() );

#pragma omp parallel num_threads(n_threads) reduction(+:plastic_count)
{

Expand All @@ -26,22 +22,22 @@ void Simulation::G2P(){
TV flipp = TV::Zero();
TM Bp = TM::Zero();
TM v_grad = TM::Zero();
unsigned int i_base = std::floor((xp(0)-grid.xc)*one_over_dx) - 1; // the subtraction of one is valid for both quadratic and cubic splines
unsigned int j_base = std::floor((xp(1)-grid.yc)*one_over_dx) - 1;
int i_base = stencilbase((xp(0)-grid.xc)*one_over_dx);
int j_base = stencilbase((xp(1)-grid.yc)*one_over_dx);
#ifdef THREEDIM
unsigned int k_base = std::floor((xp(2)-grid.zc)*one_over_dx) - 1;
int k_base = stencilbase((xp(2)-grid.zc)*one_over_dx);
#endif

for(int i = i_base; i < i_base+4; i++){
for(int i = i_base; i < i_base+STENCILWIDTH; i++){
T xi = grid.x[i];
T wi = N((xp(0)-xi)*one_over_dx);
T wi_grad = dNdu((xp(0) - xi) * one_over_dx) * one_over_dx;
for(int j = j_base; j < j_base+4; j++){
for(int j = j_base; j < j_base+STENCILWIDTH; j++){
T yi = grid.y[j];
T wj = N((xp(1) - yi)*one_over_dx);
T wj_grad = dNdu((xp(1) - yi) * one_over_dx) * one_over_dx;
#ifdef THREEDIM
for(int k = k_base; k < k_base+4; k++){
for(int k = k_base; k < k_base+STENCILWIDTH; k++){
T zi = grid.z[k];
T wk = N((xp(2) - zi)*one_over_dx);
T wk_grad = dNdu((xp(2) - zi) * one_over_dx) * one_over_dx;
Expand Down Expand Up @@ -116,7 +112,7 @@ void Simulation::G2PSparseScan() {
unsigned int plastic_count = 0;

const int node_lo_offset = 0;
const int node_hi_offset = 2;
const int node_hi_offset = STENCILWIDTH - 1;

// min/max nodal id
std::vector<int> tminx(n_threads, INT32_MAX);
Expand Down Expand Up @@ -144,9 +140,9 @@ void Simulation::G2PSparseScan() {
TV flipp = TV::Zero();
TM Bp = TM::Zero();
TM v_grad = TM::Zero();
int i_base = std::floor(xp(0)*one_over_dx - 0.5);
int j_base = std::floor(xp(1)*one_over_dx - 0.5);
int k_base = std::floor(xp(2)*one_over_dx - 0.5);
int i_base = stencilbase(xp(0)*one_over_dx);
int j_base = stencilbase(xp(1)*one_over_dx);
int k_base = stencilbase(xp(2)*one_over_dx);

// Get local min/max base
int ilo = i_base + node_lo_offset;
Expand All @@ -164,15 +160,15 @@ void Simulation::G2PSparseScan() {
lmaxz = std::max(lmaxz, khi);

// Local grid nodes loop
for(int ix = i_base+0; ix<i_base+3; ix++){
for(int ix = i_base+0; ix<i_base+STENCILWIDTH; ix++){
T xi = ix * dx;
T wi = N((xp(0)-xi)*one_over_dx);
T wi_grad = dNdu((xp(0) - xi) * one_over_dx) * one_over_dx;
for (int iy=j_base+0; iy<j_base+3; iy++) {
for (int iy=j_base+0; iy<j_base+STENCILWIDTH; iy++) {
T yj = iy * dx;
T wj = N((xp(1) - yj)*one_over_dx);
T wj_grad = dNdu((xp(1) - yj) * one_over_dx) * one_over_dx;
for (int iz=k_base+0; iz<k_base+3; iz++) {
for (int iz=k_base+0; iz<k_base+STENCILWIDTH; iz++) {
T zk = iz * dx;
T wk = N((xp(2) - zk)*one_over_dx);
T wk_grad = dNdu((xp(2) - zk) * one_over_dx) * one_over_dx;
Expand Down
12 changes: 6 additions & 6 deletions src/simulation/musl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,20 +34,20 @@ void Simulation::MUSL(){
#pragma omp for nowait
for(int p = 0; p < Np; p++){
TV xp = particles.x[p];
unsigned int i_base = std::max(0, int(std::floor((xp(0)-grid.xc)*one_over_dx)) - 1); // i_base = std::min(i_base, Nx-4); // the subtraction of one is valid for both quadratic and cubic splines
unsigned int j_base = std::max(0, int(std::floor((xp(1)-grid.yc)*one_over_dx)) - 1); // j_base = std::min(j_base, Ny-4);
int i_base = std::max(0, stencilbase((xp(0)-grid.xc)*one_over_dx));
int j_base = std::max(0, stencilbase((xp(1)-grid.yc)*one_over_dx));
#ifdef THREEDIM
unsigned int k_base = std::max(0, int(std::floor((xp(2)-grid.zc)*one_over_dx)) - 1); // k_base = std::min(k_base, Nz-4);
int k_base = std::max(0, stencilbase((xp(2)-grid.zc)*one_over_dx));
#endif

for(int i = i_base; i < i_base+4; i++){
for(int i = i_base; i < i_base+STENCILWIDTH; i++){
T xi = grid.x[i];
T wi = N((xp(0)-xi)*one_over_dx);
for(int j = j_base; j < j_base+4; j++){
for(int j = j_base; j < j_base+STENCILWIDTH; j++){
T yi = grid.y[j];
T wj = N((xp(1) - yi)*one_over_dx);
#ifdef THREEDIM
for(int k = k_base; k < k_base+4; k++){
for(int k = k_base; k < k_base+STENCILWIDTH; k++){
T zi = grid.z[k];
T wk = N((xp(2) - zi)*one_over_dx);
T weight = wi * wj * wk;
Expand Down
24 changes: 12 additions & 12 deletions src/simulation/p2g.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,10 +19,10 @@ void Simulation::P2G(){
#pragma omp for nowait
for(int p = 0; p < Np; p++){
TV xp = particles.x[p];
unsigned int i_base = std::max(0, int(std::floor((xp(0)-grid.xc)*one_over_dx)) - 1); // i_base = std::min(i_base, Nx-4); // the subtraction of one is valid for both quadratic and cubic splines
unsigned int j_base = std::max(0, int(std::floor((xp(1)-grid.yc)*one_over_dx)) - 1); // j_base = std::min(j_base, Ny-4);
int i_base = std::max(0, stencilbase((xp(0)-grid.xc)*one_over_dx));
int j_base = std::max(0, stencilbase((xp(1)-grid.yc)*one_over_dx));
#ifdef THREEDIM
unsigned int k_base = std::max(0, int(std::floor((xp(2)-grid.zc)*one_over_dx)) - 1); // k_base = std::min(k_base, Nz-4);
int k_base = std::max(0, stencilbase((xp(2)-grid.zc)*one_over_dx));
#endif

// ALT 1: Compute stress
Expand All @@ -42,16 +42,16 @@ void Simulation::P2G(){
// ALT 2: Use the saved stress
TM tau = particles.tau[p];

for(int i = i_base; i < i_base+4; i++){
for(int i = i_base; i < i_base+STENCILWIDTH; i++){
T xi = grid.x[i];
T wi = N((xp(0)-xi)*one_over_dx);
T wi_grad = dNdu((xp(0) - xi) * one_over_dx) * one_over_dx;
for(int j = j_base; j < j_base+4; j++){
for(int j = j_base; j < j_base+STENCILWIDTH; j++){
T yi = grid.y[j];
T wj = N((xp(1) - yi)*one_over_dx);
T wj_grad = dNdu((xp(1) - yi) * one_over_dx) * one_over_dx;
#ifdef THREEDIM
for(int k = k_base; k < k_base+4; k++){
for(int k = k_base; k < k_base+STENCILWIDTH; k++){
T zi = grid.z[k];
T wk = N((xp(2) - zi)*one_over_dx);
T wk_grad = dNdu((xp(2) - zi) * one_over_dx) * one_over_dx;
Expand Down Expand Up @@ -150,23 +150,23 @@ void Simulation::P2GSparseScan() {
#pragma omp for nowait
for(int p = 0; p < Np; p++){
TV xp = particles.x[p];
int i_base = std::floor(xp(0)*one_over_dx - 0.5);
int j_base = std::floor(xp(1)*one_over_dx - 0.5);
int k_base = std::floor(xp(2)*one_over_dx - 0.5);
int i_base = stencilbase(xp(0)*one_over_dx);
int j_base = stencilbase(xp(1)*one_over_dx);
int k_base = stencilbase(xp(2)*one_over_dx);

// Kirchhoff stress
TM tau = particles.tau[p];

// Grid loop for each particle
for(int ix = i_base+0; ix<i_base+3; ix++){
for(int ix = i_base+0; ix<i_base+STENCILWIDTH; ix++){
T xi = ix * dx;
T wi = N((xp(0)-xi)*one_over_dx);
T wi_grad = dNdu((xp(0) - xi) * one_over_dx) * one_over_dx;
for (int iy=j_base+0; iy<j_base+3; iy++) {
for (int iy=j_base+0; iy<j_base+STENCILWIDTH; iy++) {
T yj = iy * dx;
T wj = N((xp(1) - yj)*one_over_dx);
T wj_grad = dNdu((xp(1) - yj) * one_over_dx) * one_over_dx;
for (int iz=k_base+0; iz<k_base+3; iz++) {
for (int iz=k_base+0; iz<k_base+STENCILWIDTH; iz++) {
T zk = iz * dx;
T wk = N((xp(2) - zk)*one_over_dx);
T wk_grad = dNdu((xp(2) - zk) * one_over_dx) * one_over_dx;
Expand Down
Loading
Loading