Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/algebraic/hazmath/precond.py: 53%

1from block.algebraic import petsc # NOTE: Not used, import initializes petsc from cmdline

2from block.block_base import block_base

3from block import block_mat, supports_mpi

4from builtins import str

5from petsc4py import PETSc

6from scipy.sparse import csr_matrix

7import dolfin as df

8import numpy as np

9import haznics

11# ------------------------------------------------------------------- #

12# -------------- auxiliary functions --------------- #

13# ------------------------------------------------------------------- #

16def PETSc_to_dCSRmat(A):

17 """

18 Change data type for matrix

19 (dolfin PETScMatrix or GenericMatrix to dCSRmat pointer)

20 """

21 if not isinstance(A, df.PETScMatrix): 21 ↛ 23line 21 didn't jump to line 23, because the condition on line 21 was never false

22 A = df.as_backend_type(A)

23 petsc_mat = A.mat()

25 # NB! store copies for now

26 csr = petsc_mat.getValuesCSR()

28 return haznics.create_matrix(csr[2], csr[1], csr[0], A.size(1))

31def block_mat_to_block_dCSRmat(A):

32 """

33 Change data type for matrix

34 (block.block_mat to block_dCSRmat pointer)

35 """

36 # check type

37 assert isinstance(A, block_mat)

39 # allocate block_dCSRmat and its blocks too

40 brow, bcol = A.blocks.shape

41 Abdcsr = haznics.block_dCSRmat()

42 Abdcsr.init(brow, bcol)

44 for i in range(brow):

45 for j in range(bcol):

46 if isinstance(A[i][j], df.Matrix):

47 A[i][j] = df.as_backend_type(A[i][j])

48 if isinstance(A[i][j], df.PETScMatrix):

49 csr = A[i][j].mat().getValuesCSR() # todo: eliminate zeros!

50 mat = haznics.create_matrix(csr[2], csr[1], csr[0], A[i][j].size(1))

51 Abdcsr.set(i, j, mat)

52 elif not A[i][j]:

53 mat = haznics.dCSRmat()

54 haznics.dcsr_null(mat)

55 Abdcsr.set(i, j, mat)

56 else:

57 return NotImplemented

59 return Abdcsr

61# copied from block/algebraic/petsc/

62def discrete_gradient(mesh):

63 """ P1 to Ned1 map """

64 Ned = df.FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1)

65 P1 = df.FunctionSpace(mesh, 'CG', 1)

67 return df.as_backend_type(df.DiscreteOperators.build_gradient(Ned, P1))

70def discrete_curl(mesh):

71 """ Ned1 to RT1 map """

72 assert mesh.geometry().dim() == 3

73 assert mesh.topology().dim() == 3

75 Ned = df.FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1)

76 RT = df.FunctionSpace(mesh, 'RT', 1)

78 RT_f2dof = np.array(RT.dofmap().entity_dofs(mesh, 2))

79 Ned_e2dof = np.array(Ned.dofmap().entity_dofs(mesh, 1))

81 # Facets in terms of edges

82 mesh.init(2, 1)

83 f2e = mesh.topology()(2, 1)

85 rows = np.repeat(RT_f2dof, 3)

86 cols = Ned_e2dof[f2e()]

87 vals = np.tile(np.array([-2, 2, -2]), RT.dim())

89 Ccsr = csr_matrix((vals, (rows, cols)), shape=(RT.dim(), Ned.dim()))

91 C = PETSc.Mat().createAIJ(comm=df.MPI.comm_world,

92 size=Ccsr.shape,

93 csr=(Ccsr.indptr, Ccsr.indices, Ccsr.data))

95 return df.PETScMatrix(C)

98# NB: this is only 3d!

99def Pdiv(mesh):

100 """ nodes to faces map """

101 RT = df.FunctionSpace(mesh, 'Raviart-Thomas', 1)

102 P1 = df.VectorFunctionSpace(mesh, 'CG', 1)

103

104 RT_f2dof = np.array(RT.dofmap().entity_dofs(mesh, 2))

