Swap plane Poisson solver with slab Poisson solver in GK IWL simulations & 2x LBO energy conservation

core/unit/ctest_cudss.cu

@@ -11,6 +11,7 @@ extern "C" {
 extern "C" {
 void test_cudss_simple();
 void test_cudss_ops();
+void test_cudss_ops_update_amat();
 void test_cudss_ops_multiple_rhs();
 }
 
@@ -259,6 +260,109 @@ void test_cudss_ops()
   gkyl_culinsolver_prob_release(prob);
 }
 
+void test_cudss_ops_update_amat()
+{
+  int nfail = 0;  
+  double s, u, p, e, r, l;
+  int    nrhs, m, n;
+
+  /* Initialize matrix A. */
+  /*  A : matrix([s,0,u,u,0],[l,u,0,0,0],[0,l,p,0,0],[0,0,0,e,u],[l,l,0,0,r]); */
+  m = n = 5;
+  nrhs = 1;
+
+  s = 19.0; u = 21.0; p = 16.0; e = 5.0; r = 18.0; l = 12.0;
+  /*  A : matrix([s,0,u,u,0],[l,u,0,0,0],[0,l,p,0,0],[0,0,0,e,u],[l,l,0,0,r]); */
+  struct gkyl_mat_triples **tri_arr = (struct gkyl_mat_triples **) gkyl_malloc(sizeof(struct gkyl_mat_triples *));
+  tri_arr[0] = gkyl_mat_triples_new(m, n);
+  struct gkyl_mat_triples *tri = tri_arr[0];
+  gkyl_mat_triples_set_rowmaj_order(tri);
+  // row 0
+  gkyl_mat_triples_insert(tri, 0, 0, s);
+  gkyl_mat_triples_insert(tri, 0, 2, u);
+  gkyl_mat_triples_insert(tri, 0, 3, u);
+  // row 1
+  gkyl_mat_triples_insert(tri, 1, 0, l);
+  gkyl_mat_triples_insert(tri, 1, 1, u);
+  // row 2
+  gkyl_mat_triples_insert(tri, 2, 1, l);
+  gkyl_mat_triples_insert(tri, 2, 2, p);
+  // row 3
+  gkyl_mat_triples_insert(tri, 3, 3, e);
+  gkyl_mat_triples_insert(tri, 3, 4, u);
+  // row 4
+  gkyl_mat_triples_insert(tri, 4, 0, l);
+  gkyl_mat_triples_insert(tri, 4, 1, l);
+  gkyl_mat_triples_insert(tri, 4, 4, r);
+
+  // Create the cuSolver linear problem setup.
+  gkyl_culinsolver_prob *prob = gkyl_culinsolver_prob_new(1, m, n, nrhs);
+
+  // Allocate the A matrix from triples.
+  gkyl_culinsolver_amat_from_triples(prob, tri_arr);
+
+  // Create right-hand side matrix B = transpose([1,1,1,1,1]).
+  gkyl_mat_triples *triRHS = gkyl_mat_triples_new(m, nrhs);
+  gkyl_mat_triples_insert(triRHS, 0, 0, 1.0);
+  gkyl_mat_triples_insert(triRHS, 1, 0, 1.0);
+  gkyl_mat_triples_insert(triRHS, 2, 0, 1.0);
+  gkyl_mat_triples_insert(triRHS, 3, 0, 1.0);
+  gkyl_mat_triples_insert(triRHS, 4, 0, 1.0);
+  gkyl_culinsolver_brhs_from_triples(prob, triRHS);
+
+  gkyl_culinsolver_solve(prob);
+  gkyl_culinsolver_finish_host(prob);
+
+  // Solution is: [-1/32, 11/168, 3/224, 1/16, 11/336].
+  GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0  , gkyl_culinsolver_get_sol_lin(prob,0), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 11.0/168.0, gkyl_culinsolver_get_sol_lin(prob,1), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0 , gkyl_culinsolver_get_sol_lin(prob,2), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0  , gkyl_culinsolver_get_sol_lin(prob,3), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 11.0/336.0, gkyl_culinsolver_get_sol_lin(prob,4), 1e-14), &nfail );
+
+  // Now update the LHS matrix. Multiply it by a constant so the solution should be the same as before but divided by that constant.
+  double prob_fac = 1.3;
+
+  gkyl_mat_triples_set_rowmaj_order(tri);
+  // row 0
+  gkyl_mat_triples_insert(tri, 0, 0, s*prob_fac);
+  gkyl_mat_triples_insert(tri, 0, 2, u*prob_fac);
+  gkyl_mat_triples_insert(tri, 0, 3, u*prob_fac);
+  // row 1
+  gkyl_mat_triples_insert(tri, 1, 0, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 1, 1, u*prob_fac);
+  // row 2
+  gkyl_mat_triples_insert(tri, 2, 1, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 2, 2, p*prob_fac);
+  // row 3
+  gkyl_mat_triples_insert(tri, 3, 3, e*prob_fac);
+  gkyl_mat_triples_insert(tri, 3, 4, u*prob_fac);
+  // row 4
+  gkyl_mat_triples_insert(tri, 4, 0, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 4, 1, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 4, 4, r*prob_fac);
+
+  gkyl_culinsolver_amat_update_from_triples(prob, tri_arr);
+
+  // Reset RHS for good measure.
+  gkyl_culinsolver_brhs_from_triples(prob, triRHS);
+
+  gkyl_culinsolver_solve(prob);
+  gkyl_culinsolver_finish_host(prob);
+
+  // Solution is: (1/prob_fac)*[-1/32, 11/168, 3/224, 1/16, 11/336].
+  GKYL_CU_CHECK( gkyl_compare_double((1.0/prob_fac)*(-1.0/32.0  ), gkyl_culinsolver_get_sol_lin(prob,0), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double((1.0/prob_fac)*( 11.0/168.0), gkyl_culinsolver_get_sol_lin(prob,1), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double((1.0/prob_fac)*( 3.0/224.0 ), gkyl_culinsolver_get_sol_lin(prob,2), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double((1.0/prob_fac)*( 1.0/16.0  ), gkyl_culinsolver_get_sol_lin(prob,3), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double((1.0/prob_fac)*( 11.0/336.0), gkyl_culinsolver_get_sol_lin(prob,4), 1e-14), &nfail );
+
+  gkyl_mat_triples_release(tri);
+  gkyl_free(tri_arr);
+  gkyl_mat_triples_release(triRHS);
+  gkyl_culinsolver_prob_release(prob);
+}
+
 void test_cudss_ops_multiple_rhs()
 {
   double s, u, p, e, r, l;
@@ -320,26 +424,152 @@ void test_cudss_ops_multiple_rhs()
 
   // Solution is: [-1/32, 11/168, 3/224, 1/16, 11/336].
   // 1st problem
-  GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0,   gkyl_culinsolver_get_sol_lin(prob,0), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0  , gkyl_culinsolver_get_sol_lin(prob,0), 1e-14), &nfail );
   GKYL_CU_CHECK( gkyl_compare_double( 11.0/168.0, gkyl_culinsolver_get_sol_lin(prob,1), 1e-14), &nfail );
-  GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0,  gkyl_culinsolver_get_sol_lin(prob,2), 1e-14), &nfail );
-  GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0,   gkyl_culinsolver_get_sol_lin(prob,3), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0 , gkyl_culinsolver_get_sol_lin(prob,2), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0  , gkyl_culinsolver_get_sol_lin(prob,3), 1e-14), &nfail );
   GKYL_CU_CHECK( gkyl_compare_double( 11.0/336.0, gkyl_culinsolver_get_sol_lin(prob,4), 1e-14), &nfail );
   // 2nd problem
