|
60 | 60 | #include <unordered_map> |
61 | 61 | #include <vector> |
62 | 62 |
|
| 63 | +#include <fmt/format.h> |
| 64 | + |
63 | 65 | #include <opm/common/utility/gpuDecorators.hpp> |
64 | 66 | #if HAVE_CUDA |
65 | 67 | #if USE_HIP |
@@ -236,6 +238,8 @@ public: |
236 | 238 | { |
237 | 239 | simulatorPtr_ = nullptr; |
238 | 240 | separateSparseSourceTerms_ = Parameters::Get<Parameters::SeparateSparseSourceTerms>(); |
| 241 | + exportIndex_=-1; |
| 242 | + exportCount_=-1; |
239 | 243 | } |
240 | 244 |
|
241 | 245 | /*! |
@@ -414,6 +418,201 @@ public: |
414 | 418 | GlobalEqVector& residual() |
415 | 419 | { return residual_; } |
416 | 420 |
|
| 421 | + /*! |
| 422 | + * \brief Print first 16 block elements of a vector. |
| 423 | + */ |
| 424 | + void printVector(GlobalEqVector &x, const char *name="x") |
| 425 | + { |
| 426 | + int count = 1; |
| 427 | + fmt::print("{} =\n[\n",name); |
| 428 | + for (auto block = x.begin(); block != x.end(); block++) |
| 429 | + { |
| 430 | + for (auto i=block->begin(); i!=block->end(); i++) |
| 431 | + { |
| 432 | + fmt::print(" {:+.4e}",*i); |
| 433 | + } |
| 434 | + fmt::print("\n"); |
| 435 | + count++; |
| 436 | + if(count>16) break; |
| 437 | + } |
| 438 | + fmt::print("]\n"); |
| 439 | + } |
| 440 | + |
| 441 | + /*! |
| 442 | + * \brief Print first 16 block elements of residual. |
| 443 | + */ |
| 444 | + void printResidual(const char *name="r") |
| 445 | + { |
| 446 | + printVector(residual_, name); |
| 447 | + } |
| 448 | + |
| 449 | + /*! |
| 450 | + * \brief Print sparsity pattern of the first 16 rows |
| 451 | + * of the jacobian block matrix |
| 452 | + */ |
| 453 | + void printSparsity(const char *name="s") |
| 454 | + { |
| 455 | + auto& A = jacobian_->istlMatrix(); |
| 456 | + |
| 457 | + fmt::print("nrows = {:d}\n",A.N()); |
| 458 | + fmt::print("ncols = {:d}\n",A.M()); |
| 459 | + fmt::print("nnz = {:d}\n",A.nonzeroes()); |
| 460 | + |
| 461 | + fmt::print("{} =\n[\n",name); |
| 462 | + int count=1; |
| 463 | + int offset=0; |
| 464 | + for(auto row=A.begin(); row!=A.end(); row++) |
| 465 | + { |
| 466 | + fmt::print("{:4d}: ",offset); |
| 467 | + for(unsigned int i=0;i<row->getsize();i++) |
| 468 | + { |
| 469 | + fmt::print(" {:4d}",row->getindexptr()[i]); |
| 470 | + } |
| 471 | + fmt::print("\n"); |
| 472 | + offset+=row->getsize(); |
| 473 | + count++; |
| 474 | + if(count>16) break; |
| 475 | + } |
| 476 | + fmt::print("]\n"); |
| 477 | + } |
| 478 | + |
| 479 | + /*! |
| 480 | + * \brief Print block elements of the first 6 rows of the |
| 481 | + * j acobian block matrix |
| 482 | + */ |
| 483 | + void printNonzeros(const char *name="d") |
| 484 | + { |
| 485 | + auto& A = jacobian_->istlMatrix(); |
| 486 | + fmt::print("{} =\n[\n",name); |
| 487 | + int count=1; |
| 488 | + for(auto row=A.begin();row!=A.end();row++) |
| 489 | + { |
| 490 | + for(unsigned int j=0;j<row->getsize();j++) |
| 491 | + { |
| 492 | + fmt::print("|"); |
| 493 | + auto mat = row->getptr()[j]; |
| 494 | + for(auto vec=mat.begin();vec!=mat.end();vec++) |
| 495 | + { |
| 496 | + for(auto k=vec->begin();k!=vec->end();k++) |
| 497 | + { |
| 498 | + fmt::print(" {:+.4e}",*k); |
| 499 | + } |
| 500 | + fmt::print(" |"); |
| 501 | + } |
| 502 | + fmt::print("\n"); |
| 503 | + } |
| 504 | + count++; |
| 505 | + if(count>6) break; |
| 506 | + fmt::print("\n"); |
| 507 | + } |
| 508 | + fmt::print("]\n"); |
| 509 | + } |
| 510 | + |
| 511 | + /*! |
| 512 | + * \brief Print sparsity pattern and nonzeros of jacobian block matrix |
| 513 | + */ |
| 514 | + void printJacobian() |
| 515 | + { |
| 516 | + printSparsity(); |
| 517 | + printNonzeros(); |
| 518 | + } |
| 519 | + |
| 520 | + /*! |
| 521 | + * \brief Export blocks-sparse linear system. |
| 522 | + */ |
| 523 | + void exportSystem(int idx, char *tag, const char *path="export") |
| 524 | + { |
| 525 | + // export sparsity only once |
| 526 | + if(exportIndex_==-1) exportSparsity(path); |
| 527 | + |
| 528 | + // increment indices and generate tag |
| 529 | + exportCount_ = exportIndex_==idx ? ++exportCount_ : 0; |
| 530 | + exportIndex_ = idx; |
| 531 | + sprintf(tag,"_%03d_%02d",exportIndex_, exportCount_); |
| 532 | + |
| 533 | + fmt::print("index = {:d}\n", exportIndex_); |
| 534 | + fmt::print("count = {:d}\n", exportCount_); |
| 535 | + |
| 536 | + // export matrix |
| 537 | + exportNonzeros(tag,path); |
| 538 | + |
| 539 | + // export residual |
| 540 | + char name[256]; |
| 541 | + sprintf(name,"%s/r",path); |
| 542 | + exportVector(residual_,tag,name); |
| 543 | + } |
| 544 | + |
| 545 | + /*! |
| 546 | + * \brief Export block vector. |
| 547 | + */ |
| 548 | + void exportVector(GlobalEqVector &x, const char *tag="", const char *name="export/x") |
| 549 | + { |
| 550 | + // assume double precision and contiguous data |
| 551 | + const double *data = &x[0][0]; |
| 552 | + |
| 553 | + char filename[512]; |
| 554 | + sprintf(filename,"%s%s.f64",name,tag); |
| 555 | + FILE *out =fopen(filename,"w"); |
| 556 | + fwrite(data, sizeof(double), x.dim(),out); |
| 557 | + fclose(out); |
| 558 | + } |
| 559 | + |
| 560 | + /*! |
| 561 | + * \brief Export nonzero blocks of jacobian block-sparse matrix |
| 562 | + */ |
| 563 | + void exportNonzeros(const char *tag="", const char *path=".") |
| 564 | + { |
| 565 | + auto& A = jacobian_->istlMatrix(); |
| 566 | + |
| 567 | + // assume double precision and contiguous data |
| 568 | + const double *data = &A[0][0][0][0]; |
| 569 | + size_t dim = A[0][0].N()*A[0][0].M()*A.nonzeroes(); |
| 570 | + |
| 571 | + char filename[256]; |
| 572 | + sprintf(filename,"%s/data%s.f64",path,tag); |
| 573 | + FILE *out =fopen(filename,"w"); |
| 574 | + fwrite(data, sizeof(double), dim,out); |
| 575 | + fclose(out); |
| 576 | + } |
| 577 | + |
| 578 | + /*! |
| 579 | + * \brief Export sparsity pattern of jacobian block-sparse matrix |
| 580 | + */ |
| 581 | + void exportSparsity(const char *path=".") |
| 582 | + { |
| 583 | + //assemble csr graph |
| 584 | + auto& A = jacobian_->istlMatrix(); |
| 585 | + auto rows = std::make_unique<int[]>(A.N()+1); |
| 586 | + auto cols = std::make_unique<int[]>(A.nonzeroes()); |
| 587 | + |
| 588 | + int irow=0; |
| 589 | + int icol=0; |
| 590 | + rows[0]=0; |
| 591 | + for(auto row=A.begin(); row!=A.end(); row++) |
| 592 | + { |
| 593 | + for(unsigned int i=0;i<row->getsize();i++) |
| 594 | + { |
| 595 | + cols[icol++]=row->getindexptr()[i]; |
| 596 | + } |
| 597 | + rows[irow+1]= rows[irow]+row->getsize(); |
| 598 | + irow++; |
| 599 | + } |
| 600 | + |
| 601 | + //export arrays |
| 602 | + FILE *out; |
| 603 | + char filename[256]; |
| 604 | + |
| 605 | + sprintf(filename,"%s/rows.i32",path); |
| 606 | + out=fopen(filename,"w"); |
| 607 | + fwrite(rows, sizeof(int), A.N()+1,out); |
| 608 | + fclose(out); |
| 609 | + |
| 610 | + sprintf(filename,"%s/cols.i32",path); |
| 611 | + out=fopen(filename,"w"); |
| 612 | + fwrite(cols, sizeof(int), A.nonzeroes(),out); |
| 613 | + fclose(out); |
| 614 | + } |
| 615 | + |
417 | 616 | void setLinearizationType(LinearizationType linearizationType) |
418 | 617 | { linearizationType_ = linearizationType; } |
419 | 618 |
|
@@ -1067,6 +1266,9 @@ private: |
1067 | 1266 | bool separateSparseSourceTerms_ = false; |
1068 | 1267 |
|
1069 | 1268 | FullDomain<> fullDomain_; |
| 1269 | + |
| 1270 | + int exportIndex_; |
| 1271 | + int exportCount_; |
1070 | 1272 | }; |
1071 | 1273 | } // namespace Opm |
1072 | 1274 |
|
|
0 commit comments