105 P1_n2dof_x = np.array(P1.sub(0).dofmap().entity_dofs(mesh, 0))

106 P1_n2dof_y = np.array(P1.sub(1).dofmap().entity_dofs(mesh, 0))

107 P1_n2dof_z = np.array(P1.sub(2).dofmap().entity_dofs(mesh, 0))

108

109 # Facets in terms of nodes

110 mesh.init(2, 0)

111 f2n = mesh.topology()(2, 0)

112 coordinates = mesh.coordinates()

113

114 rows_x, rows_y, rows_z = np.repeat(RT_f2dof, 3), np.repeat(RT_f2dof, 3), np.repeat(RT_f2dof, 3)

115 cols_x, cols_y, cols_z = P1_n2dof_x[f2n()], P1_n2dof_y[f2n()], P1_n2dof_z[f2n()]

116 vals_x, vals_y, vals_z = np.zeros(3 * RT.dim()), np.zeros(3 * RT.dim()), np.zeros(3 * RT.dim())

117

118 for facet, row in enumerate(RT_f2dof):

119 vertices = f2n(facet)

120 n1, n2, n3 = coordinates[vertices]

121

122 facet_normal = np.cross(n3 - n1, n2 - n1)

123

124 indices = np.arange(3) + 3 * facet

125 vals_x[indices] = facet_normal[0] / 3

126 vals_y[indices] = facet_normal[1] / 3

127 vals_z[indices] = facet_normal[2] / 3

128

129 rows = np.concatenate((rows_x, rows_y, rows_z))

130 cols = np.concatenate((cols_x, cols_y, cols_z))

131 vals = np.concatenate((vals_x, vals_y, vals_z))

132

133 # assemble Pdiv as from xyz components

134 Pdivcsr = csr_matrix((vals, (rows, cols)), shape=(RT.dim(), P1.dim()))

135 Pdivcsr.eliminate_zeros()

136

137 Pdiv = PETSc.Mat().createAIJ(comm=df.MPI.comm_world,

138 size=Pdivcsr.shape,

139 csr=(Pdivcsr.indptr, Pdivcsr.indices, Pdivcsr.data))

140

141 return df.PETScMatrix(Pdiv)

142

143

144def Pcurl(mesh):

145 """ nodes to edges map """

146 assert mesh.geometry().dim() == 3

147 assert mesh.topology().dim() == 3

148

149 Ned = df.FunctionSpace(mesh, 'Nedelec 1st kind H(curl)', 1)

150 P1 = df.VectorFunctionSpace(mesh, 'CG', 1)

151

152 Ned_e2dof = np.array(Ned.dofmap().entity_dofs(mesh, 1))

153 P1_n2dof_x = np.array(P1.sub(0).dofmap().entity_dofs(mesh, 0))

154 P1_n2dof_y = np.array(P1.sub(1).dofmap().entity_dofs(mesh, 0))

155 P1_n2dof_z = np.array(P1.sub(2).dofmap().entity_dofs(mesh, 0))

156

157 # Facets in terms of edges

158 mesh.init(1, 0)

159 e2n = mesh.topology()(1, 0)

160 coordinates = mesh.coordinates()

161

162 rows_x, rows_y, rows_z = np.repeat(Ned_e2dof, 2), np.repeat(Ned_e2dof, 2), np.repeat(Ned_e2dof, 2)

163 cols_x, cols_y, cols_z = P1_n2dof_x[e2n()], P1_n2dof_y[e2n()], P1_n2dof_z[e2n()]

164 vals_x, vals_y, vals_z = np.zeros(2 * Ned.dim()), np.zeros(2 * Ned.dim()), np.zeros(2 * Ned.dim())

165

166 for edge, row in enumerate(Ned_e2dof):

167 vertices = e2n(edge)

168 edge_tangent = coordinates[vertices[1]] - coordinates[vertices[0]]

169

170 indices = np.arange(2) + 2 * edge

171 vals_x[indices] = 0.5 * edge_tangent[0]

172 vals_y[indices] = 0.5 * edge_tangent[1]