-  GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0,   gkyl_culinsolver_get_sol_lin(prob,5), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0  , gkyl_culinsolver_get_sol_lin(prob,5), 1e-14), &nfail );
   GKYL_CU_CHECK( gkyl_compare_double( 11.0/168.0, gkyl_culinsolver_get_sol_lin(prob,6), 1e-14), &nfail );
-  GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0,  gkyl_culinsolver_get_sol_lin(prob,7), 1e-14), &nfail );
-  GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0,   gkyl_culinsolver_get_sol_lin(prob,8), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0 , gkyl_culinsolver_get_sol_lin(prob,7), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0  , gkyl_culinsolver_get_sol_lin(prob,8), 1e-14), &nfail );
   GKYL_CU_CHECK( gkyl_compare_double( 11.0/336.0, gkyl_culinsolver_get_sol_lin(prob,9), 1e-14), &nfail );
   // 3rd problem
-  GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0,   gkyl_culinsolver_get_sol_lin(prob,10), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0  , gkyl_culinsolver_get_sol_lin(prob,10), 1e-14), &nfail );
   GKYL_CU_CHECK( gkyl_compare_double( 11.0/168.0, gkyl_culinsolver_get_sol_lin(prob,11), 1e-14), &nfail );
-  GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0,  gkyl_culinsolver_get_sol_lin(prob,12), 1e-14), &nfail );
-  GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0,   gkyl_culinsolver_get_sol_lin(prob,13), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0 , gkyl_culinsolver_get_sol_lin(prob,12), 1e-14), &nfail );
+  GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0  , gkyl_culinsolver_get_sol_lin(prob,13), 1e-14), &nfail );
   GKYL_CU_CHECK( gkyl_compare_double( 11.0/336.0, gkyl_culinsolver_get_sol_lin(prob,14), 1e-14), &nfail );
 
   gkyl_culinsolver_prob_release(prob);
 }
 
