# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import os
import numpy as np
import dolfin as df
import time
from custEM.misc import logger_print as lp
from petsc4py import PETSc
from dolfin import as_backend_type
import math
[docs]
class Solver:
"""
Solver class called internally from MOD instance.
Methods
-------
- solve_var_form_default()
solve variaitonal formulation using default FEniCS solver,
no support of symmetry, memory limits
- solve_system_default()
solve assembled system of equations with default FEniCS solver,
no support of symmetry, memory limits
- solve_system_mumps()
use MUMPS solver to solve system of equations,
no memory limits and symmetry support
- call_mumps()
call MUMPS solver within *solve_system_mumps()* method
- solve_var_form_iter()
use iteative solvers to solve variational formulation, in development
- solve_system_iter()
use iteative solvers to solve system of equations, in development
BoundaryConditions
------------------
- default choice of boundary conditions are implicit Neumann BC
- for zero Dirichlet BC, change *bc* argument in the 'solve' functions
or during assembly for all total field approaches
"""
def __init__(self, FS, FE, mumps_debug=False, serial_ordering=False):
"""
Initializes Solver instance for assembling and solution of FE systems.
Required arguments
------------------
- FS, type class
FunctionSpaces instance
- FE, type class
FE-approach (e.g., **E_vector**) instance
Keyword arguments
-----------------
- mumps_debug = False, type bool
set **True** during **MOD** class initialization for enabling
debugging MUMPS by adjusting the MUMPS internal log level
- serial_ordering = False, type bool
use serial ordering of *METIS* instead parallel ordering of
*SCOTCH*; leads to slight increase of solution time for *MUMPS*
but significantly reduces chance of potential *MUMPS* crashes
"""
self.FS = FS
self.FE = FE
#df.PETScOptions.set("mat_mumps_icntl_13", 1)
df.PETScOptions.set("mat_mumps_icntl_10", 0)
if serial_ordering:
df.PETScOptions.set("mat_mumps_icntl_28", 1)
else:
df.PETScOptions.set("mat_mumps_icntl_28", 2)
df.PETScOptions.set("mat_mumps_cntl_1", 1e-8)
# df.PETScOptions.set("mat_mumps_cntl_4", 0.)
# df.PETScOptions.set("mat_mumps_cntl_7", 1e-10)
# df.PETScOptions.set("mat_mumps_icntl_35", 2) # BLR test
# choose SCOTCH explicitly, not necessary as icntl 28 forces SCOTCH
# df.PETScOptions.set("mat_mumps_icntl_29", 1)
if mumps_debug:
df.PETScOptions.set("mat_mumps_icntl_4", 3)
[docs]
def solve_system_default(self, A, b):
"""
Use default petsc_lu solver to solve an assembled system of equations.
Required arguments
------------------
- A = left-hand-side matrix, type PETScMatrix
assembled LHS matrix
- b = right-hand-side vector, type PETScVector
assembled RHS vector
"""
t_0 = time.time()
lp(self.FE.MP.logger, 20,
' - default petsc LU solver initialized - ')
lp(self.FE.MP.logger, 20, '... solving system ...', pre_dash=False)
df.solve(A, self.FS.U.vector(), b.vector())
self.FS.solution_time = time.time() - t_0
lp(self.FE.MP.logger, 20,
' - solution time for main system of equations (s) --> ' +
str(self.FS.solution_time) + ' - ', pre_dash=False)
[docs]
def solve_system_mumps(self, A, b, fs=None, sym=False, out_of_core=False,
first_call=True, tx_selection=None):
"""
Use MUMPS solver to solve an assembled linear system of equations.
Required arguments
------------------
- A = left-hand-side matrix, type PETScMatrix
assembled LHS matrix
- b = right-hand-side vector, type PETScVector
assembled RHS vector
Keyword arguments
-----------------
- fs = None, type str or dolfin FunctionSpace
define FunctionSpace for solution function(s)
- sym = False, type int or bool
set to **1** for positive definite symmetric matrix
set to **2** or **True** for general structurally symmetric matrix
- out_of_core = False, type bool
set to **True** if MUMPS should be allowed to use disc space for
the solution process if the memory is completely exploited
- tx_selection = None, type list of int
specify if system should be solved for not all but only selected
right-hand sides; so far used to solve DC system only for grounded
Tx if there are also ungrounded ones in the same model
- first_call = True, type bool
specify if this is the first call of a solution based on the same
system matrix, re-use factorization if first_call=**False**
"""
print_dummy = False
if not first_call:
print_dummy = True
if fs is None:
self.U = [df.Function(self.FS.M) for ti in range(self.FE.n_tx)]
elif fs == 'scalar':
self.U = [df.Function(self.FS.S) for ti in range(self.FE.n_tx)]
else:
self.U = [df.Function(fs) for ti in range(self.FE.n_tx)]
counter = -1
second_call = True
for ti in range(self.FE.n_tx):
if not first_call:
if second_call:
lp(self.FE.MP.logger, 20,
'... solving additional right-hand sides ...',
pre_dash=False)
t2 = time.time()
second_call = False
try:
to_solve = b[ti].vector()
except AttributeError:
to_solve = b[ti]
if tx_selection is None:
if df.MPI.rank(df.MPI.comm_world) == 0:
self.U[ti].vector()[:]
self.call_mumps(A, to_solve, self.U[ti].vector(), sym,
out_of_core, first_call)
first_call = False
else:
if tx_selection[ti]:
self.call_mumps(A, to_solve, self.U[ti].vector(), sym,
out_of_core, first_call)
first_call = False
counter += 1
if counter > 0 and counter % 10 == 0:
lp(self.FE.MP.logger, 20,
' - solved right-hand '
'sides counter --> ' + str(counter) + ' - ')
if self.FE.n_tx > 1:
if print_dummy:
lp(self.FE.MP.logger, 20,
' - solution time for ' + str(self.FE.n_tx) +
' additional right-hand sides (s) --> ' +
str(int(time.time() - t2)) + ' - ', pre_dash=False)
else:
lp(self.FE.MP.logger, 20,
' - solution time for ' + str(self.FE.n_tx - 1) +
' additional right-hand sides (s) --> ' +
str(int(time.time() - t2)) + ' - ', pre_dash=False)
if fs is None:
self.FS.U = self.U
else:
return(self.U)
[docs]
def call_mumps(self, A, b, u, sym, out_of_core, first_call=True):
"""
Call MUMPS solver.
Required arguments
------------------
- A = left-hand-side matrix, type PETScMatrix
assembled LHS matrix
- b = right-hand-side vector, type PETScVector
assembled RHS vector
- u = solution vector, type PETScVector
empty solution vector to be filled
- see the definitions of the *solve_system_mumps()* method for the
description of the other input parameters
"""
t_0 = time.time()
if first_call:
if out_of_core:
df.PETScOptions.set("mat_mumps_icntl_22", 1)
lp(self.FE.MP.logger, 30,
'Warning! OUT_OF_CORE option applied which can makes the '
'solution process very slow if the available RAM is '
'exceeded. Continuing ...')
if type(sym) is bool and sym:
A.mat().setOption(PETSc.Mat.Option.SYMMETRIC, True)
lp(self.FE.MP.logger, 20,
' - general symmetric MUMPS solver initialized - ',
pre_dash=False)
elif sym == 1:
A.mat().setOption(PETSc.Mat.Option.SPD, True)
lp(self.FE.MP.logger, 20,
' - positive definite symmetric MUMPS solver '
'initialized - ', pre_dash=False)
elif sym == 2:
A.mat().setOption(PETSc.Mat.Option.SYMMETRIC, True)
lp(self.FE.MP.logger, 20,
' - general symmetric MUMPS solver initialized - ',
pre_dash=False)
else:
A.mat().setOption(PETSc.Mat.Option.SYMMETRIC, False)
A.mat().setOption(PETSc.Mat.Option.SPD, False)
lp(self.FE.MP.logger, 20,
' - unsymmetric MUMPS LU solver initialized - ',
pre_dash=False)
self.solver = df.PETScLUSolver(self.FS.mesh.mpi_comm(), A, "mumps")
lp(self.FE.MP.logger, 20,
'... solving system ...', pre_dash=False)
self.A = A
self.solver.solve(u, b)
self.FS.solution_time += time.time() - t_0
if first_call:
lp(self.FE.MP.logger, 20,
' - accumulated solution time for FE system (s) --> ' +
str(int(self.FS.solution_time)) + ' - ', pre_dash=False)
[docs]
def solve_system_iter_PRESB(self, Asys, Hh, b, B, A, inner, method,
flag_mt, fs=None, serial_ordering=False,
outer_tol=1e-8, inner_tol=1e-3,
maxiter=2000):
"""
Use generalised conjugate residual (GCR) method preconditioned by PRESB
to solve the system of equation (based on Weiss et al. (2023) Iterative solution
methods for 3D controlled-source electromagnetic modelling of geophysical
exploration scenarios)
The procedure is of "inner-outer" type, where the inner solver is a FGMRES method
preconditioned with the multi-grid based auxiliary-space preconditioner AMS.
The implementation is largely based on available functions of the PETSc and hypre
libraries.
The convergence history of the outer solver (i.e. the residual for each iteration)
is written to outerrelres.txt.
Implemented and commented by M. Weiss (January 2024)
Required arguments (assembly of matrices in .py)
------------------
- Asys = 2 by 2 system block matrix, type PETScMatrix
- H = matrix to be solved by inner solver, type PETScMatrix
- b = right-hand-side vector, type PETScVector
- A = diagonal block of system matrix, type PETScMatrix
- B = off-diagonal block of system matrix, type PETScMatrix
- inner = 'iterative' or 'directly' depending on whether the inner systems
are to be solved by an iterative or direct method
(first argument in bundled input variable solvingapproach required to
initialise a custEM model)
- method = 'fgmres' by default, alternatively 'cg' or 'gcr'
(set when calling solve_main_problem)
- flag_mt = False (default) or True (automatically set if MT problem is solved)
- fs = None, type str or dolfin FunctionSpace
define FunctionSpace for solution function(s)
- serial_ordering = False (default) or True
to use serial ordering of *METIS* instead of parallel ordering by
*SCOTCH*; leads to slight increase of solution time for *MUMPS*
but significantly reduces chance of potential *MUMPS* crashes
Keyword arguments
-----------------
- outer_tol = 1e-8, type float
tolerance of outer GCR solver
- inner_tol = 1e-3, type float
tolerance of inner solver
- maxiter = 2000, type int
maximum number of iterations
"""
# Matrices
Asystem=as_backend_type(Asys).mat()
A11=as_backend_type(A).mat()
A12 = as_backend_type(B).mat()
H = as_backend_type(Hh).mat()
# Vectors
r=PETSc.Vec()
f1, f2 = PETSc.Vec(), PETSc.Vec()
g, h = PETSc.Vec(), PETSc.Vec()
Ag = PETSc.Vec()
Aw = PETSc.Vec()
w=PETSc.Vec()
# Inner iterative solver and preconditioner
ksp=PETSc.KSP()
ksp.create(PETSc.COMM_WORLD)
ksp.setOperators(H,H)
if(inner=='iterative'):
ksp.setType(method)
ksp.setTolerances(rtol=inner_tol,max_it=200)
res_inner=ksp.setConvergenceHistory(200,True)
pcc=ksp.getPC()
pcc.setType('hypre')
pcc.setHYPREType('ams')
V=df.FunctionSpace(self.FS.mesh,'N1curl',1)
P1=df.FunctionSpace(self.FS.mesh,'Lagrange',1)
G=df.DiscreteOperators.build_gradient(V,P1)
pcc.setHYPREDiscreteGradient(as_backend_type(G).mat())
constants=[df.Function(V) for i in range(3)]
for i, c in enumerate(constants):
direction = [1.0 if i==j else 0.0 for j in range(3)]
c.interpolate(df.Constant(direction))
cvecs = [as_backend_type(constant.vector()).vec() for constant in constants]
pcc.setHYPRESetEdgeConstantVectors(cvecs[0],cvecs[1],cvecs[2])
elif(inner=='directly'):
ksp.setType('preonly')
pcc=ksp.getPC()
pcc.setType('lu')
pcc.setFactorSolverType('mumps')
pcc.setFactorSetUpSolverType()
pcc.getFactorMatrix()
opts = PETSc.Options()
if serial_ordering:
opts["mat_mumps_icntl_28"] = 1
else:
opts["mat_mumps_icntl_28"] = 2
opts["mat_mumps_icntl_10"] = 0
opts["mat_mumps_cntl_1"] = 1e-8
else:
lp(self.FE.MP.logger, 50,
'Error! Unknown flag set for inner. '
'Aborting ...')
raise SystemExit
# Additional preparation
H.getColumnVector(1,f1)
H.getColumnVector(1,f2)
H.getColumnVector(1,g)
H.getColumnVector(1,h)
H.getColumnVector(1,Ag)
f1.zeroEntries()
f2.zeroEntries()
g.zeroEntries()
h.zeroEntries()
Ag.zeroEntries()
if fs is None:
self.U = [df.Function(self.FS.M) for ti in range(self.FE.n_tx)]
elif fs == 'scalar':
self.U = [df.Function(self.FS.S) for ti in range(self.FE.n_tx)]
else:
self.U = [df.Function(fs) for ti in range(self.FE.n_tx)]
for ti in range(self.FE.n_tx):
try:
to_solve = b[ti].vector()
except AttributeError:
to_solve = b[ti]
#if tx_selection is None:
if df.MPI.rank(df.MPI.comm_world) == 0:
self.U[ti].vector()[:]
t_0 = time.time()
restart=100
# Initialise required arrays, vectors, matrices and other variables
y=np.zeros(restart)
ADd=np.zeros(restart)
resvec=np.zeros(maxiter+1)
AD=[]
DD=[]
u=as_backend_type(self.U[ti].vector()).vec()
u.zeroEntries()
u.assemble()
# Scalars
itout=0
it2=0
its=0
relres=1.0
resold=2.0
tol=outer_tol
# writing out matrix or vector using PetscViewer
#Mat = as_backend_type(mat).mat()
#myviewer=PETSc.Viewer().createBinary('f.dat','w')
#myviewer(f)
# Set initial residual to RHS
f = as_backend_type(to_solve).vec()
r=f.copy()
w=r.copy()
Aw=w.copy()
# Copmute initial norm
nf=r.norm()
resvec[0]=nf
# Implementation of the PRESB-preconditioned GCR solver
while (relres >= tol) and (itout < maxiter) and ((resold-relres) > tol):
# START OF INNER SOLVER
# Extract RHS from residual vector
istart,iend=f1.getOwnershipRange()
j=f1.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=r.getValue(ig*2+1)
f1.setValue(ig,af1,addv=False)
istart,iend=f2.getOwnershipRange()
j=f2.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=r.getValue(ig*2)
f2.setValue(ig,af1,addv=False)
f1.assemble()
f2.assemble()
f2.aypx(1,f1)
ksp.setUp()
ksp.solve(f2,g)
its=ksp.getIterationNumber()
res_inner=ksp.getConvergenceHistory()
#if df.MPI.comm_world.rank == 0 and inner=='iterative':
# with open('innerrelres_s1.txt','ab') as m:
# np.savetxt(m,[itout],fmt='%d')
# np.savetxt(m, [res_inner[0:its]], fmt='%1.9e')
A11.mult(g,Ag)
f1.axpy(-1,Ag)
ksp.setUp()
ksp.solve(f1,h)
its=ksp.getIterationNumber()
res_inner=ksp.getConvergenceHistory()
#if df.MPI.comm_world.rank == 0 and inner=='iterative':
# with open('innerrelres_s2.txt','ab') as n:
# np.savetxt(n,[itout],fmt='%d')
# np.savetxt(n, [res_inner[0:its]], fmt='%1.9e')
# Build solution w of inner solves
g.aypx(1,h)
h.scale(-1)
# END OF INNER SOLVER
istart,iend=g.getOwnershipRange()
j=g.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=g.getValue(ig)
w.setValue(ig*2+1,af1,addv=False)
istart,iend=h.getOwnershipRange()
j=h.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=h.getValue(ig)
w.setValue(ig*2,af1,addv=False)
w.assemble()
Asystem.mult(w,Aw)
itout=itout+1
it2=it2+1
if (it2>1):
t=Aw.copy()
t.zeroEntries()
wt=w.copy()
wt.zeroEntries()
At=Aw.copy()
At.zeroEntries()
for i in range(0,it2-1):
At=AD[i].copy()
y[i]=At.dot(Aw)
y[i]=y[i]/ADd[i]
y[i]=y[i]*(-1)
t.axpy(y[i],AD[i])
wt.axpy(y[i],DD[i])
Aw.axpy(1,t)
At=Aw.copy()
AD.append(At)
nw2=Aw.norm()
nw2=nw2*nw2
ADd[it2-1]=nw2
w.axpy(1,wt)
wt=w.copy()
DD.append(wt)
else:
At=Aw.copy()
AD.append(At)
nw2=Aw.norm()
nw2=nw2*nw2
ADd[it2-1]=nw2
wt=w.copy()
DD.append(wt)
al=r.dot(Aw)
al=al/nw2
u.axpy(al,w)
al=al*(-1.0)
r.axpy(al,Aw)
nr=r.norm()
resold=relres
relres=nr/nf
resvec[itout]=nr
#restart condition
if it2 >= restart:
it2=0
DD=[]
AD=[]
y=np.zeros(restart)
ADd=np.zeros(restart)
ut=u.copy();
u.zeroEntries()
istart,iend=ut.getOwnershipRange()
j=ut.getLocalSize()
af1=0
if flag_mt==False:
for i in range(0,j):
ig=i+istart
af1=0
af1=ut.getValue(ig)
if (ig%2==0):
af1=af1*(-1)
u.setValue(ig+1,af1,addv=False)
else:
u.setValue(ig-1,af1,addv=False)
u.assemble()
u.scale(-1)
elif flag_mt==True:
for i in range(0,j):
ig=i+istart
af1=0
af1=ut.getValue(ig)
if (ig%2==0):
#af1=af1*(-1)
u.setValue(ig,af1,addv=False)
else:
#af1=af1*(-1)
u.setValue(ig,af1,addv=False)
u.assemble()
else:
lp(self.logger, 50, "Unknown flag. Abort simulation!")
raise SystemExit
if df.MPI.comm_world.rank == 0:
with open('outerrelres.txt','ab') as f:
np.savetxt(f, [resvec[0:itout]], fmt='%1.9e')
self.U[ti].vector()[:]=df.PETScVector(u)
#self.FS.U = self.U
if fs is None:
self.FS.U = self.U
else:
return(self.U)
stop = time.time() - t_0
lp(self.FE.MP.logger, 20,
' - solution time for main system of equations --> ' +
str(stop) + ' - ',pre_dash=False, post_dash=False)
[docs]
def solve_system_iter_BD(self, Asys, Hh, b, B, A, inner, method, flag_mt,
fs=None, serial_ordering=False,outer_tol=1e-8,
inner_tol=1e-3, maxiter=2000):
"""
Use generalised conjugate residual (GCR) method preconditioned by block-diagonal
preconditioner to solve the system of equation (based on Grayver & Buerg (2014)
Robust and scalable 3-D geo-electromagnetic modelling approach
using the finite element method and Grayver & Kolev (2015)
Large-scale 3D geoelectromagnetic modeling using parallel
adaptive high-order finite element method)
The procedure is of "inner-outer" type, where the inner solver is a FGMRES method
preconditioned with the multi-grid based auxiliary-space preconditioner AMS.
The implementation is largely based on available functions of the PETSc and hypre
libraries.
The convergence history of the outer solver (i.e. the residual for each iteration)
is written to outerrelres.txt.
Implemented and commented by M. Weiss (January 2024)
Required arguments (assembly of matrices in .py)
------------------
- Asys = 2 by 2 system block matrix, type PETScMatrix
- H = matrix to be solved by inner solver, type PETScMatrix
- b = right-hand-side vector, type PETScVector
- A = diagonal block of system matrix, type PETScMatrix
- B = off-diagonal block of system matrix, type PETScMatrix
- inner = 'iterative' or 'directly' depending on whether the inner systems
are to be solved by an iterative or direct method
(first argument in bundled input variable solvingapproach required to
initialise a custEM model)
- method = 'fgmres' by default, alternatively 'cg' or 'gcr'
(set when calling solve_main_problem)
- flag_mt = False (default) or True (automatically set if MT problem is solved)
- fs = None, type str or dolfin FunctionSpace
define FunctionSpace for solution function(s)
- serial_ordering = False (default) or True
to use serial ordering of *METIS* instead of parallel ordering by
*SCOTCH*; leads to slight increase of solution time for *MUMPS*
but significantly reduces chance of potential *MUMPS* crashes
Keyword arguments
-----------------
- outer_tol = 1e-8, type float
tolerance of outer GCR solver
- inner_tol = 1e-3, type float
tolerance of inner solver
- maxiter = 2000, type int
maximum number of iterations
"""
# Matrices
Asystem=as_backend_type(Asys).mat()
A11=as_backend_type(A).mat()
A12 = as_backend_type(B).mat()
H = as_backend_type(Hh).mat()
# Vectors
r=PETSc.Vec()
f1, f2 = PETSc.Vec(), PETSc.Vec()
g, h = PETSc.Vec(), PETSc.Vec()
Aw = PETSc.Vec()
w=PETSc.Vec()
# Inner iterative solver and preconditioner
ksp=PETSc.KSP()
ksp.create(PETSc.COMM_WORLD)
ksp.setOperators(H,H)
if(inner=='iterative'):
ksp.setType(method)
ksp.setTolerances(rtol=inner_tol,max_it=200)
res_inner=ksp.setConvergenceHistory(200,True)
pcc=ksp.getPC()
pcc.setType('hypre')
pcc.setHYPREType('ams')
V=df.FunctionSpace(self.FS.mesh,'N1curl',1)
P1=df.FunctionSpace(self.FS.mesh,'Lagrange',1)
G=df.DiscreteOperators.build_gradient(V,P1)
pcc.setHYPREDiscreteGradient(as_backend_type(G).mat())
constants=[df.Function(V) for i in range(3)]
for i, c in enumerate(constants):
direction = [1.0 if i==j else 0.0 for j in range(3)]
c.interpolate(df.Constant(direction))
cvecs = [as_backend_type(constant.vector()).vec() for constant in constants]
pcc.setHYPRESetEdgeConstantVectors(cvecs[0],cvecs[1],cvecs[2])
elif(inner=='directly'):
ksp.setType('preonly')
pcc=ksp.getPC()
pcc.setType('lu')
pcc.setFactorSolverType('mumps')
pcc.setFactorSetUpSolverType()
pcc.getFactorMatrix()
opts = PETSc.Options()
if serial_ordering:
opts["mat_mumps_icntl_28"] = 1
else:
opts["mat_mumps_icntl_28"] = 2
opts["mat_mumps_icntl_10"] = 0
opts["mat_mumps_cntl_1"] = 1e-8
else:
lp(self.FE.MP.logger, 50,
'Error! Unknown flag set for inner. '
'Aborting ...')
raise SystemExit
# Additional preparation
H.getColumnVector(1,f1)
H.getColumnVector(1,f2)
H.getColumnVector(1,g)
H.getColumnVector(1,h)
f1.zeroEntries()
f2.zeroEntries()
g.zeroEntries()
h.zeroEntries()
if fs is None:
self.U = [df.Function(self.FS.M) for ti in range(self.FE.n_tx)]
elif fs == 'scalar':
self.U = [df.Function(self.FS.S) for ti in range(self.FE.n_tx)]
else:
self.U = [df.Function(fs) for ti in range(self.FE.n_tx)]
for ti in range(self.FE.n_tx):
try:
to_solve = b[ti].vector()
except AttributeError:
to_solve = b[ti]
#if tx_selection is None:
if df.MPI.rank(df.MPI.comm_world) == 0:
self.U[ti].vector()[:]
# Implementation of the PRESB-preconditioned GCR solver
t_0 = time.time()
restart=100
# Initialise required arrays, vectors, matrices and other variables
y=np.zeros(restart)
ADd=np.zeros(restart)
resvec=np.zeros(maxiter+1)
AD=[]
DD=[]
u=as_backend_type(self.U[ti].vector()).vec()
u.zeroEntries()
# Scalars
itout=0
it2=0
its=0
relres=1.0
resold=2.0
tol=outer_tol
# Set initial residual to RHS
f = as_backend_type(to_solve).vec()
r=f.copy()
w=r.copy()
Aw=w.copy()
# Copmute initial norm
nf=r.norm()
resvec[0]=nf
while (relres >= tol) and (itout < maxiter) and ((resold-relres) > tol):
# START OF INNER SOLVER
# Extract RHS from residual vector
istart,iend=f1.getOwnershipRange()
j=f1.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=r.getValue(ig*2+1)
f1.setValue(ig,af1,addv=False)
istart,iend=f2.getOwnershipRange()
j=f2.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=r.getValue(ig*2)
f2.setValue(ig,af1,addv=False)
f1.assemble()
f2.assemble()
ksp.setUp()
ksp.solve(f2,g)
its=ksp.getIterationNumber()
res_inner=ksp.getConvergenceHistory()
#if df.MPI.comm_world.rank == 0 and inner=='iterative':
# with open('innerrelres_s1.txt','ab') as m:
# np.savetxt(m,[itout],fmt='%d')
# np.savetxt(m, [res_inner[0:its]], fmt='%1.9e')
ksp.setUp()
ksp.solve(f1,h)
its=ksp.getIterationNumber()
res_inner=ksp.getConvergenceHistory()
#if df.MPI.comm_world.rank == 0 and inner=='iterative':
# with open('innerrelres_s2.txt','ab') as n:
# np.savetxt(n,[itout],fmt='%d')
# np.savetxt(n, [res_inner[0:its]], fmt='%1.9e')
# END OF INNER SOLVER
istart,iend=g.getOwnershipRange()
j=g.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=g.getValue(ig)
w.setValue(ig*2+1,af1,addv=False)
istart,iend=h.getOwnershipRange()
j=h.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=h.getValue(ig)
w.setValue(ig*2,af1,addv=False)
w.assemble()
Asystem.mult(w,Aw)
itout=itout+1
it2=it2+1
if (it2>1):
t=Aw.copy()
t.zeroEntries()
wt=w.copy()
wt.zeroEntries()
At=Aw.copy()
At.zeroEntries()
for i in range(0,it2-1):
At=AD[i].copy()
y[i]=At.dot(Aw)
y[i]=y[i]/ADd[i]
y[i]=y[i]*(-1)
t.axpy(y[i],AD[i])
wt.axpy(y[i],DD[i])
Aw.axpy(1,t)
At=Aw.copy()
AD.append(At)
nw2=Aw.norm()
nw2=nw2*nw2
ADd[it2-1]=nw2
w.axpy(1,wt)
wt=w.copy()
DD.append(wt)
else:
At=Aw.copy()
AD.append(At)
nw2=Aw.norm()
nw2=nw2*nw2
ADd[it2-1]=nw2
wt=w.copy()
DD.append(wt)
al=r.dot(Aw)
al=al/nw2
u.axpy(al,w)
al=al*(-1.0)
r.axpy(al,Aw)
nr=r.norm()
resold=relres
relres=nr/nf
resvec[itout]=nr
#restart condition
if it2 >= restart:
it2=0
DD=[]
AD=[]
y=np.zeros(restart)
ADd=np.zeros(restart)
ut=u.copy();
u.zeroEntries()
istart,iend=ut.getOwnershipRange()
j=ut.getLocalSize()
af1=0
if flag_mt==False:
for i in range(0,j):
ig=i+istart
af1=0
af1=ut.getValue(ig)
if (ig%2==0):
af1=af1*(-1)
u.setValue(ig+1,af1,addv=False)
else:
u.setValue(ig-1,af1,addv=False)
u.assemble()
u.scale(-1)
elif flag_mt==True:
for i in range(0,j):
ig=i+istart
af1=0
af1=ut.getValue(ig)
if (ig%2==0):
#af1=af1*(-1)
u.setValue(ig,af1,addv=False)
else:
#af1=af1*(-1)
u.setValue(ig,af1,addv=False)
u.assemble()
else:
lp(self.logger, 50, "Unknown flag. Abort simulation!")
raise SystemExit
if df.MPI.comm_world.rank == 0:
with open('outerrelres.txt','ab') as f:
np.savetxt(f, [resvec[0:itout]], fmt='%1.9e')
self.U[ti].vector()[:]=df.PETScVector(u)
#self.FS.U = self.U
if fs is None:
self.FS.U = self.U
else:
return(self.U)
stop = time.time() - t_0
lp(self.FE.MP.logger, 20,
' - solution time for main system of equations --> ' +
str(stop) + ' - ',pre_dash=False, post_dash=False)
[docs]
def solve_system_iter_Hfield(self, H, RHS, fs=None, method='fgmres'):
"""
Computes the magnetic fields from the electric fields using the FGMRES
method preconditioned with the auxiliary-space solver AMS.
Presents an alternative to deriving the magnetic fields by the direct
solver MUMPS.
Can be employed by setting the 4th argument of the bundled input variable
solvingapproach (when initialising a custem model) to 'H_iterative'.
The convergence history of the solver (i.e. the residual for each iteration)
is written to Hfielditer.txt.
Implemented and commented by M. Weiss (January 2024)
"""
# Matrices
Asystem=as_backend_type(H).mat()
# Set up solver
ksp=PETSc.KSP()
ksp.create(PETSc.COMM_WORLD)
ksp.setOperators(Asystem,Asystem)
ksp.setType(method)
ksp.setTolerances(rtol=1e-8,max_it=2000)
res_inner=ksp.setConvergenceHistory(2000,True)
pcc=ksp.getPC()
pcc.setType('hypre')
pcc.setHYPREType('ams')
V=df.FunctionSpace(self.FS.mesh,'N1curl',1)
P1=df.FunctionSpace(self.FS.mesh,'Lagrange',1)
G=df.DiscreteOperators.build_gradient(V,P1)
pcc.setHYPREDiscreteGradient(as_backend_type(G).mat())
constants=[df.Function(V) for i in range(3)]
for i, c in enumerate(constants):
direction = [1.0 if i==j else 0.0 for j in range(3)]
c.interpolate(df.Constant(direction))
cvecs = [as_backend_type(constant.vector()).vec() for constant in constants]
pcc.setHYPRESetEdgeConstantVectors(cvecs[0],cvecs[1],cvecs[2])
if fs is None:
self.UV = [df.Function(self.FS.M1) for ti in range(self.FE.n_tx)]
elif fs == 'scalar':
self.UV = [df.Function(self.FS.S) for ti in range(self.FE.n_tx)]
else:
self.UV = [df.Function(fs) for ti in range(self.FE.n_tx)]
for ti in range(self.FE.n_tx):
try:
to_solve = RHS[ti].vector()
except AttributeError:
to_solve = RHS[ti]
t_0 = time.time()
resvec=np.zeros(2000)
# Vectors
u=as_backend_type(self.UV[ti].vector()).vec()
u.zeroEntries()
r=as_backend_type(to_solve).vec()
ksp.setUp()
ksp.solve(r,u)
its=ksp.getIterationNumber()
res_inner=ksp.getConvergenceHistory()
if df.MPI.comm_world.rank == 0:
with open('Hfielditer.txt','ab') as m:
np.savetxt(m, [res_inner[0:its]], fmt='%1.9e')
self.UV[ti].vector()[:]=df.PETScVector(u)
#self.FS.U = self.U
if fs is None:
self.FS.V = self.UV
else:
return(self.UV)
stop = time.time() - t_0
lp(self.FE.MP.logger, 20,
' - solution time for obtaining H-fields --> ' +
str(stop) + ' - ',pre_dash=False, post_dash=False)
[docs]
def build_PRESB_Jacobian(self, Asys, Hh, B, A, inner_tol=1e-1):
"""
Sets up the iterative solver and system to be solved.
(GCR preconditioned with PRESB adjusted for the computation of
entries of the Jacobian)
inner_tol = 1e-1 (default): set tolerance for inner solver
(in inv_base.py)
Function returns required variables (matrices, vectors and inner solver)
for GCR method
Implemented and commented by M. Weiss (January 2024)
"""
# Matrices
Asystem=as_backend_type(Asys).mat()
A11=as_backend_type(A).mat()
A12 = as_backend_type(B).mat()
H = as_backend_type(Hh).mat()
# Vectors
r=PETSc.Vec()
f1, f2 = PETSc.Vec(), PETSc.Vec()
g, h = PETSc.Vec(), PETSc.Vec()
Ag = PETSc.Vec()
Aw = PETSc.Vec()
w=PETSc.Vec()
AD=[]
DD=[]
# Inner iterative solver and preconditioner
ksp=PETSc.KSP()
ksp.create(PETSc.COMM_WORLD)
ksp.setOperators(H,H)
ksp.setType('fgmres')
ksp.setTolerances(rtol=inner_tol,max_it=200)
res_inner=ksp.setConvergenceHistory(200,True)
pcc=ksp.getPC()
pcc.setType('hypre')
pcc.setHYPREType('ams')
V=df.FunctionSpace(self.FS.mesh,'N1curl',1)
P1=df.FunctionSpace(self.FS.mesh,'Lagrange',1)
G=df.DiscreteOperators.build_gradient(V,P1)
pcc.setHYPREDiscreteGradient(as_backend_type(G).mat())
constants=[df.Function(V) for i in range(3)]
for i, c in enumerate(constants):
direction = [1.0 if i==j else 0.0 for j in range(3)]
c.interpolate(df.Constant(direction))
cvecs = [as_backend_type(constant.vector()).vec() for constant in constants]
pcc.setHYPRESetEdgeConstantVectors(cvecs[0],cvecs[1],cvecs[2])
# Additional preparation
H.getColumnVector(1,f1)
H.getColumnVector(1,f2)
H.getColumnVector(1,g)
H.getColumnVector(1,h)
H.getColumnVector(1,Ag)
f1.zeroEntries()
f2.zeroEntries()
g.zeroEntries()
h.zeroEntries()
Ag.zeroEntries()
return ksp,f1,f2,g,h,Ag,H,A11,Asystem,Aw,w,r,DD,AD
[docs]
def solve_PRESB_Jacobian(self,ksp,f1,f2,g,h,Ag,H,A11,Asystem,Aw,w,r,DD,AD,
b,fs=None, flag_mt=False, outer_tol=1e-4,
sol=None):
"""
Solves the system for multiple RHS using GCR preconditioned with PRESB
to obtain the entries of the Jacobian (sensitivities)
outer_tol = 1e-4 (default): set tolerance for outer solver
(in inv_base.py)
Function returns sensitivites for one receiver
(equivalent to a column of the Jacobian)
Function writes out convergence history (i.e. the residual for each iteration)
to outerrelresPID.txt, where PID is the process ID
Implemented and commented by M. Weiss (January 2024)
"""
# Scalars
itout=0
it2=0
its=0
relres=1.0
resold=2.0
tol=outer_tol
AD=[]
DD=[]
restart=100
y=np.zeros(restart)
ADd=np.zeros(restart)
self.U = df.Function(fs)
#if df.MPI.rank(df.MPI.comm_world) == 0:
# self.U.vector()[:]
#u=as_backend_type(self.U.vector()).vec()
#u.zeroEntries()
if sol==None:
u=as_backend_type(self.U.vector()).vec()
u.zeroEntries()
elif sol==0:
u=as_backend_type(self.U.vector()).vec()
u.zeroEntries()
else:
u=as_backend_type(sol.vector()).vec()
to_solve = b
resvec=np.zeros(2000+1)
# Set initial residual to RHS
f = as_backend_type(to_solve).vec()
r=f.copy()
w=r.copy()
Aw=w.copy()
# Copmute initial norm
nf=r.norm()
resvec[0]=nf
# Implementation of the PRESB-preconditioned GCR solver
while (relres >= tol) and (itout < 2000) and ((resold-relres) > tol):
# START OF INNER SOLVER
# Extract RHS from residual vector
istart,iend=f1.getOwnershipRange()
j=f1.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=r.getValue(ig*2+1)
f1.setValue(ig,af1,addv=False)
istart,iend=f2.getOwnershipRange()
j=f2.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=r.getValue(ig*2)
f2.setValue(ig,af1,addv=False)
f1.assemble()
f2.assemble()
f2.aypx(1,f1)
#g.zeroEntries()
#g.assemble()
ksp.setUp()
ksp.solve(f2,g)
#its=ksp.getIterationNumber()
#res_inner=ksp.getConvergenceHistory()
#if df.MPI.comm_world.rank == 0 and inner=='iterative':
# with open('innerrelres_s1.txt','ab') as m:
# np.savetxt(m,[itout],fmt='%d')
# np.savetxt(m, [res_inner[0:its]], fmt='%1.9e')
#Ag.zeroEntries()
#Ag.assemble()
A11.mult(g,Ag)
f1.axpy(-1,Ag)
#h.zeroEntries()
#h.assemble()
ksp.setUp()
ksp.solve(f1,h)
#its=ksp.getIterationNumber()
#res_inner=ksp.getConvergenceHistory()
#if df.MPI.comm_world.rank == 0 and inner=='iterative':
# with open('innerrelres_s2.txt','ab') as n:
# np.savetxt(n,[itout],fmt='%d')
# np.savetxt(n, [res_inner[0:its]], fmt='%1.9e')
# Build solution w of inner solves
g.aypx(1,h)
h.scale(-1)
#w=r.copy()
#w.zeroEntries()
#w.assemble()
istart,iend=g.getOwnershipRange()
j=g.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=g.getValue(ig)
w.setValue(ig*2+1,af1,addv=False)
istart,iend=h.getOwnershipRange()
j=h.getLocalSize()
af1=0
for i in range(0,j):
ig=i+istart
af1=h.getValue(ig)
w.setValue(ig*2,af1,addv=False)
w.assemble()
#Aw=w.copy()
#Aw.zeroEntries()
#Aw.assemble()
Asystem.mult(w,Aw)
itout=itout+1
it2=it2+1
if (it2>1):
t=Aw.copy()
t.zeroEntries()
wt=w.copy()
wt.zeroEntries()
At=Aw.copy()
At.zeroEntries()
for i in range(0,it2-1):
At=AD[i].copy()
y[i]=At.dot(Aw)
y[i]=y[i]/ADd[i]
y[i]=y[i]*(-1)
t.axpy(y[i],AD[i])
wt.axpy(y[i],DD[i])
Aw.axpy(1,t)
At=Aw.copy()
AD.append(At)
nw2=Aw.norm()
nw2=nw2*nw2
ADd[it2-1]=nw2
w.axpy(1,wt)
wt=w.copy()
DD.append(wt)
else:
At=Aw.copy()
AD.append(At)
nw2=Aw.norm()
nw2=nw2*nw2
ADd[it2-1]=nw2
wt=w.copy()
DD.append(wt)
al=r.dot(Aw)
al=al/nw2
u.axpy(al,w)
al=al*(-1.0)
r.axpy(al,Aw)
nr=r.norm()
resold=relres
relres=nr/nf
resvec[itout]=nr
#restart condition
if it2 >= restart:
it2=0
DD=[]
AD=[]
y=np.zeros(restart)
ADd=np.zeros(restart)
ut=u.copy();
u.zeroEntries()
istart,iend=ut.getOwnershipRange()
j=ut.getLocalSize()
af1=0
if flag_mt==False:
for i in range(0,j):
ig=i+istart
af1=0
af1=ut.getValue(ig)
if (ig%2==0):
af1=af1*(-1)
u.setValue(ig+1,af1,addv=False)
else:
u.setValue(ig-1,af1,addv=False)
u.assemble()
u.scale(-1)
elif flag_mt==True:
for i in range(0,j):
ig=i+istart
af1=0
af1=ut.getValue(ig)
if (ig%2==0):
#af1=af1*(-1)
u.setValue(ig,af1,addv=False)
else:
#af1=af1*(-1)
u.setValue(ig,af1,addv=False)
u.assemble()
else:
lp(self.logger, 50, "Unknown flag. Abort simulation!")
raise SystemExit
if df.MPI.comm_world.rank == 0:
PID = os.getpid()
with open('outerrelres' + str(PID) +'.txt','ab') as fil:
np.savetxt(fil, [resvec[0:itout]], fmt='%1.9e')
u.scale(-1)
self.U.vector()[:]=df.PETScVector(u)
#self.FS.U = self.U
if fs is None:
self.FS.U = self.U
else:
return(self.U)