173 vals_z[indices] = 0.5 * edge_tangent[2]

174

175 rows = np.concatenate((rows_x, rows_y, rows_z))

176 cols = np.concatenate((cols_x, cols_y, cols_z))

177 vals = np.concatenate((vals_x, vals_y, vals_z))

178

179 # assemble Pcurl from xyz components

180 Pcurlcsr = csr_matrix((vals, (rows, cols)), shape=(Ned.dim(), P1.dim()))

181 Pcurlcsr.eliminate_zeros()

182

183 Pcurl = PETSc.Mat().createAIJ(comm=df.MPI.comm_world,

184 size=Pcurlcsr.shape,

185 csr=(Pcurlcsr.indptr, Pcurlcsr.indices, Pcurlcsr.data))

186

187 return df.PETScMatrix(Pcurl)

188

189

190# ------------------------------------------------------------------- #

191# -------------- preconditioners --------------- #

192# ------------------------------------------------------------------- #

193

194

195class Precond(block_base):

196 """

197 Class of general preconditioners from HAZmath using SWIG

198

199 Note: Parallel execution (MPI) not supported.

200 """

201

202 def __init__(self, A, prectype=None, parameters=None, amg_parameters=None, precond=None):

203 supports_mpi(False, 'HAZMath does not work in parallel')

204

205 # haznics.dCSRmat* type (assert?)

206 self.A = A

207 # python dictionary of parameters

208 self.__parameters = parameters if parameters else {}

209 self.__amg_parameters = amg_parameters

210

211 # init and set preconditioner (precond *)

212 if precond: 212 ↛ 215line 212 didn't jump to line 215, because the condition on line 212 was never false

213 self.precond = precond

214 else:

215 import warnings

216 warnings.warn(

217 "!! Preconditioner not specified !! Creating default UA-AMG "

218 "precond... ",

219 RuntimeWarning)

220 # change data type for the matrix (to dCSRmat pointer)

221 A_ptr = PETSc_to_dCSRmat(A)

222

223 # initialize amg parameters (AMG_param pointer)

224 amgparam = haznics.AMG_param()

225

226 # set extra amg parameters

227 if parameters:

228 haznics.param_amg_set_dict(parameters, amgparam)

229

230 # print (relevant) amg parameters

231 haznics.param_amg_print(amgparam)

232 self.__amg_parameters = amgparam

233

234 self.precond = haznics.create_precond(A_ptr, amgparam)

235

236 # if fail, setup returns null

237 if not precond:

238 raise RuntimeError(

239 "AMG levels failed to set up (null pointer returned) ")

240

241 # setup time

242 self.setup_time = precond.setup_time if precond.setup_time else 0.

243

244 # preconditioner type (string)

245 self.prectype = prectype if prectype else "AMG"

246

247 def matvec(self, b):

248 if not isinstance(b, df.GenericVector): 248 ↛ 249line 248 didn't jump to line 249, because the condition on line 248 was never true

249 return NotImplemented

250

251 x = self.A.create_vec(dim=1)

252 x = df.Vector(df.MPI.comm_self, x.size())

253

254 if len(x) != len(b): 254 ↛ 255line 254 didn't jump to line 255, because the condition on line 254 was never true

255 raise RuntimeError('incompatible dimensions for matvec, %d != %d'

256 % (len(x), len(b)))

257

258 # convert rhs and dx to numpy arrays

259 b_np = b[:]

260 x_np = x[:]

261

262 # apply the preconditioner (solution dx saved in x_np)

263 haznics.apply_precond(b_np, x_np, self.precond)

264 # convert dx to GenericVector

265 x.set_local(x_np)

266

267 return x

268

269 # noinspection PyMethodMayBeStatic

270 def down_cast(self):

271 return NotImplemented

272

273 def __str__(self):

274 return '<%s prec of %s>' % (self.__class__.__name__, str(self.A))

275

276 def create_vec(self, dim=1):

277 return self.A.create_vec(dim)

278

279 def print_amg_parameters(self):

280 haznics.param_amg_print(self.__amg_parameters)

281

