OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-08-04 18:43:21 -05:00
Code Issues Packages Projects Releases Wiki Activity

Removed CPU implementation of WellContributions

Browse Source
This commit is contained in:
T.D. (Tongdong) Qiu
2020-03-18 15:08:48 +01:00
parent ce73337769
commit 009b5076a4
4 changed files with 12 additions and 239 deletions
Download Patch File Download Diff File Expand all files Collapse all files
opm/simulators/linalg/bda/BdaBridge.cpp
-1
View File
@@ -41,7 +41,6 @@ BdaBridge::BdaBridge(bool use_gpu_, int linear_solver_verbosity, int maxit, doub
{
if (use_gpu) {
backend.reset(new cusparseSolverBackend(linear_solver_verbosity, maxit, tolerance));
WellContributions::setMode(use_gpu);
}
}
#else
opm/simulators/linalg/bda/WellContributions.cu
+10 -208
View File
@@ -129,32 +129,8 @@ namespace Opm
}
#endif
void WellContributions::alloc_all(){
#if HAVE_CUDA
if(gpu_mode){
alloc_gpu();
}else{
alloc_cpu();
}
#else
alloc_cpu();
#endif
allocated = true;
}
void WellContributions::alloc_cpu(){
Cnnzs = new double[num_blocks * dim * dim_wells];
Dnnzs = new double[num_wells * dim_wells * dim_wells];
Bnnzs = new double[num_blocks * dim * dim_wells];
Ccols = new int[num_blocks];
Bcols = new int[num_blocks];
val_pointers = new unsigned int[num_wells + 1];
z1 = new double[dim_wells]; // B * x
z2 = new double[dim_wells]; // D^-1 * B * x
}
#if HAVE_CUDA
void WellContributions::alloc_gpu(){
void WellContributions::alloc(){
cudaMalloc((void**)&d_Cnnzs, sizeof(double) * num_blocks * dim * dim_wells);
cudaMalloc((void**)&d_Dnnzs, sizeof(double) * num_wells * dim_wells * dim_wells);
cudaMalloc((void**)&d_Bnnzs, sizeof(double) * num_blocks * dim * dim_wells);
@@ -163,36 +139,11 @@ namespace Opm
val_pointers = new unsigned int[num_wells + 1];
cudaMalloc((void**)&d_val_pointers, sizeof(int) * (num_wells + 1));
cudaCheckLastError("apply_gpu malloc failed");
allocated = true;
}
#endif
WellContributions::~WellContributions()
{
#if HAVE_CUDA
if(gpu_mode){
free_gpu();
}else{
free_cpu();
}
#else
free_cpu();
#endif
}
void WellContributions::free_cpu(){
delete[] Cnnzs;
delete[] Dnnzs;
delete[] Bnnzs;
delete[] Ccols;
delete[] Bcols;
delete[] val_pointers;
delete[] z1;
delete[] z2;
//delete[] Mnnzs;
}
#if HAVE_CUDA
void WellContributions::free_gpu(){
cudaFree(d_Cnnzs);
cudaFree(d_Dnnzs);
cudaFree(d_Bnnzs);
@@ -200,166 +151,20 @@ namespace Opm
cudaFree(d_Bcols);
delete[] val_pointers;
cudaFree(d_val_pointers);
// cudaFree(d_z1);
// cudaFree(d_z2);
}
#endif
void WellContributions::apply(double *x, double *y){
#if HAVE_CUDA
if (gpu_mode){
apply_gpu(x, y);
}else{
apply_cpu(x, y);
}
#else
apply_cpu(x, y);
#endif
}
// Apply the WellContributions, similar to StandardWell::apply()
// y -= (C^T *(D^-1*( B*x)))
void WellContributions::apply_cpu(double *x, double *y)
{
#if 0
// Mnnzs contains a sparse matrix with a symmetric pattern
// Mrows would contain 'i*val_size' for every entry i, since every row has the same number of blocks
// Mcols are the same as Ccols, normally, there is an entry for every block, but since all rows have the same sparsity pattern, we only have to store 1 row
bool dbg = false;
for(int i = 0; i < dim*dim*val_size*val_size; ++i){
if(dbg)printf("Mnnzs[%d]: %.5e\n", i, Mnnzs[i]);
}
if(dbg)printf("row_size: %u, val_size: %u\n", row_size, val_size);
for(int r = 0; r < val_size; ++r){
for(int c = 0; c < val_size; ++c){
int colIdx = Ccols[c];
if(dbg)printf("colIdx: %d\n", colIdx);
for(int i = 0; i < dim; ++i){
double sum = 0.0;
for(int j = 0; j < dim; ++j){
sum += Mnnzs[r * dim * dim * val_size + c * dim * dim + i * dim + j] * x[colIdx * dim + j];
}
if(dbg)printf("sum: %f\n", sum);
y[colIdx * dim + i] -= sum;
}
}
}
if(dbg)exit(0);
#else
for(int wellID = 0; wellID < num_wells; ++wellID){
unsigned int val_size = val_pointers[wellID+1] - val_pointers[wellID];
// B * x
for (unsigned int i = 0; i < dim_wells; ++i) {
z1[i] = 0.0;
}
for (unsigned int i = 0; i < val_size; ++i) {
unsigned int blockID = i + val_pointers[wellID];
int colIdx = Bcols[blockID];
for (unsigned int j = 0; j < dim_wells; ++j) {
double temp = 0.0;
for (unsigned int k = 0; k < dim; ++k) {
temp += Bnnzs[blockID * dim * dim_wells + j * dim + k] * x[colIdx * dim + k];
}
z1[j] += temp;
}
}
// D^-1 * B * x
for (unsigned int i = 0; i < dim_wells; ++i) {
z2[i] = 0.0;
}
for (unsigned int j = 0; j < dim_wells; ++j) {
double temp = 0.0;
for (unsigned int k = 0; k < dim_wells; ++k) {
temp += Dnnzs[wellID * dim_wells * dim_wells + j * dim_wells + k] * z1[k];
}
z2[j] += temp;
}
// C^T * D^-1 * B * x
for (unsigned int i = 0; i < val_size; ++i) {
unsigned int blockID = i + val_pointers[wellID];
int colIdx = Ccols[blockID];
for (unsigned int j = 0; j < dim; ++j) {
double temp = 0.0;
for (unsigned int k = 0; k < dim_wells; ++k) {
temp += Cnnzs[blockID * dim * dim_wells + j + k * dim] * z2[k];
}
y[colIdx * dim + j] -= temp;
}
}
}
#endif
}
// Apply the WellContributions, similar to StandardWell::apply()
// y -= (C^T *(D^-1*( B*x)))
#if HAVE_CUDA
void WellContributions::apply_gpu(double *d_x, double *d_y)
void WellContributions::apply(double *d_x, double *d_y)
{
int smem_size = 2 * sizeof(double) * dim_wells;
apply_well_contributions<<<num_wells, 32, smem_size, stream>>>(d_Cnnzs, d_Dnnzs, d_Bnnzs, d_Ccols, d_Bcols, d_x, d_y, dim, dim_wells, d_val_pointers);
}
#endif
void WellContributions::addMatrix(int idx, int *colIndices, double *values, unsigned int val_size)
{
#if HAVE_CUDA
if(gpu_mode){
addMatrix_gpu(idx, colIndices, values, val_size);
}else{
addMatrix_cpu(idx, colIndices, values, val_size);
}
#else
addMatrix_cpu(idx, colIndices, values, val_size);
#endif
if(idx == 2){
num_blocks_so_far += val_size;
}
if(idx == 2){
num_wells_so_far++;
}
}
void WellContributions::addMatrix_cpu(int idx, int *colIndices, double *values, unsigned int val_size)
{
switch (idx) {
case 0:
memcpy(Cnnzs + num_blocks_so_far * dim * dim_wells, values, sizeof(double) * val_size * dim * dim_wells);
memcpy(Ccols + num_blocks_so_far, colIndices, sizeof(int) * val_size);
break;
case 1:
memcpy(Dnnzs + num_wells_so_far * dim_wells * dim_wells, values, sizeof(double) * dim_wells * dim_wells);
break;
case 2:
memcpy(Bnnzs + num_blocks_so_far * dim * dim_wells, values, sizeof(double) * val_size * dim * dim_wells);
memcpy(Bcols + num_blocks_so_far, colIndices, sizeof(int) * val_size);
val_pointers[num_wells_so_far] = num_blocks_so_far;
if(num_wells_so_far == num_wells - 1){
val_pointers[num_wells] = num_blocks;
}
break;
case 3:
// store (C*D*B)
printf("ERROR unsupported matrix ID for WellContributions::addMatrix()\n");
exit(1);
// memcpy(Mnnzs, values, sizeof(double) * dim * dim * val_size * val_size);
// memcpy(Ccols, colIndices, sizeof(int) * val_size);
break;
default:
printf("ERROR unknown matrix ID for WellContributions::addMatrix()\n");
exit(1);
}
}
#if HAVE_CUDA
void WellContributions::addMatrix_gpu(int idx, int *colIndices, double *values, unsigned int val_size)
{
switch (idx) {
case 0:
cudaMemcpy(d_Cnnzs + num_blocks_so_far * dim * dim_wells, values, sizeof(double) * val_size * dim * dim_wells, cudaMemcpyHostToDevice);
@@ -387,6 +192,12 @@ namespace Opm
exit(1);
}
cudaCheckLastError("WellContributions::addMatrix() failed");
if(idx == 2){
num_blocks_so_far += val_size;
}
if(idx == 2){
num_wells_so_far++;
}
}
void WellContributions::setCudaStream(cudaStream_t stream_)
@@ -394,7 +205,6 @@ namespace Opm
this->stream = stream_;
}
#endif
void WellContributions::addSizes(unsigned int nnz, unsigned int numEq, unsigned int numWellEq)
{
@@ -408,13 +218,5 @@ namespace Opm
num_wells++;
}
// Default value
bool WellContributions::gpu_mode = false;
// If HAVE_CUDA is false, use_gpu must be false too
void WellContributions::setMode(bool use_gpu){
gpu_mode = use_gpu;
}
} //namespace Opm
opm/simulators/linalg/bda/WellContributions.hpp
+1 -29
View File
@@ -37,7 +37,6 @@ namespace Opm
{
private:
static bool gpu_mode; // gpu_mode should be initialized in the ISTLSolverEbos constructor, and its value must not change afterwards
unsigned int num_blocks = 0; // total number of blocks in all wells
unsigned int dim;
unsigned int dim_wells;
@@ -68,33 +67,6 @@ namespace Opm
double *z1 = nullptr; // z1 = B * x
double *z2 = nullptr; // z2 = D^-1 * B * x
/// Apply the WellContributions on CPU
void apply_cpu(double *x, double *y);
/// Allocate memory on the CPU
void alloc_cpu();
/// Free memory on the CPU
void free_cpu();
/// Same as addMatrix(), stores matrix on CPU
void addMatrix_cpu(int idx, int *colIndices, double *values, unsigned int val_size);
#if HAVE_CUDA
/// Apply all wellcontributions on GPU, performs y -= C^T * (D^-1 * (B * x))
/// Kernel is asynchronous
void apply_gpu(double *d_x, double *d_y);
/// Allocate memory on the GPU
void alloc_gpu();
/// Free memory on the GPU
void free_gpu();
/// Same as addMatrix(), stores matrix on GPU
void addMatrix_gpu(int idx, int *colIndices, double *values, unsigned int val_size);
#endif
public:
#if HAVE_CUDA
/// Set a cudaStream to be used
@@ -112,7 +84,7 @@ namespace Opm
void apply(double *x, double *y);
/// Allocate memory for the wellcontributions
void alloc_all();
void alloc();
/// Indicate how large the next wellcontributions are, this function cannot be called after alloc_all() is called
void addSizes(unsigned int nnz, unsigned int numEq, unsigned int numWellEq);
opm/simulators/wells/BlackoilWellModel_impl.hpp
+1 -1
View File
@@ -931,7 +931,7 @@ namespace Opm {
derived->getWellSizes(nnz, numEq, numWellEq);
wellContribs.addSizes(nnz, numEq, numWellEq);
}
wellContribs.alloc_all();
wellContribs.alloc();
for(unsigned int i = 0; i < well_container_.size(); i++){
auto& well = well_container_[i];
std::shared_ptr<StandardWell<TypeTag> > derived = std::dynamic_pointer_cast<StandardWell<TypeTag> >(well);