# -*- coding: utf-8 -*-
#
# This file is subject to the terms and conditions defined in
# file 'LICENSE.txt', which is part of this source code package.
#
from typing import Dict, Any, Tuple
import numpy as np
from scipy.sparse import linalg as splinalg
from scipy import sparse
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from Muscat.Types import MuscatFloat, MuscatIndex, ArrayLike
from Muscat.Helpers.Logger import Debug, Info
import Muscat.LinAlg.NativeEigenSolver as NativeEigenSolver
methodFactory: Dict[str, Any] = {}
defaultMethod = "Projection"
[docs]
class ConstraintsHolder:
"""
Constrained linear problem: Class to store and to enforce constraints for a
quadratic system of the form:
minimize 1/2 u.K.u - u.f
subject to:
A.u = b
This class can store, and manipulate the system Au=b and enforce the constraints
by different methods (penalization, substitution, lagrange multiplier)
the first set of function is used to fill this class with the system Au=b:
NextEquation must be called after the call of a set of : AddFactor, AddConstant
One per constraint
AddEquationSparse, AddEquation, AddEquationsFromIJV to add a full
constraint at once
Compact will clean (zeros, and duplicated) entries
"""
def __init__(self):
super().__init__()
self.constraints = []
self.methodName = "Projection"
self.Reset()
self.SetConstraintsMethod(defaultMethod)
self.SetNumberOfDofs(0)
self.tol = 1e-8 # tolerance for the detection of redundant constraints
# ----------------------- for population of the class --------------------------
[docs]
def Reset(self):
self.numberOfEquations = 0
self.rows = []
self.cols = []
self.vals = []
self.op = None
self.rhs = []
self.SetConstraintsMethod(self.methodName)
[docs]
def AddConstraint(self, constraints):
self.constraints.append(constraints)
[docs]
def ComputeConstraintsEquations(self, mesh, unknownFields, solutionFields=None, reactionFields=None):
self.Reset()
for cons in self.constraints:
cons.GenerateEquations(mesh, unknownFields, self, solutionFields=solutionFields, reactionFields=reactionFields)
Info(str(self.nbdof))
Info(len(self.constraints))
Info(self.numberOfEquations)
[docs]
def SetNumberOfDofs(self, nbdof):
self.nbdof = nbdof
[docs]
def AddFactor(self, ddl, factor):
if factor == 0:
return
self.rows.append(self.numberOfEquations)
self.cols.append(ddl)
self.vals.append(factor)
[docs]
def AddConstant(self, constant):
if constant == 0:
return
self.rows.append(self.numberOfEquations)
self.cols.append(self.nbdof)
self.vals.append(constant)
[docs]
def AddEquation(self, vals, const=0):
for j, val in enumerate(vals):
self.AddFactor(j, val)
if const:
self.AddConstant(const)
self.NextEquation()
[docs]
def AddEquationSparse(self, index, vals, const):
for j, v in zip(index, vals):
self.AddFactor(j, v)
if const:
self.AddConstant(const)
self.NextEquation()
[docs]
def AddEquationsFromIJV(self, ei, ej, ev):
ei = np.require(ei, dtype=MuscatIndex)
ej = np.require(ej, dtype=MuscatIndex)
ev = np.require(ev, dtype=MuscatFloat)
s = np.unique(ei)
for i in s:
mask = np.equal(ei, i)
self.AddEquationSparse(ej[mask], ev[mask], 0)
[docs]
def NextEquation(self):
self.numberOfEquations += 1
[docs]
def Compact(self):
r = self.numberOfEquations
c = self.nbdof + 1
res = coo_matrix((self.vals, (self.rows, self.cols)), shape=((r, c)), copy=False, dtype=MuscatFloat)
self.rows = res.row
self.cols = res.col
self.vals = res.data
[docs]
def AddEquationsFromOperatorAAndb(self, A, b=None):
self.rows.extend(A.row + self.numberOfEquations)
self.cols.extend(A.col)
self.vals.extend(A.data)
if b is not None:
self.rows.extend(np.arange(A.shape[0], dtype=MuscatIndex) + self.numberOfEquations)
self.cols.extend(np.full(A.shape[0], self.nbdof, dtype=MuscatIndex))
self.vals.extend(b)
self.numberOfEquations += A.shape[0]
[docs]
def SetOp(self, op, rhs=None):
self.op = sparse.csr_matrix(op)
self.rhs = rhs
self.nbdof = self.op.shape[1]
[docs]
def Residual(self, u: ArrayLike) -> Tuple[np.ndarray, np.ndarray]:
"""Recover A.u and b
Parameters
----------
u : ArrayLike
Vector of unknowns
Returns
-------
Tuple[np.ndarray,np.ndarray]
se result of
"""
self.Compact()
Ab = self.ToSparseFull()
utemp = np.zeros(len(u) + 1, dtype=MuscatFloat)
utemp[0:-1] = u
res = Ab.dot(utemp)
utemp[0:-1] = 0
utemp[-1] = 1
return res, Ab.dot(utemp)
# ----------------------- To recover the system of equations -------------------
[docs]
def ToSparse(self):
"""
Function to convert to the constraints system Ax=B to a sparse matrix
this function returns:
- res : a sparse matrix of `[A'|B]` , `A'` has only the non zero columns of `A`
- usedDofs : the indices of the columns present in `A'`
"""
return self.CleanEmptyColumns(self.ToSparseFull())
[docs]
def ToSparseFull(self):
"""
Function to convert to the constraints system Ax=B to a sparse matrix
this function returns:
res : a sparse matrix of [A|B] ,
"""
mat = coo_matrix((self.vals, (self.rows, self.cols)), shape=((self.numberOfEquations, self.nbdof + 1)), copy=False)
return mat
[docs]
def CleanEmptyColumns(self, mat):
"""
Function to clean empty columns for a matrix,
this function returns:
res : a sparse matrix with only non zero columns of mat (the last
column is always present, event if it is full of zeros)
usedDofs : the indices of the columns prensent in res
"""
mat = mat.tocoo()
mat.eliminate_zeros()
nbdof = mat.shape[1] - 1
usedDofs = np.unique(mat.col).astype(dtype=MuscatIndex)
if len(usedDofs) == 0:
return coo_matrix(([], ([], [])), shape=((0, 0)), copy=True), usedDofs
if usedDofs[-1] != nbdof:
# we add the independent term in the case is all zeros (last column = 0)
usedDofs = np.append(usedDofs, [nbdof])
nbUsedDofs = len(usedDofs)
# map fron global to the reduced matrix
mapGToR = np.zeros(mat.shape[1], dtype=MuscatIndex)
mapGToR[usedDofs] = range(nbUsedDofs)
col = mapGToR[mat.col]
res = coo_matrix((mat.data, (mat.row, col)), shape=((mat.shape[0], nbUsedDofs)), copy=True)
return res, usedDofs
[docs]
def GetNumberOfConstraints(self):
return self.numberOfEquations
[docs]
def GetNumberOfDofsOnOriginalSystem(self):
return self.nbdof
[docs]
def GetCleanConstraintBase(self):
"""
This function compute the QR decomposition of the augmented system [A|b]
Il will use the self.tol value to detect redundant constraints (rank revealing)
this functions returns:
res : a parse matrix q' with only the relevant part of the QR
decomposition.
Q' constraint only the non zero columns of [A|b] and always the last
column (the b)
usedDofs : the indices of the columns present in Q
"""
M, usedDofs = self.ToSparse()
Mcsc = M.tocsc()
am = abs(Mcsc[:, :-1])
am.data[am.data != 0] = 1
adjacency_matrix = am.dot(am.T)
del am
Mcsr = M.tocsr()
ncc, cc = connected_components(adjacency_matrix, directed=False)
res_vals = []
res_rows = []
res_cols = []
rankcpt = 0
numberOfDofPerEquation = adjacency_matrix.diagonal()
if ncc == M.shape[0] and np.all(numberOfDofPerEquation == 1):
# shortcut to check if we have independent 1 dof equations
# this means all dofs are slaves
Info(f"using shortcut to treat KCs")
res_vals = M.data
res_rows = M.row
res_cols = M.col
rankcpt = M.shape[0]
slaves = np.arange(M.shape[1] - 1, dtype=MuscatIndex)
else:
LS = NativeEigenSolver.CEigenSolvers()
LS.SetTolerance(self.tol)
LS.SetSolverType("SPQR")
atonce = np.zeros(M.shape[0], dtype=bool)
slaves = []
for incc in range(ncc):
mask = cc == incc
nl = np.count_nonzero(mask)
if nl == 1 and np.sum(numberOfDofPerEquation[mask]) == 1:
# we have one independent equation
np.logical_or(atonce, mask, out=atonce)
continue
# extraction of the equation to treat (rows) and filter the matrix
# to have only the used dofs (cols)
subsubM, subusedDofs = self.CleanEmptyColumns(Mcsr[mask, :])
LS.SetOp(subsubM.T)
rank = LS.GetSPQRRank()
slaves.extend(subusedDofs[0:rank])
Q = LS.GetSPQR_Q()
realQ = Q.tocsr()[:, 0:rank].tocoo().T
res_vals.extend(realQ.data)
res_rows.extend(realQ.row + rankcpt)
res_cols.extend(subusedDofs[realQ.col])
rankcpt += rank
# we treat all the single equation constraint at once
nb_single_constraint = np.count_nonzero(atonce)
if nb_single_constraint:
subM = Mcsr[atonce, :].tocoo()
res_vals.extend(subM.data)
res_rows.extend(subM.row + rankcpt)
res_cols.extend(subM.col)
rankcpt += nb_single_constraint
cols = np.unique(subM.col)
if cols[-1] == M.shape[1] - 1:
slaves.extend(cols[0:-1])
else:
slaves.extend(cols)
# normalisation
res = coo_matrix((res_vals, (res_rows, res_cols)), dtype=MuscatFloat, shape=(rankcpt, len(usedDofs)))
rescsr = res.tocsr()
norm = 1 / splinalg.norm(rescsr[:, :-1], axis=1)
res.data *= norm[res.row]
return res, usedDofs, slaves
# ----------------------- External API ------------------
[docs]
def SetConstraintsMethod(self, method: str):
MethodClass = methodFactory.get(method.lower(), None)
if MethodClass is None: # pragma: no cover
raise (Exception(f"Method '{method}' to available: possible options are {list(methodFactory.keys())}"))
self.method = MethodClass()
self.methodName = method
[docs]
def UpdateCSystem(self):
self.method.UpdateSystem(self)
[docs]
def GetNumberOfDofsOnCSystem(self):
return self.method.GetNumberOfDofs()
[docs]
def GetCOp(self):
return self.method.GetCOp(self.op)
[docs]
def matvec(self, arg):
return self.method.matvec(self.op, arg)
[docs]
def rmatvec(self, arg):
return self.method.rmatvec(self.op, arg)
[docs]
def GetCRhs(self, rhs=None):
if rhs is not None:
self.rhs = rhs
return self.method.GetCRhs(self.op, self.rhs)
[docs]
def RestoreSolution(self, arg):
return self.method.RestoreSolution(self.op, self.rhs, arg)
[docs]
def RestrictSolution(self, arg):
return self.method.RestrictSolution(self.op, self.rhs, arg)
def __str__(self):
res = "Constraints Holder: \n"
res = f" Number of contraints: {self.GetNumberOfConstraints()}\n"
res += str(self.method)
return res
[docs]
def ExpandMatrix(op, mattoglobal, nbdofs, treatrows=True, treatcols=True):
op = op.tocoo()
data = op.data
if treatrows:
rows = mattoglobal[op.row]
rs = nbdofs[0]
else:
rows = op.row
rs = op.shape[0]
if treatcols:
cols = mattoglobal[op.col]
cs = nbdofs[1]
else:
cols = op.col
cs = op.shape[1]
return sparse.coo_matrix((data, (rows, cols)), shape=(rs, cs)).tocsr()
# --------------------------------------------------------
[docs]
class Ainsworth:
"""
Method from https://doi.org/10.1016/S0045-7825(01)00236-5
Because the CleanConstraint(C) matrix is orthonormal then C.dot(C.T) = I
this make the P Q R expression simpler
"""
def __init__(self):
super().__init__()
self.P = None
self.Q = None
self.R = None
self.g = None
self.mattoglobal = None
[docs]
def UpdateSystem(self, CH):
self.nbdof = CH.GetNumberOfDofsOnOriginalSystem()
mat, self.mattoglobal, slaves = CH.GetCleanConstraintBase()
matcsc = mat.tocsc()
C = matcsc[:, 0:-1]
nbdofs = (self.nbdof, self.nbdof)
self.g = matcsc[:, -1]
CCT_1 = splinalg.inv(C.dot(C.T))
r = C.T.dot(CCT_1)
self.R = ExpandMatrix(r, self.mattoglobal, nbdofs, treatcols=False).tocsc()
p = r.dot(C)
self.P = ExpandMatrix(p, self.mattoglobal, nbdofs).tocsc()
self.Q = (sparse.eye(self.nbdof) - self.P).tocsc()
[docs]
def GetCOp(self, op):
return self.P + self.Q.T.dot(op.dot(self.Q))
[docs]
def GetCRhs(self, op, rhs):
a = self.R.dot(self.g).toarray()[:, 0]
b = self.Q.T.dot(rhs - op.dot(self.R.dot(self.g).toarray()[:, 0]))
return (a + b).flatten()
[docs]
def GetNumberOfDofs(self):
return self.nbdof
[docs]
def matvec(self, op, arg):
return self.P.dot(arg) + self.Q.T.dot(op.dot(self.Q.dot(arg)))
[docs]
def RestoreSolution(self, op, rhs, arg):
return arg
[docs]
def RestrictSolution(self, op, rhs, arg):
return arg
methodFactory["Ainsworth".lower()] = Ainsworth
[docs]
class Penalisation:
def __init__(self):
super().__init__()
self.factor = 1e8
self.maxdiag = 1.0
[docs]
def UpdateSystem(self, CH):
self.nbdof = CH.GetNumberOfDofsOnOriginalSystem()
mat, self.mattoglobal, slaves = CH.GetCleanConstraintBase()
matcsc = mat.tocsc()
A = matcsc[:, 0:-1]
b = matcsc[:, -1]
op = A.T.dot(A).tocoo() # A^T A
data = op.data
rows = self.mattoglobal[op.row]
cols = self.mattoglobal[op.col]
self.penalOp = sparse.coo_matrix((data, (rows, cols)), shape=(self.nbdof, self.nbdof)).tocsr()
self.penalRhs = np.zeros(self.nbdof)
rhs = A.T.dot(b).toarray().ravel()
self.penalRhs[self.mattoglobal[:-1]] = rhs
[docs]
def GetCOp(self, op):
self.maxdiag = max(op.diagonal())
return op + self.maxdiag * self.factor * self.penalOp
[docs]
def GetCRhs(self, op, rhs):
return rhs + self.maxdiag * self.factor * self.penalRhs
[docs]
def GetNumberOfDofs(self):
return self.nbdof
[docs]
def matvec(self, op, arg):
return op.dot(arg) + self.factor * self.maxdiag * self.penalOp.dot(arg)
[docs]
def RestoreSolution(self, op, rhs, arg):
return arg
[docs]
def RestrictSolution(self, op, rhs, arg):
return arg
methodFactory["Penalisation".lower()] = Penalisation
[docs]
class Lagrange:
def __init__(self):
super().__init__()
def __str__(self):
res = "Lagrange method:\n"
return res
[docs]
def UpdateSystem(self, CH):
mat, self.mattoglobal, slaves = CH.GetCleanConstraintBase()
self.mat = mat.tocsr()[:, 0:-1].tocoo()
self.rhs = mat.tocsr()[:, -1].toarray().ravel()
self.nbdof = CH.GetNumberOfDofsOnOriginalSystem()
[docs]
def GetCOp(self, op):
nbdof = op.shape[0]
nbcons = self.mat.shape[0]
op = op.tocoo()
data = np.hstack((op.data, self.mat.data, self.mat.data))
rows = np.hstack((op.row, self.mat.row + nbdof, self.mattoglobal[self.mat.col]))
cols = np.hstack((op.col, self.mattoglobal[self.mat.col], self.mat.row + nbdof))
return sparse.coo_matrix((data, (rows, cols)), shape=(nbdof + nbcons, nbdof + nbcons))
[docs]
def GetCRhs(self, op, rhs):
return np.hstack((rhs.ravel(), self.rhs.ravel()))
[docs]
def GetNumberOfDofs(self):
nbcons = self.mat.shape[0]
return self.nbdof + nbcons
[docs]
def matvec(self, op, arg):
u = np.zeros(self.GetNumberOfDofs())
u[0 : self.nbdof] += op.dot(arg[0 : self.nbdof])
u[self.nbdof :] += self.mat.dot(arg[self.mattoglobal[:-1]])
u[self.mattoglobal[:-1]] += self.mat.T.dot(arg[self.nbdof :])
return u
[docs]
def RestoreSolution(self, op, rhs, arg):
return arg[0 : self.nbdof]
[docs]
def RestrictSolution(self, op, rhs, arg):
res = np.zeros(self.GetNumberOfDofs())
res[0 : self.nbdof] = arg
res[self.nbdof :] = self.mat.dot((rhs - op.dot(arg))[self.mattoglobal[:-1]])
return res
methodFactory["Lagrange".lower()] = Lagrange
[docs]
class Projection:
def __init__(self):
super().__init__()
self.slaves = None
self.masters = None
self.X = None
self.D = None
def __str__(self):
res = "Projection method:\n"
res += "Slaves " + str(self.slaves) + "\n"
res += "Masters " + str(self.masters) + "\n"
if self.X is not None:
res += "X " + str(self.X.toarray()) + "\n"
else:
res += "X None\n"
res += "D " + str(self.D) + "\n"
return res
[docs]
def GetNumberOfDofs(self)-> MuscatIndex:
return self.nbMaster
[docs]
def matvec(self, op, arg):
return self.X.T.dot(op.dot(self.X.dot(arg)))
[docs]
def GetCOp(self, op):
return self.X.T.dot(op.dot(self.X)).tocsc() # X^T*Op*X
[docs]
def GetCRhs(self, op, rhs) -> np.ndarray:
return self.X.T.dot((rhs) - op.dot(self.D))
[docs]
def RestoreSolution(self, op, rhs, arg) -> np.ndarray:
return self.X.dot(arg) + self.D
[docs]
def RestrictSolution(self, op, rhs, arg):
return arg[self.masters]
[docs]
def UpdateSystem(self, CH: ConstraintsHolder):
Info("Compute Factorisation")
nbdof = CH.GetNumberOfDofsOnOriginalSystem()
mat, self.mattoglobal, slaves = CH.GetCleanConstraintBase()
slavesInMat = np.require(slaves, dtype=MuscatIndex)
slaveIds = self.mattoglobal[slavesInMat]
mat = mat.tocsc()
self.slaves = np.require(slaveIds, dtype=MuscatIndex)
masterMatMask = np.ones(mat.shape[1] - 1, dtype=bool)
masterMatMask[slavesInMat] = False
masterInMat = np.where(masterMatMask)[0]
masterMask = np.ones((nbdof), dtype=bool)
masterMask[self.slaves] = False
masterMask[self.mattoglobal[masterInMat]] = False
#we need to put the constrained master first in the list, then the free masters
self.masters = np.hstack((self.mattoglobal[masterInMat],np.where( masterMask)[0]))
self.nbMaster = len(self.masters)
Ms = mat[:, slavesInMat]
Mm = mat[:, masterInMat]
self.Ms_1 = splinalg.splu(Ms)
Ms_1Mm = sparse.coo_matrix(self.Ms_1.solve(Mm.toarray()))
Xdata = np.hstack((np.ones(self.nbMaster, dtype=MuscatFloat), -Ms_1Mm.data))
Xrow = np.hstack((self.masters, slaveIds[Ms_1Mm.row]))
Xcol = np.hstack((np.arange(self.nbMaster, dtype=MuscatIndex), Ms_1Mm.col))
self.X = sparse.coo_matrix((Xdata, (Xrow, Xcol)), shape=(nbdof, self.nbMaster))
self.D = np.zeros((nbdof), dtype=MuscatFloat)
v = mat[:, -1].toarray()
self.D[self.slaves] = self.Ms_1.solve(v).ravel()
methodFactory["Projection".lower()] = Projection
[docs]
def CheckIntegrityMasterSlaveDetection(GUI: bool = False):
def CheckSlaveDofs(A, b, checkSlaves):
CH = ConstraintsHolder()
CH.AddEquationsFromOperatorAAndb(A, b)
CH.SetNumberOfDofs(A.shape[1])
CH.Compact()
mat, mattoglobal, slaves = CH.GetCleanConstraintBase()
slaves = np.unique(slaves)
checkSlaves = np.unique(checkSlaves)
print(A.todense())
print(checkSlaves, slaves)
np.testing.assert_equal(slaves, checkSlaves)
A = coo_matrix(np.eye(3))
b = np.zeros(3)
slaves = [0, 1, 2]
CheckSlaveDofs(A, b, slaves)
A = coo_matrix(np.array([[1, 0, 0], [0, 1, 1]]))
b = np.zeros(2)
slaves = [0, 1]
CheckSlaveDofs(A, b, slaves)
A = coo_matrix(np.array([[1, 0, 0, 0], [0, 1, 1, 0], [0, 0, 0, 1]]))
b = np.zeros(3)
slaves = [0, 1, 3]
CheckSlaveDofs(A, b, slaves)
return "ok"
[docs]
def CheckIntegrityTTC(ttc, GUI=False):
print(f"------------------------------- {ttc} ------------------------------")
CH = ConstraintsHolder()
CH.SetConstraintsMethod(ttc)
CH.SetNumberOfDofs(4)
# Reference solution
refsol = np.array([0, 1, 2, 3], dtype=MuscatFloat)
sys = 100 * np.array([[1, -1, 0, 0], [-1, 1, 0, 0], [0, 0, 1, -1], [0, 0, -1, 1]], dtype=MuscatFloat)
rhs = np.array([0, 0, 0, 0], dtype=MuscatFloat)
# Dirichlet zero left
CH.AddEquation([1, 0, 0, 0], 0)
CH.AddEquationsFromIJV([0], [0], [1])
# kinematic relation
CH.AddEquation([0.0, -1.0, 1.0, 0.0], 1.0)
# to test the elimination of redundant equations
CH.AddEquation([0.0, -2.0, 2.0, 0.0], 2.0)
# Dirichlet right side
CH.AddEquation([0.0, 0.0, 0.0, 1.0], 3.0)
CH.AddEquationSparse([3], [1.0], 3.0)
CH.AddConstant(0)
CH.Compact()
np.set_printoptions(threshold=np.inf, linewidth=np.inf)
print("System to Solve ")
print("Number Of constraints: ", CH.GetNumberOfConstraints())
print("K ")
print(sys)
print("F ")
print(rhs)
print("Under Constraints ")
Ab = CH.ToSparse()[0].toarray()
A = Ab[:, :-1]
b = Ab[:, -1]
print("A ")
print(A)
print("b ")
print(b)
print("Reference solution ")
print(refsol)
CH.UpdateCSystem()
Mv = CH.GetCleanConstraintBase()[0].toarray()
M = Mv[:, :-1]
v = Mv[:, -1]
print("Clean Contraints Base ")
print("M ")
print(M)
print("v ")
print(v)
np.testing.assert_equal(len(v), 3)
CH.SetOp(sys, rhs)
Kc = CH.GetCOp().tocsc()
Fc = CH.GetCRhs()
print("\nConstrained operator : \n", Kc.toarray())
conditionNumber = np.linalg.cond(Kc.toarray(), p="fro")
print("\nCondition number of Constrained operator : \n", conditionNumber)
print("\nConstrained rhs : \n", Fc)
Uc = CH.RestrictSolution(refsol)
print("\nRestricted solution (from refsol) : \n", Uc)
bb, au = CH.Residual(refsol)
np.testing.assert_equal(bb, au)
contrainedSolutionError = np.linalg.norm(Kc.dot(Uc) - Fc) / np.linalg.norm(Fc)
print("\nContraines solution error (|Kc(Uc)-Fc|/|Fc|) : \n", contrainedSolutionError)
if contrainedSolutionError > 1e-8: # pragma: no cover
raise Exception("Error in CheckIntegrity")
# try to solve the constrained system Using the matvec
# we use the gmres to correctrly treat non defined matrices
from scipy.sparse import linalg
contrainedSysNbDofs = CH.GetNumberOfDofsOnCSystem()
op = linalg.LinearOperator((contrainedSysNbDofs, contrainedSysNbDofs), matvec=CH.matvec, rmatvec=CH.rmatvec, dtype=MuscatFloat)
# the penalisation is not a good method with a iterative solver so we need
# a simple precoditioning to help the solver
diag = Kc.diagonal()
diag[diag == 0] = 1.0
M = sparse.dia_matrix((1.0 / diag, 0), shape=Kc.shape)
sol, gmres_info = linalg.gmres(op, Fc, atol=1e-15, rtol=1e-15, M=M)
print(f"gmres_info = {gmres_info}")
print("---------- Solution using gmres -------------- \n", sol)
errorA = np.linalg.norm(Kc.dot(sol) - Fc) / np.linalg.norm(Fc)
print("Error on the computed solution (using matricexs) |Kc(Uc)-Fc|/|Fc| \n", errorA)
if errorA > 1e-7: # pragma: no cover
raise Exception("Error in CheckIntegrity")
errorB = np.linalg.norm(CH.matvec(sol) - Fc) / np.linalg.norm(Fc)
print("Error on the computed solution (using matvec) | Op(Uc)-Fc |/|Fc| \n", errorB)
if errorB > 1e-7: # pragma: no cover
raise Exception("Error in CheckIntegrity")
computed_solution = CH.RestoreSolution(sol)
print("--------------------------")
print("Reference solution \n", refsol)
print("Computed solution \n", computed_solution)
errorC = np.linalg.norm(computed_solution - refsol) / np.linalg.norm(refsol)
print("relative error |ref_sol - computed_sol|/|ref_sol|", errorC)
if errorC > 1e-8: # pragma: no cover
raise Exception("Error in CheckIntegrity")
print("---------- Solution using spsolve -------------- \n", sol)
sol = sparse.linalg.spsolve(Kc, Fc)
# sol = np.linalg.inv(K).dot(CH.GetCRhs())
dfsol = CH.RestoreSolution(sol)
print("Final direct solution", dfsol)
errorD = np.linalg.norm(dfsol - refsol) / np.linalg.norm(refsol)
print("Error on the solution |ref_sol - computed_sol|/|ref_sol| \n ", errorD)
if errorD > 1e-8: # pragma: no cover
raise Exception("Error in CheckIntegrity")
print(CH)
print("------------------- " + ttc + " DONE -----------------")
return "OK"
[docs]
def CheckIntegrity(GUI: bool = False):
from Muscat.Helpers.CheckTools import RunListOfCheckIntegrities, CheckIntegrityPresent
toTest = [CheckIntegrityMasterSlaveDetection]
typeToCheck = list(methodFactory.keys())
import functools
for ttc in typeToCheck:
toTest.append(functools.partial(CheckIntegrityTTC, ttc))
return RunListOfCheckIntegrities(toTest)
if __name__ == "__main__":
print(CheckIntegrity(GUI=True)) # pragma: no cover