Skip to content
Snippets Groups Projects
Commit 2c7e7c5e authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

removing lagrange stuff

parent 8ddbd0fa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 7 additions and 564 deletions
bin/args.py
+ 0
− 9
View file @ 2c7e7c5e
......@@ -50,15 +50,6 @@ def parse():
help="set constant number NAME=VALUE")
parser.add_argument("-quiet", action='store_true',
help="should I be quiet?")
parser.add_argument("-non_holomorphic", action='store_true',
help='start a fixed-point loop for non-holomorphic '+
'functions')
parser.add_argument("-nhMaxIt", type=int, default=5,
help='maximum iteration of non-holomorphic solver')
parser.add_argument("-nhMinRratio", type=float, default=10,
help='minimum ratio between non-holomorphic and ' +
'holomorphic residuals for ending fixed-point ' +
'iterations')
# Done
return parser.parse_args()
......
bin/cim.py
+ 7
− 18
View file @ 2c7e7c5e
......@@ -13,12 +13,11 @@ try:
except ImportError:
from mpiseq import MPI # If not available use a dummy mpiseq
from solver import GetDPWave
from beyn import Beyn
from fixed_point import FixedPoint
import args as args
import numpy as np
import time as timer
from solver import GetDPWave
from beyn import Beyn
import args as args
import numpy as np
import time as timer
# Start timer
tStart = timer.time()
......@@ -46,18 +45,8 @@ beyn = Beyn(operator, arg.origin, arg.radius,
not arg.quiet)
# Compute
if(not arg.non_holomorphic):
if(not arg.quiet): beyn.display()
l, V, r = beyn.simple()
else:
fixed = FixedPoint(beyn, arg.nhMaxIt, arg.nhMinRratio, not arg.quiet)
# This variable ('fixed') is meant to be called
if(verbose): # elsewhere (i.e. Python[] call in GetDP) thanks
fixed.display() # to the python interpreter.
beyn.display() # Please don't change its name!
l, V, r = fixed.solve()
if(not arg.quiet): beyn.display()
l, V, r = beyn.simple()
# Display the computed eigenvalues
if(verbose):
......
bin/fixed_point.py deleted 100644 → 0
+ 0
− 196
View file @ 8ddbd0fa
"""
cim.py, a non-linear eigenvalue solver.
Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
Institut fuer Theorie Elektromagnetischer Felder (TEMF),
Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
"""
TODO: cluster identification / can a spline be holomorphic?
"""
try:
from mpi4py import MPI # Try to import mpi4pi
except ImportError:
from mpiseq import MPI # If not available use a dummy mpiseq
import sys as sys
import numpy as np
class FixedPoint:
"""A fixed-point solver for handling non-holomorphic operator
This solver allows to approximate a non-holomorphic operator by an
holomorphic one. This approximation is refined by a fixed-point iteration,
until the residual of the original (non-holomorphic) problem reaches the
residual of the corrected (holomorphic) problem. The holomorphic correction
is preformed by calling a python handler from GetDP (using the Python[]
function), which is responsible for implementing the holomorphic correction.
The handler can communicate with this fixed-point solver. By convention,
this solver is refred to as 'fixed', and is accessible from the handler.
An example is provided in: ./example/maxwell_sibc_lagrange/
"""
def __init__(self, solver, maxIt=5, minNHHRatio=10, verbose=True):
"""Instantiates a new FixedPoint non-holomorphic solver
Keyword arguments:
solver -- the solver used to compute holomorphic solutions
maxIt -- the maximum number of fixed-point iteration (optional)
minNHHRatio -- the solver stopping criterion (see below, optional)
verbose -- should I be verbose? (optional)
What is minNHHRatio?
The FixedPoint solver computes two residuals for each solution it finds:
1. rH, computed with respect to the corrected, holomorphic, operator;
2. rNH, computed with respect to the original, non-holomorphic, operator
Ultimately, the FixedPoint solver tries to impose:
rNH = rH,
for all eigenpairs. Thus, we use max(rNH/rH) as an indicator of the
overall accuracy. Therefore, the FixedPoint solver will stop as soon as:
max(rNH/rH) > minNHHRatio,
or if the number of iteration exceeds maxIt. In other words, the closer
is minNHHRatio to 1, the sharper the solution. Please note that it does
not make sense to have minNHHRatio < 1. That is, the corrected
holomorphic solution cannot be more precise that the original
non-holomorphic one.
"""
# Save
self.__solver = solver
self.__maxIt = maxIt
self.__minNHHRatio = minNHHRatio
self.__myProc = MPI.COMM_WORLD.Get_rank()
self.__verbose = verbose and (self.__myProc == 0)
# State
self.__cluster = []
self.__holomorphic = True
def solve(self):
"""Solves the eigenvalue problem by a fixed-point iteration
It returns the computed eigenvalues, eigenvectors and the norm of the
residual associated to the original, non-holomorphic, operator
"""
# Pretty copy
verbose = self.__verbose
# (Re)initialize
self.__cluster = []
self.__holomorphic = True
# First run with corrector initial guess
if(verbose): print ">> Non-holomorphic fixed-point solver: initial run"
l, V, rH = self.__solver.simple() # Holomorphic solution
rNH = self.__nhNorm(l, V) # Non-holomorphic residual
maxRNHH = self.__maxRNHH(rNH, rH) # Maximim rNH/rH
if(verbose): print rH
if(verbose): print rNH
if(verbose): print ">> Initial max(rNH/rH): " + format(maxRNHH, ".2e")
# Fixed-point iteration
it = 0
if(verbose): print ">> Non-holomorphic fixed-point solver: fixed-point"
while(maxRNHH > self.__minNHHRatio and it != self.__maxIt):
self.__cluster = self.__findCluster(l) # Clusters for correction
l, V, rH = self.__solver.simple() # Solve corrected system
rNH = self.__nhNorm(l, V) # Norm of NH residual
maxRNHH = self.__maxRNHH(rNH, rH) # Maximim rNH/rH
if(verbose): print rH
if(verbose): print rNH
if(verbose): print ">> Fixed-point iteration " + str(it+1),
if(verbose): print "done: " + format(maxRNHH, ".2e")
it = it + 1 # One more iteration :-)
# Warn if too many fixed-point iterations
if(it == self.__maxIt):
if(verbose): ">> Maximum number of fixed-point iterations reached"
# Return last solution with the norm of the non-holomorphic residual
return l, V, rNH
def cluster(self):
"""Returns the eigenvalues (by cluster) found by the last iteration
This method returns the eigenvalues "by cluster". That is, it some
eigenvalues are too close from each other, only one will be considered
"""
return self.__cluster
def needNonHolomorphic(self):
"""This method returns True if this FixedPoint needs the true,
non-holomorphic, operator and not its holomorphic approximation"""
return not self.__holomorphic
def display(self):
"""Prints this FixedPoint solver parameters"""
print "Fixed-point non-holomorphic correction "
print "---------------------------------------"
print ">> Maximum number of interation: "+ " " + str(self.__maxIt)
print ">> Minimum NH-H residual ratio: ",
print format(self.__minNHHRatio, ".2g")
print "---------------------------------------"
def __nhNorm(self, l, V):
"""Returns the norm of the non-holomorphic residual
for all eigenpairs (l, V)"""
operator = self.__solver.operator() # Get operator
nFound = l.shape[0] # Number of found eigenpair
nhNorm = np.empty(nFound) # Non-holomorphic residual
self.__holomorphic = False # Non-holomorphic (NH) operator
for n in range(nFound): # Iterate on solutions
operator.apply(V[:, n], l[n]) # -> Apply solution on NH operator
nhRes = operator.solution(0) # -> residual
nhNorm[n] = np.linalg.norm(nhRes) # -> norm
self.__holomorphic = True # Back to holomorphic operator
return nhNorm # Done
@staticmethod
def __findCluster(data):
"""Retruns the cluster contained in the given *sorted* list"""
rank = MPI.COMM_WORLD.Get_rank()
# Get clusters
cluster = []
delta = 0.01
i = 1
while i < len(data):
j = i
sub = [data[i-1]]
while j < len(data):
d = (data[i-1] - data[j]) / data[j]
if abs(d.real) < delta and abs(d.imag) < delta:
sub.append(data[j])
j = j + 1
else:
break
cluster.append(sub)
i = j + 1
# Mean of each cluster
out = []
for i in range(len(cluster)):
out.append(np.mean(cluster.pop()))
out.reverse()
if(rank == 0): print "Cluster:",
if(rank == 0): print out
return out
# return [data[0], data[3], data[6], data[11]]
@staticmethod
def __maxRNHH(rNH, rH):
""" Returns the maximal rNH/rH"""
return max([a/b for a, b in zip(rNH, rH)])
bin/interpolator.py deleted 100644 → 0
+ 0
− 33
View file @ 8ddbd0fa
"""
cim.py, a non-linear eigenvalue solver.
Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
Institut fuer Theorie Elektromagnetischer Felder (TEMF),
Technische Universitaet Darmstadt.
See the LICENSE.txt and README.md for more license and copyright information.
"""
from numpy import poly1d
from scipy import interpolate
from abc import ABCMeta, abstractmethod
class Interpolator:
"""A generic class for interpolation"""
@abstractmethod
def __init__(self):
pass
def at(self, x):
"""Returns the value of this Interpolation at a given point"""
return self._p(x) # We use poly1d from scipy
__metaclass__ = ABCMeta
class Lagrange(Interpolator):
"""A class for Lagrange interpolation"""
def __init__(self, xdata, ydata):
"""Instantiates a new Interpolator, such that:
p(xdata) = ydata for all xdata, where p is build with a Lagrange basis
"""
self._p = interpolate.lagrange(xdata, ydata)
example/maxwell_sibc_lagrange/Makefile deleted 100644 → 0
+ 0
− 21
View file @ 8ddbd0fa
BIN=../../bin
sphere: init cim
init:
gmsh sphere.geo -3
cim:
# python ${BIN}/cim.py sphere.pro sphere.msh Solve 1e10+1e9j 4e9 -nodes 50 -lStart 20 -lStep 30 -setnumber isSC 0 -non_holomorphic
python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0 -non_holomorphic
#7.474355560164e+09 + 7.742536715699e+08j
# python ${BIN}/cim.py sphere.pro sphere.msh Solve 1.0e10+1.1e9j 1e9 -nodes 10 -lStart 6 -setnumber isSC 0 -non_holomorphic
#1.050336615299e+10 + 1.186115651023e+09j
clean:
rm -f *.pos
rm -f *.msh
rm -f *.pre
rm -f *.res
example/maxwell_sibc_lagrange/cimParameters.pro deleted 100644 → 0
+ 0
− 26
View file @ 8ddbd0fa
/////////////////////////////////////////////////////////////////
// Parameters for contour integral method. //
// The master file should inclue "cimResolution.pro" later on. //
/////////////////////////////////////////////////////////////////
Function{
// Physical data //
DefineConstant[angularFreqRe = 1, // [rad/s]
angularFreqIm = 0]; // [rad/s]
Puls[] = Complex[angularFreqRe,
angularFreqIm]; // [rad/m]
// Algebraic data //
DefineConstant[nRHS = 1]; // Number of RHS for this run
For i In {0:nRHS-1}
DefineConstant[x~{i}() = {}, // Solution
b~{i}() = {}]; // Right hand side
EndFor
// Control data //
DefineConstant[doInit = 0, // Should I initialize some stuff?
doSolve = 0, // Should I solve Ax = b?
doPostpro = 0, // Should I only create a view for x()?
doApply = 0, // Should I only apply x(): x <- Ax?
fileName = "eig.pos"]; // Postpro file name
}
example/maxwell_sibc_lagrange/cimResolution.pro deleted 100644 → 0
+ 0
− 47
View file @ 8ddbd0fa
/////////////////////////////////////////////////////////////////////
// Resolution for contour integral method. //
// Warning, the master file should include "cimParameters.pro", //
// and should also define "FormulationCIM" the formulation to use. //
/////////////////////////////////////////////////////////////////////
Resolution{
{ Name Solve;
System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
Operation{
MPI_SetCommSelf;
If(!doInit)
Evaluate[Python[]{"real_corrector_init.py"}];
EndIf
Generate[A];
If(doInit)
CopySolution[A, x~{0}()];
EndIf
If(doSolve)
// Full solve for first RHS
CopyRightHandSide[b~{0}(), A];
Solve[A];
CopySolution[A, x~{0}()];
// SolveAgain for other RHSs
For i In {1:nRHS-1}
CopyRightHandSide[b~{i}(), A];
SolveAgain[A];
CopySolution[A, x~{i}()];
EndFor
EndIf
If(doApply)
CopySolution[x~{0}(), A];
Apply[A];
CopySolution[A, x~{0}()];
EndIf
If(doPostpro)
CopySolution[x~{0}(), A];
PostOperation[Post];
EndIf
}
}
}
example/maxwell_sibc_lagrange/real_corrector_handl.py deleted 100644 → 0
+ 0
− 22
View file @ 8ddbd0fa
"""Real-part operator corrector
Approximates the real-part operator by an holomorphic function, constructed with
a Lagrange polynomial defined at some correction points. Please make sure that
GetDP calls real_corrector_init.py before calling this file.
The naming convention below must be followed:
-- a FixedPoint object named 'fixed' is available and prepared by cim.py
-- a Lagrange Interpolator object named 'lagrange' coming from initialization
-- the variable 'output' is read back by GetDP as the Pyhton[] call output
"""
# Input angular frequency
omega = input[0]
# Compute holomorphic approximation of the real part operator
if(fixed.needNonHolomorphic()):
output = np.real(omega) # If need non-holomorphic, do it
elif(not fixed.cluster()):
output = omega # If no cluster: Re(omega) ~ omega
else:
output = (omega + lagrange.at(omega)) / 2 # Else: correct
example/maxwell_sibc_lagrange/real_corrector_init.py deleted 100644 → 0
+ 0
− 13
View file @ 8ddbd0fa
"""Real-part operator corrector: initialization
See: real_corrector_handl.py for more information
"""
from interpolator import Lagrange
# Previous eigenvalue clusters
cluster = fixed.cluster()
# Initialize interpolator
if(cluster):
lagrange = Lagrange(cluster, np.conj(cluster).tolist())
example/maxwell_sibc_lagrange/sphere.dat deleted 100644 → 0
+ 0
− 22
View file @ 8ddbd0fa
// Geometry //
DefineConstant[R = { 100e-3, Name "Input/00Geometry/00Radius [m]" }];
// Mesh //
DefineConstant[md = { 10, Name "Input/01Mesh/00Density [Radius^-1]" }];
cl = R/md;
// Material //
DefineConstant[isSC = { 0,
Name "Input/02Material/00Is superconductor?",
GmshOption "Reset",
Choices {0 = "No", 1 = "Yes"} }];
If(isSC == 0)
DefineConstant[Sigma = { 1e00,
Name "Input/02Material/01Conductivity [S m^-1]" }];
EndIf
If(isSC == 1)
DefineConstant[Sigma = { 1e9,
Name "Input/02Material/01Conductivity [S m^-1]" }];
DefineConstant[Lambda = { 1e-8,
Name "Input/02Material/02London depth [m]" }];
EndIf
example/maxwell_sibc_lagrange/sphere.geo deleted 100644 → 0
+ 0
− 50
View file @ 8ddbd0fa
Include "sphere.dat";
Point(0) = { 0, 0, 0, cl};
Point(1) = {+R, 0, 0, cl};
Point(2) = { 0, +R, 0, cl};
Point(3) = {-R, 0, 0, cl};
Point(4) = { 0, -R, 0, cl};
Point(5) = { 0, 0, +R, cl};
Point(6) = { 0, 0, -R, cl};
Circle(01) = {1, 0, 2};
Circle(02) = {2, 0, 3};
Circle(03) = {3, 0, 4};
Circle(04) = {4, 0, 1};
Circle(05) = {5, 0, 1};
Circle(06) = {1, 0, 6};
Circle(07) = {6, 0, 3};
Circle(08) = {3, 0, 5};
Circle(09) = {2, 0, 6};
Circle(10) = {6, 0, 4};
Circle(11) = {4, 0, 5};
Circle(12) = {5, 0, 2};
Line Loop(1) = {1, -12, 5};
Line Loop(2) = {1, 9, -6};
Line Loop(3) = {9, 7, -2};
Line Loop(4) = {2, 8, 12};
Line Loop(5) = {11, 5, -4};
Line Loop(6) = {4, 6, 10};
Line Loop(7) = {10, -3, -7};
Line Loop(8) = {3, 11, -8};
Ruled Surface(1) = {1};
Ruled Surface(2) = {2};
Ruled Surface(3) = {3};
Ruled Surface(4) = {4};
Ruled Surface(5) = {5};
Ruled Surface(6) = {6};
Ruled Surface(7) = {7};
Ruled Surface(8) = {8};
Surface Loop(1) = {3, 2, 1, 4, 8, 7, 6, 5};
Volume(1) = {1};
Physical Volume(1) = {1};
Physical Surface(2) = {1, 2, 3, 4, 5, 6, 7, 8};
Physical Point(3) = {0};
example/maxwell_sibc_lagrange/sphere.pro deleted 100644 → 0
+ 0
− 107
View file @ 8ddbd0fa
Include "sphere.dat";
Group{
Omega = Region[1];
Bndry = Region[2];
}
Include "./cimParameters.pro"; // Functions for cim.py
Function{
// Imaginary Unit //
J[] = Complex[0, 1];
// Physical Constants //
C0 = 299792458; // [m/s]
Mu0 = 4*Pi*1e-7; // [H/m]
Eps0 = 1/(Mu0 * C0^2); // [F/m]
// Correction for real-part operator //
If(doInit)
PyPulsRe[] = 1;
EndIf
If(!doInit)
PyPulsRe[] = Python[Puls[]]{"real_corrector_handl.py"};
EndIf
// Conductivity (coming from sphere.dat) //
If(isSC == 1) // Superconductor
Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * PyPulsRe[]]); // [S/m]
EndIf
If(isSC == 0) // Normal conductor
Sigma[] = Sigma; // [S/m]
EndIf
// Leontovich SIBC (H-Formulation) //
SL[] = Sqrt[(J[] * PyPulsRe[] * Mu0) / Sigma[]]; // H-Formulation
SIbc[] = -1.0 / SL[]; // E-Formulation
}
Integration{
{ Name I1;
Case{
{ Type Gauss;
Case{
{ GeoElement Tetrahedron; NumberOfPoints 4; } // O1 * O1
{ GeoElement Triangle; NumberOfPoints 3; } // O1 * O1
{ GeoElement Line; NumberOfPoints 2; } // O1 * O1
}
}
}
}
}
Jacobian{
{ Name Jac;
Case{
{ Region Omega; Jacobian Vol; }
{ Region Bndry; Jacobian Sur; }
}
}
}
FunctionSpace{
{ Name HCurl; Type Form1;
BasisFunction{
{ Name sx; NameOfCoef ux; Function BF_Edge;
Support Region[{Omega, Bndry}]; Entity EdgesOf[All]; }
}
}
}
Formulation{
{ Name FormulationCIM; Type FemEquation;
Quantity{
{ Name e; Type Local; NameOfSpace HCurl; }
}
Equation{
// Maxwell
Galerkin{ [ 1/Mu0 * Dof{d e}, {d e}];
In Omega; Jacobian Jac; Integration I1; }
Galerkin{ [ -Puls[]^2 * Eps0 * Dof{ e}, { e}];
In Omega; Jacobian Jac; Integration I1; }
// IBC
Galerkin{ [-J[] * Puls[] * SIbc[] * Dof{ e}, { e}];
In Bndry; Jacobian Jac; Integration I1; }
}
}
}
Include "./cimResolution.pro"; // Resolution for cim.py
PostProcessing{
{ Name Post; NameOfFormulation FormulationCIM;
Quantity{
{ Name E; Value{ Term{ [{e}]; In Omega; Jacobian Jac; } } }
}
}
}
PostOperation{
{ Name Post; NameOfPostProcessing Post;
Operation{
Print[E, OnElementsOf Omega, File fileName];
}
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment