Can’t get a solution on a large dimension

Now I am trying to work with your library and I have a problem. I solve the diffusion equation by the implicit Galerkin method. I want to get a solution on a triangular grid. When I take 50 cells and less, I get a completely correct solution. But when I increase the number of cells, I immediately get only the zero solution.

The calculations are organized in such a way that the matrix of the system is static, and the right-hand side changes at each iteration. I store the matrix coefficients in CSR format and consider the matrix scalar, as you do in the documentation. That is, I use 1x1 blocks, 1 element in each block.

I tried to use different configurations: V.json, FGMRES.json and others. Some do not give a solution even on small grids. But FGMRES.json works fine when the number of cells is less than 50.

I think my matrix is correct. I believe that I am working incorrectly with your library. And I need help to figure it out. I attach the code in which I work with your library. What am I doing wrong?

I’m using CUDA 10.0 and CUDA Runtime 10.0 I tried CUDA 9, but it also does not work.

My code:

`float* b_data = new float[A_block_size * cellsCount];
float* lastLayerSolution = new float[A_block_size * cellsCount];
float* currentLayerSolution = new float[A_block_size * cellsCount];

memset(b_data, 0.0, sizeof(float) * A_block_size * cellsCount);
memset(lastLayerSolution, 0.0, sizeof(float) * A_block_size * cellsCount);
memset(currentLayerSolution, 0.0, sizeof(float) * A_block_size * cellsCount);

for (int iCell = 0; iCell < cellsCount; iCell++) {
	lastLayerSolution[A_block_size * iCell + 6] = sin(cellCx[iCell] * M_PI) * sin(cellCy[iCell] * M_PI);
}

// begin AMGX init
AMGX_initialize();
AMGX_initialize_plugins();

AMGX_config_handle config;
AMGX_config_create_from_file(&config, "configs/config.json");

AMGX_resources_handle rsrc;
AMGX_resources_create_simple(&rsrc, config);

AMGX_solver_handle solver;
AMGX_matrix_handle A_amgx;
AMGX_vector_handle b_amgx;
AMGX_vector_handle solution_amgx;

AMGX_solver_create(&solver, rsrc, AMGX_mode_dFFI, config);
AMGX_matrix_create(&A_amgx, rsrc, AMGX_mode_dFFI);
AMGX_vector_create(&b_amgx, rsrc, AMGX_mode_dFFI);
AMGX_vector_create(&solution_amgx, rsrc, AMGX_mode_dFFI);
// end AMGX init

CSRMatrix matrix(A_block_size * cellsCount);

initRightPart(lastLayerSolution, b_data);

initMatrix(matrix);

matrix.printToFile("resMatrix.txt");

int n_amgx = A_block_size * cellsCount;
int nnz_amgx = matrix.na;

AMGX_pin_memory(currentLayerSolution, n_amgx * sizeof(float));
AMGX_pin_memory(b_data, n_amgx * sizeof(float));

AMGX_matrix_upload_all(A_amgx, n_amgx, nnz_amgx, 1, 1, matrix.ia, matrix.ja, matrix.a, 0);
AMGX_vector_upload(b_amgx, n_amgx, 1, b_data);
AMGX_vector_set_zero(solution_amgx, n_amgx, 1);
AMGX_solver_setup(solver, A_amgx);

writeVTK_AMGX((char*)TASK_NAME, step, lastLayerSolution);

while (t < TMAX) {
	t += TAU;
	step++;

	// ************ begin AMGX solver ************

	AMGX_solver_solve_with_0_initial_guess(solver, b_amgx, solution_amgx);

	AMGX_vector_download(solution_amgx, currentLayerSolution);

	copyToOld(lastLayerSolution, currentLayerSolution, n_amgx);

	initRightPart(lastLayerSolution, b_data);

	AMGX_vector_upload(b_amgx, n_amgx, 1, b_data);

	if (step % SAVE_STEP == 0)
	{
		writeVTK_AMGX((char*)TASK_NAME, step, currentLayerSolution);

		double L2 = 0.0;

		for (int i = 0; i < cellsCount; i++) {
			L2 += cellS[i] * (currentLayerSolution[i * A_block_size + 6] - exact(t, cellCx[i], cellCy[i])) * (currentLayerSolution[i * A_block_size + 6] - exact(t, cellCx[i], cellCy[i]));
		}

		OutPut("L2 = %15.10f\n", sqrt(L2));
	}
	if (step % 100 == 0)
	{
		OutPut("%8d | time: %15.8e | U: %15.8e | WX: %15.8e | WY: %15.8e |\n", step, t, resU, resWx, resWy);
	}

	if (step == 700)
	{
		printf("runtime = %f", (clock() / 1000.0));
	}

	memset(currentLayerSolution, 0.0, sizeof(float) * A_block_size * cellsCount);
	AMGX_vector_set_zero(solution_amgx, n_amgx, 1);
}

AMGX_solver_destroy(solver);
AMGX_vector_destroy(b_amgx);
AMGX_vector_destroy(solution_amgx);
AMGX_matrix_destroy(A_amgx);
AMGX_resources_destroy(rsrc);

AMGX_finalize_plugins();
AMGX_finalize();`

Start by reviewing configs/config.json and the attached code around AMGX_matrix_upload_all, AMGX_vector_upload, AMGX_solver_setup, and AMGX_solver_solve_with_0_initial_guess. Reproduce the behavior with the reported grid sizes and CUDA versions, then identify why larger systems produce a zero solution and verify a correct nonzero solution on those grids.