282 def print_all_parameters(self):

283 if self.prectype == "AMG": 283 ↛ 284line 283 didn't jump to line 284, because the condition on line 283 was never true

284 self.print_amg_parameters()

285 elif self.prectype == "FAMG": 285 ↛ 286line 285 didn't jump to line 286, because the condition on line 285 was never true

286 self.print_amg_parameters()

287 print(" Other parameters:")

288 print("-----------------------------------------------")

289 print("Fractional exponent: ", self.__parameters['fpwr']) if self.__parameters['fpwr'] \

290 else print()

291 print("-----------------------------------------------")

292 elif self.prectype == "RA": 292 ↛ 300line 292 didn't jump to line 300, because the condition on line 292 was never false

293 print(" Other parameters:")

294 print("-----------------------------------------------")

295 print("Fractional exponents: ", self.__parameters['pwrs']) if self.__parameters['pwrs'] \

296 else print()

297 print("Fractional weights: ", self.__parameters['coefs']) if self.__parameters['coefs'] \

298 else print()

299 print("-----------------------------------------------")

300 elif self.prectype == "HXCurl_add":

301 print(" Other parameters:")

302 print("-----------------------------------------------")

303 print("HXCurl type: ", "Additive") if self.__parameters['prectype'] else print()

304 print("-----------------------------------------------")

305 elif self.prectype == "HXDiv_add":

306 print(" Other parameters:")

307 print("-----------------------------------------------")

308 print("HXDiv type: ", "Additive") if self.__parameters['prectype'] else print()

309 print("-----------------------------------------------")

310 else:

311 NotImplemented

312

313

314class AMG(Precond):

315 """

316 AMG preconditioner from the HAZmath Library with SWIG

317

318 """

319

320 def __init__(self, A, parameters=None):

321 supports_mpi(False, 'HAZMath does not work in parallel')

322

323 # change data type for the matrix (to dCSRmat pointer)

324 A_ptr = PETSc_to_dCSRmat(A)

325

326 # initialize amg parameters (AMG_param pointer)

327 amgparam = haznics.AMG_param()

328

329 # set extra amg parameters

330 if parameters:

331 haznics.param_amg_set_dict(parameters, amgparam)

332

333 # print (relevant) amg parameters

334 haznics.param_amg_print(amgparam)

335

336 # set AMG preconditioner

337 precond = haznics.create_precond_amg(A_ptr, amgparam)

338

339 # if fail, setup returns null

340 if not precond: 340 ↛ 341line 340 didn't jump to line 341, because the condition on line 340 was never true

341 raise RuntimeError(

342 "AMG levels failed to set up (null pointer returned) ")

343

344 Precond.__init__(self, A, "AMG", parameters, amgparam, precond)

345

346

347class FAMG(Precond):

348 """

349 Fractional AMG preconditioner from the HAZmath Library

350

351 """

352

353 def __init__(self, A, M, parameters=None):

354 supports_mpi(False, 'HAZMath does not work in parallel')

355

356 # change data type for the matrices (to dCSRmat pointer)

357 A_ptr = PETSc_to_dCSRmat(A)

358 M_ptr = PETSc_to_dCSRmat(M)

359

360 # initialize amg parameters (AMG_param pointer)

361 amgparam = haznics.AMG_param()

362

363 # set extra amg parameters

364 parameters = parameters if (parameters and isinstance(parameters, dict)) \

365 else {'fpwr': 0.5, 'smoother': haznics.SMOOTHER_FJACOBI}

366 haznics.param_amg_set_dict(parameters, amgparam)

367

368 # print (relevant) amg parameters

369 haznics.param_amg_print(amgparam)

370

371 # set FAMG preconditioner

372 precond = haznics.create_precond_famg(A_ptr, M_ptr, amgparam)

373

374 # if fail, setup returns null

375 if not precond:

376 raise RuntimeError(

377 "FAMG levels failed to set up (null pointer returned) ")

378

379 Precond.__init__(self, A, "FAMG", parameters, amgparam, precond)

380

381

382class RA(Precond):

