OpenCAEPoro  v0.5.0
A simulator for multicomponent porous media flow
Decoupling.cpp
Go to the documentation of this file.
1 
12 #include "DenseMat.hpp"
13 #include "FaspSolver.hpp"
14 
15 // ABF decoupling method
16 static void decouple_abf(dBSRmat* A, REAL* diaginv, dBSRmat* B)
17 {
18  // members of A
19  const INT ROW = A->ROW;
20  const INT NNZ = A->NNZ;
21  const INT nb = A->nb;
22  const INT nb2 = nb * nb;
23  const INT* IA = A->IA;
24  const INT* JA = A->JA;
25  REAL* val = A->val;
26 
27  INT i, k, m;
28 
29  // Create a link to dBSRmat 'B'
30  INT* IAb = B->IA;
31  INT* JAb = B->JA;
32  REAL* valb = B->val;
33  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
34  memcpy(JAb, JA, NNZ * sizeof(INT));
35 
36  for (i = 0; i < ROW; ++i) {
37  // get the diagonal sub-blocks D and their inverse
38  for (k = IA[i]; k < IA[i + 1]; ++k) {
39  if (JA[k] == i) {
40  memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
41  fasp_smat_inv(diaginv + i * nb2, nb);
42  }
43  }
44  // compute D^{-1}*A
45  for (k = IA[i]; k < IA[i + 1]; ++k) {
46  m = k * nb2;
47  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
48  }
49  } // end of main loop
50 }
51 
52 // Analytical decoupling method
53 static void decouple_anl(dBSRmat* A, REAL* diaginv, dBSRmat* B)
54 {
55  // members of A
56  const INT ROW = A->ROW;
57  const INT NNZ = A->NNZ;
58  const INT nb = A->nb;
59  const INT nb2 = nb * nb;
60  const INT* IA = A->IA;
61  const INT* JA = A->JA;
62  REAL* val = A->val;
63 
64  INT i, k, m;
65 
66  // Create a dBSRmat 'B'
67  INT* IAb = B->IA;
68  INT* JAb = B->JA;
69  REAL* valb = B->val;
70  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
71  memcpy(JAb, JA, NNZ * sizeof(INT));
72 
73  for (i = 0; i < ROW; ++i) {
74  // form the diagonal sub-blocks for analytical decoupling
75  for (k = IA[i]; k < IA[i + 1]; ++k) {
76  if (JA[k] == i) {
77  m = k * nb2;
78  fasp_smat_identity(diaginv + i * nb2, nb, nb2);
79  for (int l = 0; l < nb - 1; l++)
80  diaginv[i * nb2 + 1 + l] = -val[m + 1 + l];
81  }
82  }
83 
84  // compute D^{-1}*A
85  for (k = IA[i]; k < IA[i + 1]; ++k) {
86  m = k * nb2;
87  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
88  }
89  } // end of main loop
90 }
91 
92 static void decouple_truetrans(dBSRmat* A, REAL* diaginv, dBSRmat* B)
93 {
94  // members of A
95  const INT ROW = A->ROW;
96  const INT NNZ = A->NNZ;
97  const INT nb = A->nb;
98  const INT nb2 = nb * nb;
99  const INT* IA = A->IA;
100  const INT* JA = A->JA;
101  REAL* val = A->val;
102 
103  INT i, k, l, m;
104 
105  // Create a dBSRmat 'B'
106  INT* IAb = B->IA;
107  INT* JAb = B->JA;
108  REAL* valb = B->val;
109  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
110  memcpy(JAb, JA, NNZ * sizeof(INT));
111 
112  REAL *mat1, *mat2;
113  mat1 = new REAL[nb2];
114  mat2 = new REAL[nb2];
115 
116  // Dset0(nb2*ROW, diaginv);
117  for (i = 0; i < ROW; ++i) {
118  double Tt = 0.0;
119  // get the diagonal sub-blocks
120  // mat2 is OCP_TRUE-IMPES matrix.
121  fasp_smat_identity(mat2, nb, nb2);
122  // mat1 is the mobility ratio matrix
123  fasp_smat_identity(mat1, nb, nb2);
124 
125  for (l = 1; l < nb; ++l) {
126  mat2[l] = diaginv[i * nb2 + l];
127  Tt += diaginv[i * nb2 + l] * diaginv[i * nb2 + l * nb];
128  }
129  for (k = IA[i]; k < IA[i + 1]; ++k) {
130  if (JA[k] == i) break;
131  }
132  for (l = 0; l < nb - 1; ++l) {
133  if (val[k * nb2 + (l + 1) * nb] > 0)
134  mat1[(l + 1) * nb] = -diaginv[i * nb2 + (l + 1) * nb] / Tt; // diag(l)
135  }
136  DaABpbC(nb, nb, nb, 1, mat1, mat2, 0, diaginv + i * nb2);
137 
138  // compute D^{-1}*A
139  for (k = IA[i]; k < IA[i + 1]; ++k) {
140  m = k * nb2;
141  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
142  }
143  } // end of main loop
144 
145  delete mat1;
146  delete mat2;
147 }
148 
149 static void decouple_truetrans_alg(dBSRmat* A, REAL* diaginv, dBSRmat* B)
150 {
151  // members of A
152  const INT ROW = A->ROW;
153  const INT NNZ = A->NNZ;
154  const INT nb = A->nb;
155  const INT nb2 = nb * nb;
156  const INT* IA = A->IA;
157  const INT* JA = A->JA;
158  REAL* val = A->val;
159 
160  INT i, k, l, m;
161 
162  // Create a link to dBSRmat 'B'
163  INT* IAb = B->IA;
164  INT* JAb = B->JA;
165  REAL* valb = B->val;
166  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
167  memcpy(JAb, JA, NNZ * sizeof(INT));
168 
169  REAL *mat1, *mat2;
170  mat1 = new REAL[nb2];
171  mat2 = new REAL[nb2];
172 
173  // Dset0(nb2*ROW, diaginv);
174  for (i = 0; i < ROW; ++i) {
175  for (k = IA[i]; k < IA[i + 1]; ++k) {
176  double Tt = 0.0;
177  if (JA[k] == i) {
178  m = k * nb2;
179  fasp_smat_identity(mat2, nb, nb2);
180  for (l = 1; l < nb; ++l) {
181  mat2[l] = -val[m + l];
182  if (val[m + l * nb] > 0) {
183  Tt += -val[m + l] * val[m + l * nb];
184  }
185  }
186 
187  fasp_smat_identity(mat1, nb, nb2);
188  for (l = 1; l < nb; ++l) {
189  if (val[m + l * nb] > 0) {
190  mat1[l * nb] = -val[m + l * nb] / Tt; // diag(l)
191  }
192  }
193 
194  DaABpbC(nb, nb, nb, 1, mat1, mat2, 0, diaginv + i * nb2);
195  }
196  }
197 
198  // compute D^{-1}*A
199  for (k = IA[i]; k < IA[i + 1]; ++k) {
200  m = k * nb2;
201  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
202  }
203  } // end of main loop
204 
205  delete mat1;
206  delete mat2;
207 }
208 
209 // Semi-analytical decoupling method
210 static void decouple_abftrue(dBSRmat* A, REAL* diaginv, dBSRmat* B)
211 {
212  // members of A
213  const INT ROW = A->ROW;
214  const INT NNZ = A->NNZ;
215  const INT nb = A->nb;
216  const INT nb2 = nb * nb;
217  const INT* IA = A->IA;
218  const INT* JA = A->JA;
219  REAL* val = A->val;
220 
221  INT i, k, m;
222 
223  // Create a dBSRmat 'B'
224  INT* IAb = B->IA;
225  INT* JAb = B->JA;
226  REAL* valb = B->val;
227  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
228  memcpy(JAb, JA, NNZ * sizeof(INT));
229 
230  for (i = 0; i < ROW; ++i) {
231  // get the diagonal sub-blocks
232  for (k = IA[i]; k < IA[i + 1]; ++k) {
233  if (JA[k] == i) {
234  // Form the ABF part first
235  memcpy(diaginv + i * nb2, val + k * nb2, nb2 * sizeof(REAL));
236  fasp_smat_inv(diaginv + i * nb2, nb);
237  // Replace the first line with analytical decoupling
238  m = k * nb2;
239  diaginv[i * nb2] = 1;
240  for (int l = 0; l < nb - 1; l++)
241  diaginv[i * nb2 + 1 + l] = -val[m + 1 + l];
242  }
243  }
244 
245  // compute D^{-1}*A
246  for (k = IA[i]; k < IA[i + 1]; ++k) {
247  m = k * nb2;
248  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
249  }
250  } // end of main loop
251 }
252 
253 static void decouple_true_scale(dBSRmat* A, REAL* diaginv, dBSRmat* B)
254 {
255  // members of A
256  const INT ROW = A->ROW;
257  const INT NNZ = A->NNZ;
258  const INT nb = A->nb;
259  const INT nb2 = nb * nb;
260  const INT* IA = A->IA;
261  const INT* JA = A->JA;
262  REAL* val = A->val;
263 
264  INT i, k, m;
265 
266  // Create a dBSRmat 'B'
267  INT* IAb = B->IA;
268  INT* JAb = B->JA;
269  REAL* valb = B->val;
270  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
271  memcpy(JAb, JA, NNZ * sizeof(INT));
272 
273  for (i = 0; i < ROW; ++i) {
274  // get the diagonal sub-blocks
275  for (k = IA[i]; k < IA[i + 1]; ++k) {
276  if (JA[k] == i) {
277  m = k * nb2;
278  fasp_smat_identity(diaginv + i * nb2, nb, nb2);
279  diaginv[i * nb2] = 1 / val[m];
280  for (int l = 0; l < nb - 1; l++)
281  diaginv[i * nb2 + 1 + l] = -val[m + 1 + l] / val[m];
282  }
283  }
284 
285  // compute D^{-1}*A
286  for (k = IA[i]; k < IA[i + 1]; ++k) {
287  m = k * nb2;
288  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
289  }
290  } // end of main loop
291 }
292 
293 static void decouple_rotate(dBSRmat* A, REAL* diaginv, dBSRmat* B)
294 {
295  // members of A
296  const INT ROW = A->ROW;
297  const INT NNZ = A->NNZ;
298  const INT nb = A->nb;
299  const INT nb2 = nb * nb;
300  const INT* IA = A->IA;
301  const INT* JA = A->JA;
302  REAL* val = A->val;
303 
304  INT i, k, m;
305 
306  // Create a dBSRmat 'B'
307  INT* IAb = B->IA;
308  INT* JAb = B->JA;
309  REAL* valb = B->val;
310  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
311  memcpy(JAb, JA, NNZ * sizeof(INT));
312 
313  for (i = 0; i < ROW; ++i) {
314  // get the diagonal sub-blocks
315  for (k = IA[i]; k < IA[i + 1]; ++k) {
316  if (JA[k] == i) {
317  m = k * nb2;
318  for (int li = 0; li < nb; li++) {
319  for (int lj = 0; lj < nb; lj++) {
320  diaginv[i * nb2 + li * nb + lj] = 0;
321  if (lj - li == 1) {
322  diaginv[i * nb2 + li * nb + lj] = 1;
323  }
324  if (lj == 0 && li == nb - 1) {
325  diaginv[i * nb2 + li * nb + lj] = 1;
326  }
327  }
328  }
329  }
330  }
331 
332  // compute D^{-1}*A
333  for (k = IA[i]; k < IA[i + 1]; ++k) {
334  m = k * nb2;
335  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
336  }
337  } // end of main loop
338 }
339 
340 // Quasi-IMPES decoupling method
341 static void decouple_quasi(dBSRmat* A, REAL* diaginv, dBSRmat* B)
342 {
343  // members of A
344  const INT ROW = A->ROW;
345  const INT NNZ = A->NNZ;
346  const INT nb = A->nb;
347  const INT nb2 = nb * nb;
348  const INT* IA = A->IA;
349  const INT* JA = A->JA;
350  REAL* val = A->val;
351 
352  INT i, k, m;
353 
354  INT* IAb = B->IA;
355  INT* JAb = B->JA;
356  REAL* valb = B->val;
357  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
358  memcpy(JAb, JA, NNZ * sizeof(INT));
359 
360  for (i = 0; i < ROW; ++i) {
361  // get the diagonal sub-blocks
362  for (k = IA[i]; k < IA[i + 1]; ++k) {
363  if (JA[k] == i) {
364  m = k * nb2;
365  fasp_smat_identity(diaginv + i * nb2, nb, nb2);
366  for (int l = 0; l < nb - 1; l++)
367  diaginv[i * nb2 + 1 + l] =
368  -val[m + 1 + l] / val[m + (l + 1) * nb + l + 1];
369  }
370  }
371 
372  // compute D^{-1}*A
373  for (k = IA[i]; k < IA[i + 1]; ++k) {
374  m = k * nb2;
375  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
376  }
377  } // end of main loop
378 }
379 
380 // What's the difference with decouple_anl?
381 static void decouple_trueabf(dBSRmat* A, REAL* diaginv, dBSRmat* B)
382 {
383  // members of A
384  const INT ROW = A->ROW;
385  const INT NNZ = A->NNZ;
386  const INT nb = A->nb;
387  const INT nb2 = nb * nb;
388  const INT* IA = A->IA;
389  const INT* JA = A->JA;
390  REAL* val = A->val;
391 
392  INT i, k, m;
393 
394  INT* IAb = B->IA;
395  INT* JAb = B->JA;
396  REAL* valb = B->val;
397  memcpy(IAb, IA, (ROW + 1) * sizeof(INT));
398  memcpy(JAb, JA, NNZ * sizeof(INT));
399 
400  for (i = 0; i < ROW; ++i) {
401  // get the diagonal sub-blocks
402  for (k = IA[i]; k < IA[i + 1]; ++k) {
403  if (JA[k] == i) {
404  m = k * nb2;
405  fasp_smat_identity(diaginv + i * nb2, nb, nb2);
406  for (int l = 0; l < nb; l++) diaginv[i * nb2 + l] = val[m + l];
407  fasp_smat_inv(diaginv + i * nb2, nb);
408  }
409  }
410 
411  // compute D^{-1}*A
412  for (k = IA[i]; k < IA[i + 1]; ++k) {
413  m = k * nb2;
414  fasp_blas_smat_mul(diaginv + i * nb2, val + m, valb + m, nb);
415  }
416  } // end of main loop
417 }
418 
420 void VectorFaspSolver::Decoupling(dBSRmat* Absr,
421  dvector* b,
422  dBSRmat* Asc,
423  dvector* fsc,
424  ivector* order,
425  double* Dmatvec,
426  int decoupleType)
427 {
428  int nrow = Absr->ROW;
429  int nb = Absr->nb;
430  double* Dmat = Dmatvec;
431  precond_diag_bsr diagA;
432 
433  // Natural ordering
434  for (int i = 0; i < nrow; ++i) order->val[i] = i;
435 
436  // Without decoupling
437  if (decoupleType == 0) {
438  fasp_dbsr_cp(Absr, Asc); // Asc = Absr;
439  fasp_dvec_cp(b, fsc); // fsc = b;
440  return;
441  }
442 
443  // With decoupling
444  switch (decoupleType) {
445  case 2:
446  decouple_anl(Absr, Dmat, Asc);
447  break;
448  case 3:
449  decouple_quasi(Absr, Dmat, Asc);
450  break;
451  case 4:
452  decouple_trueabf(Absr, Dmat, Asc);
453  break;
454  case 5:
455  decouple_abftrue(Absr, Dmat, Asc);
456  break;
457  case 6:
458  decouple_abftrue(Absr, Dmat, Asc);
459  break;
460  case 7:
461  decouple_truetrans_alg(Absr, Dmat, Asc);
462  break;
463  case 8:
464  decouple_truetrans(Absr, Dmat, Asc);
465  break;
466  case 9:
467  decouple_true_scale(Absr, Dmat, Asc);
468  break;
469  case 10:
470  decouple_rotate(Absr, Dmat, Asc);
471  break;
472  default: // case 1:
473  decouple_abf(Absr, Dmat, Asc);
474  }
475 
476  diagA.diag.row = nrow * nb * nb;
477  diagA.diag.val = Dmat;
478  diagA.nb = nb;
479  fasp_precond_dbsr_diag(b->val, fsc->val, &diagA);
480 }
481 
482 /*----------------------------------------------------------------------------*/
483 /* Brief Change History of This File */
484 /*----------------------------------------------------------------------------*/
485 /* Author Date Actions */
486 /*----------------------------------------------------------------------------*/
487 /* Shizhe Li Oct/15/2021 Create file */
488 /* Chensong Zhang Nov/09/2021 Restruct decoupling methods */
489 /* Chensong Zhang Nov/30/2021 Add null decoupling */
490 /*----------------------------------------------------------------------------*/
DenseMat.hpp
Operations about small dense mat.
DaABpbC
void DaABpbC(const int &m, const int &n, const int &k, const double &alpha, const double *A, const double *B, const double &beta, double *C)
Computes C' = alpha B'A' + beta C', all matrices are column-major.
Definition: DenseMat.cpp:75
FaspSolver.hpp
Declaration of classes interfacing to the FASP solvers.
Generated on Tue Feb 7 2023 01:49:50 for OpenCAEPoro by doxygen 1.9.1