+void test_cudss_ops_update_amat_multiple_rhs()
+{
+  double s, u, p, e, r, l;
+  int    nrhs, m, n;
+
+  /* Initialize matrix A. */
+  /*  A : matrix([s,0,u,u,0],[l,u,0,0,0],[0,l,p,0,0],[0,0,0,e,u],[l,l,0,0,r]); */
+  m = n = 5;
+  nrhs = 3;
+
+  s = 19.0; u = 21.0; p = 16.0; e = 5.0; r = 18.0; l = 12.0;
+  /*  A : matrix([s,0,u,u,0],[l,u,0,0,0],[0,l,p,0,0],[0,0,0,e,u],[l,l,0,0,r]); */
+  struct gkyl_mat_triples **tri_arr = (struct gkyl_mat_triples **) gkyl_malloc(sizeof(struct gkyl_mat_triples *));
+  tri_arr[0] = gkyl_mat_triples_new(m, n);
+  struct gkyl_mat_triples *tri = tri_arr[0];
+  gkyl_mat_triples_set_rowmaj_order(tri);
+  // row 0
+  gkyl_mat_triples_insert(tri, 0, 0, s);
+  gkyl_mat_triples_insert(tri, 0, 2, u);
+  gkyl_mat_triples_insert(tri, 0, 3, u);
+  // row 1
+  gkyl_mat_triples_insert(tri, 1, 0, l);
+  gkyl_mat_triples_insert(tri, 1, 1, u);
+  // row 2
+  gkyl_mat_triples_insert(tri, 2, 1, l);
+  gkyl_mat_triples_insert(tri, 2, 2, p);
+  // row 3
+  gkyl_mat_triples_insert(tri, 3, 3, e);
+  gkyl_mat_triples_insert(tri, 3, 4, u);
+  // row 4
+  gkyl_mat_triples_insert(tri, 4, 0, l);
+  gkyl_mat_triples_insert(tri, 4, 1, l);
+  gkyl_mat_triples_insert(tri, 4, 4, r);
+
+  // Create the cuSolver linear problem setup.
+  gkyl_culinsolver_prob *prob = gkyl_culinsolver_prob_new(1, m, n, nrhs);
+
+  // Allocate the A matrix from triples.
+  gkyl_culinsolver_amat_from_triples(prob, tri_arr);
+
+  // Create right-hand side matrix B = transpose([1,1,1,1,1]).
+  gkyl_mat_triples *triRHS = gkyl_mat_triples_new(m, nrhs);
+  for (int k=0; k<nrhs; k++) {
+    gkyl_mat_triples_insert(triRHS, 0, k, 1.0);
+    gkyl_mat_triples_insert(triRHS, 1, k, 1.0);
+    gkyl_mat_triples_insert(triRHS, 2, k, 1.0);
+    gkyl_mat_triples_insert(triRHS, 3, k, 1.0);
+    gkyl_mat_triples_insert(triRHS, 4, k, 1.0);
+  }
+  gkyl_culinsolver_brhs_from_triples(prob, triRHS);
+
+  gkyl_culinsolver_solve(prob);
+  gkyl_culinsolver_finish_host(prob);
+
+  int nfail = 0;
+
+  // Solution is: [-1/32, 11/168, 3/224, 1/16, 11/336].
+  for (int k=0; k<nrhs; k++) {
+    int off = m*k;
+    GKYL_CU_CHECK( gkyl_compare_double(-1.0/32.0  , gkyl_culinsolver_get_sol_lin(prob,off+0), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double( 11.0/168.0, gkyl_culinsolver_get_sol_lin(prob,off+1), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double( 3.0/224.0 , gkyl_culinsolver_get_sol_lin(prob,off+2), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double( 1.0/16.0  , gkyl_culinsolver_get_sol_lin(prob,off+3), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double( 11.0/336.0, gkyl_culinsolver_get_sol_lin(prob,off+4), 1e-14), &nfail );
+  }
+
+  // Now update the LHS matrix. Multiply it by a constant, and multiply the RHS
+  // by 2X that constant, so the solution should be the same but divided by 2.
+  double prob_fac = 1.3;
+  double *prob_fac_per_rhs = (double *) gkyl_malloc(nrhs*sizeof(double));
+  for (int k=0; k<nrhs; k++)
+    prob_fac_per_rhs[k] = (k+1)*prob_fac;
+
+  gkyl_mat_triples_set_rowmaj_order(tri);
+  // row 0
+  gkyl_mat_triples_insert(tri, 0, 0, s*prob_fac);
+  gkyl_mat_triples_insert(tri, 0, 2, u*prob_fac);
+  gkyl_mat_triples_insert(tri, 0, 3, u*prob_fac);
+  // row 1
+  gkyl_mat_triples_insert(tri, 1, 0, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 1, 1, u*prob_fac);
+  // row 2
+  gkyl_mat_triples_insert(tri, 2, 1, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 2, 2, p*prob_fac);
+  // row 3
+  gkyl_mat_triples_insert(tri, 3, 3, e*prob_fac);
+  gkyl_mat_triples_insert(tri, 3, 4, u*prob_fac);
+  // row 4
+  gkyl_mat_triples_insert(tri, 4, 0, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 4, 1, l*prob_fac);
+  gkyl_mat_triples_insert(tri, 4, 4, r*prob_fac);
+
+  // Create the cuSolver linear problem setup.
+  gkyl_culinsolver_amat_update_from_triples(prob, tri_arr);
+
+  // Reset RHS.
+  for (int k=0; k<nrhs; k++) {
+    gkyl_mat_triples_insert(triRHS, 0, k, 1.0*prob_fac_per_rhs[k]);
+    gkyl_mat_triples_insert(triRHS, 1, k, 1.0*prob_fac_per_rhs[k]);
+    gkyl_mat_triples_insert(triRHS, 2, k, 1.0*prob_fac_per_rhs[k]);
+    gkyl_mat_triples_insert(triRHS, 3, k, 1.0*prob_fac_per_rhs[k]);
+    gkyl_mat_triples_insert(triRHS, 4, k, 1.0*prob_fac_per_rhs[k]);
+  }
+  gkyl_culinsolver_brhs_from_triples(prob, triRHS);
+
+  gkyl_culinsolver_solve(prob);
+  gkyl_culinsolver_finish_host(prob);
+
+  // Solution is: (k+1)*[-1/32, 11/168, 3/224, 1/16, 11/336].
+  for (int k=0; k<nrhs; k++) {
+    int off = m*k;
+    double fac = prob_fac_per_rhs[k]/prob_fac;
+    GKYL_CU_CHECK( gkyl_compare_double(fac*(-1.0/32.0  ), gkyl_culinsolver_get_sol_lin(prob,off+0), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double(fac*( 11.0/168.0), gkyl_culinsolver_get_sol_lin(prob,off+1), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double(fac*( 3.0/224.0 ), gkyl_culinsolver_get_sol_lin(prob,off+2), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double(fac*( 1.0/16.0  ), gkyl_culinsolver_get_sol_lin(prob,off+3), 1e-14), &nfail );
+    GKYL_CU_CHECK( gkyl_compare_double(fac*( 11.0/336.0), gkyl_culinsolver_get_sol_lin(prob,off+4), 1e-14), &nfail );
+  }
+
+  gkyl_free(prob_fac_per_rhs);
+  gkyl_mat_triples_release(tri);
+  gkyl_free(tri_arr);
+  gkyl_mat_triples_release(triRHS);
+  gkyl_culinsolver_prob_release(prob);
+}
+
 // End ifdef GKYL_HAVE_CUDSS statement.
 #endif