383 """

384 Rational approximation preconditioner from the HAZmath library

385

386 """

387

388 def __init__(self, A, M, parameters=None):

389 supports_mpi(False, 'HAZMath does not work in parallel')

390

391 # change data type for the matrices (to dCSRmat pointer)

392 A_ptr = PETSc_to_dCSRmat(A)

393 M_ptr = PETSc_to_dCSRmat(M)

394

395 # initialize amg parameters (AMG_param pointer)

396 amgparam = haznics.AMG_param()

397

398 # set extra amg parameters

399 parameters = parameters if (parameters and isinstance(parameters, dict)) \

400 else {'coefs': [1.0, 0.0], 'pwrs': [0.5, 0.0]}

401 haznics.param_amg_set_dict(parameters, amgparam)

402

403 # print (relevant) amg parameters

404 haznics.param_amg_print(amgparam)

405

406 # get scalings

407 scaling_a = 1. / A.norm("linf")

408 scaling_m = 1. / df.as_backend_type(M).mat().getDiagonal().min()[1]

409

410 # get coefs and powers

411 alpha, beta = parameters['coefs']

412 s_power, t_power = parameters['pwrs']

413

414 # set RA preconditioner #

415 precond = haznics.create_precond_ra(A_ptr, M_ptr, s_power, t_power,

416 alpha, beta, scaling_a, scaling_m,

417 ra_tol=1e-6, amgparam=amgparam)

418

419 # if fail, setup returns null

420 if not precond: 420 ↛ 421line 420 didn't jump to line 421, because the condition on line 420 was never true

421 raise RuntimeError(

422 "Rational Approximation data failed to set up (null pointer "

423 "returned) ")

424

425 Precond.__init__(self, A, "RA", parameters, amgparam, precond)

426

427

428class HXCurl(Precond):

429 """

430 HX preconditioner from the HAZmath library for the curl-curl inner product

431 NB! only for 3D problems

432 """

433

434 def __init__(self, Acurl, V, parameters=None):

435 supports_mpi(False, 'HAZMath does not work in parallel')

436

437 # get auxiliary operators

438 mesh = V.mesh()

439 Pc = Pcurl(mesh)

440 Grad = discrete_gradient(mesh)

441

442 # change data type for the matrices (to dCSRmat pointer)

443 Acurl_ptr = PETSc_to_dCSRmat(Acurl)

444 Pcurl_ptr = PETSc_to_dCSRmat(Pc)

445 Grad_ptr = PETSc_to_dCSRmat(Grad)

446

447 # initialize amg parameters (AMG_param pointer)

448 amgparam = haznics.AMG_param()

449

450 # set extra amg parameters

451 if parameters and isinstance(parameters, dict):

452 haznics.param_amg_set_dict(parameters, amgparam)

453 else:

454 parameters = {'prectype': haznics.PREC_HX_CURL_A}

455

456 # make sure coarse solver is always iterative

457 amgparam.coarse_solver = 0

458 # print (relevant) amg parameters

459 haznics.param_amg_print(amgparam)

460

461 # add or multi

462 try:

463 prectype = parameters['prectype']

464 except KeyError:

465 prectype = haznics.PREC_HX_CURL_A

466

467 # set HX CURL preconditioner (NB: this sets up both data and fct)

468 precond = haznics.create_precond_hxcurl(Acurl_ptr, Pcurl_ptr, Grad_ptr,

469 prectype, amgparam)

470

471 # if fail, setup returns null

472 if not precond:

473 raise RuntimeError(

474 "HXcurl data failed to set up (null pointer returned) ")

475 """

476 try:

477 prectype = parameters['prectype']

478 except KeyError:

479 prectype = ''

480

481 if prectype in ["add", "Add", "ADD", "additive", "ADDITIVE"]:

482 precond.fct = haznics.precond_hx_curl_additive

483 elif prectype in ["multi", "MULTI", "Multi", "multiplicative",

484 "MULTIPLICATIVE"]:

485 precond.fct = haznics.precond_hx_curl_multiplicative

486 else: # default is additive

487 precond.fct = haznics.precond_hx_curl_additive

488 """

489 Precond.__init__(self, Acurl, "HXCurl_add", parameters, amgparam, precond)

490

491

492class HXDiv(Precond):

493 """

494 HX preconditioner from the HAZmath library for the div-div inner product

495 NB! Pdiv only works for 3D problems

496 """

497

498 def __init__(self, Adiv, V, parameters=None):

499 supports_mpi(False, 'HAZMath does not work in parallel')

500

501 # get auxiliary operators

502 mesh = V.mesh()

503 Pd = Pdiv(mesh)

504 Curl = discrete_curl(mesh)

505

506 # change data type for the matrices (to dCSRmat pointer)

507 Adiv_ptr = PETSc_to_dCSRmat(Adiv)

508 Pdiv_ptr = PETSc_to_dCSRmat(Pd)

509 Curl_ptr = PETSc_to_dCSRmat(Curl)

510

511 # initialize amg parameters (AMG_param pointer)

512 amgparam = haznics.AMG_param()

513

514 # set extra amg parameters

515 if parameters and isinstance(parameters, dict):

516 haznics.param_amg_set_dict(parameters, amgparam)

517 else:

518 parameters = {'dim': mesh.topology().dim(),

519 'prectype': haznics.PREC_HX_DIV_A}

520

521 # make sure coarse solver is always iterative

522 amgparam.coarse_solver = 0

523 # print (relevant) amg parameters

524 haznics.param_amg_print(amgparam)

525

526 # get dimension and type of HX precond application

527 try:

528 dim = parameters['dim']

529 except KeyError:

530 dim = 3

531

532 # add or multi

533 try:

534 prectype = parameters['prectype']

535 except KeyError:

536 prectype = haznics.PREC_HX_DIV_A

537

538 if dim == 3:

539 Pc = Pcurl(mesh)

540

541 # change data type for the Pcurl matrix (to dCSRmat pointer)

542 Pcurl_ptr = PETSc_to_dCSRmat(Pc)

543

544 # set HX DIV preconditioner (NB: this sets up both data and fct)

545 precond = haznics.create_precond_hxdiv_3D(Adiv_ptr, Pdiv_ptr,

546 Curl_ptr, Pcurl_ptr,

547 prectype, amgparam)

548

549 # if fail, setup returns null

550 if not precond:

551 raise RuntimeError(

552 "HXdiv data failed to set up (null pointer returned) ")

553 """

554 if prectype in ["add", "Add", "ADD", "additive", "ADDITIVE"]:

555 precond.fct = haznics.precond_hx_div_additive

556

557 elif prectype in ["multi", "MULTI", "Multi", "multiplicative",

558 "MULTIPLICATIVE"]:

559 precond.fct = haznics.precond_hx_div_multiplicative

560

561 else:

562 # default is additive

563 precond.fct = haznics.precond_hx_div_additive

564 """

565 Precond.__init__(self, Adiv, "HXDiv_add", parameters, amgparam, precond)

566

567 else:

568 # set HX DIV preconditioner (NB: this sets up both data and fct)

569 precond = haznics.create_precond_hxdiv_2D(Adiv_ptr, Pdiv_ptr,

570 Curl_ptr, prectype,

571 amgparam)

572

573 # if fail, setup returns null

574 if not precond:

575 raise RuntimeError(

576 "HXdiv data failed to set up (null pointer returned) ")

577 """

578 if prectype in ["add", "Add", "ADD", "additive", "ADDITIVE"]:

579 precond.fct = haznics.precond_hx_div_additive_2D

580

581 elif prectype in ["multi", "MULTI", "Multi", "multiplicative", "MULTIPLICATIVE"]:

582 precond.fct = haznics.precond_hx_div_multiplicative_2D

583

584 else:

585 # default is additive

586 precond.fct = haznics.precond_hx_div_additive_2D

587 """

588 Precond.__init__(self, Adiv, "HXDiv_add", parameters, amgparam, precond)

589

590

591# ----------------------------------- EOF ----------------------------------- #