Skip to content
Snippets Groups Projects
Commit 31a1783c authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Adjust epsilon to 1e-8 in symbolicMinimization

parent f785b70a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Matlab-Programs/symMinimization.m
+ 15
− 10
View file @ 31a1783c
......@@ -28,6 +28,11 @@ f_minus(v1,v2,q1,q2,q3,q12,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2+2*(q1*b1
% ---- Fix parameters
% Epsilon used:
epsilon = 1e-8;
if ~exist('InputPath','var')
% third parameter does not exist, so default it to something
absPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs";
......@@ -52,7 +57,7 @@ Bmat = full(Bmat);
% --- TODO CHECK: assert if Q is orthotropic ???
% --- TODO CHECK: assert if Q is orthotropic ??? check ifq13=q31=q23=q32= 0 ?
if print_Input
......@@ -66,15 +71,16 @@ if print_Input
% Test: issymmetric(Qmat) does not work for float matrices?
% symmetric part 0.5*(Qmat+Qmat')
% anti-symmetric part 0.5*(Qmat-Qmat')
if norm(0.5*(Qmat-Qmat'),'fro') < 1e-10
if norm(0.5*(Qmat-Qmat'),'fro') < 1e-8
fprintf('Qmat (close to) symmetric \n')
norm(0.5*(Qmat-Qmat'),'fro') % TEST
else
fprintf('Qmat not symmetric \n')
end
% Check if B_eff is diagonal this is equivalent to b3 == 0
if abs(Bmat(3)) < 1e-10
if abs(Bmat(3)) < 1e-8
fprintf('B_eff is diagonal (b3 == 0) \n')
else
fprintf('B_eff is NOT diagonal (b3 != 0) \n')
......@@ -231,11 +237,10 @@ end
% SP_value = T_minus(idx) % not needed?
Matrix = sign*(SP'*SP);
if norm(double(Matrix-Minimizer),'fro') < eps %check is this sufficient here?
% both StationaryPoints correspond to the same
% (Matrix-)Minimizer
if norm(double(Matrix-Minimizer),'fro') < 1e-8 %check is this sufficient here?
% fprintf('both StationaryPoints correspond to the same(Matrix-)Minimizer')
else
% StationaryPoint corresponds to a different (Matrix-)Minimizer
% fprintf('StationaryPoint corresponds to a different (Matrix-)Minimizer')
MinimizerCount = MinimizerCount + 1;
end
end
......@@ -289,7 +294,7 @@ G = double(Minimizer);
if MinimizerCount == 1
% fprintf('Unique Minimzier')
% Check if Minimizer is axial or non-axial:
if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
Type = 3;
else % non-axial Minimizer
Type = 1;
......@@ -297,7 +302,7 @@ if MinimizerCount == 1
elseif MinimizerCount == 2
% fprintf('Two Minimziers')
% Check if Minimizer is axial or non-axial:
if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
Type = 3;
else % non-axial Minimizer
fprintf('ERROR: Two non-axial Minimizers cannot happen!')
......@@ -305,7 +310,7 @@ elseif MinimizerCount == 2
else
% fprintf('1-Parameter family of Minimziers')
% Check if Minimizer is axial or non-axial:
if (abs(angle-pi/2) < 1e-9 || abs(angle) < 1e-9) % axial Minimizer
if (abs(angle-pi/2) < 1e-8 || abs(angle) < 1e-8) % axial Minimizer
% fprintf('ERROR: axial Minimizers cannot happen for 1-Parameter Family!')
else % non-axial Minimizer
Type = 2;
......
src/CellScript.py
+ 7
− 259
View file @ 31a1783c
......@@ -9,168 +9,12 @@ import re
import matlab.engine
import time
from ClassifyMin import *
from HelperFunctions import *
# from scipy.io import loadmat #not Needed anymore?
import codecs
import sys
def ReadEffectiveQuantities(QFilePath = os.path.dirname(os.getcwd()) + '/outputs/QMatrix.txt', BFilePath = os.path.dirname(os.getcwd())+ '/outputs/BMatrix.txt'):
# Read Output Matrices (effective quantities)
# From Cell-Problem output Files : ../outputs/Qmatrix.txt , ../outputs/Bmatrix.txt
# -- Read Matrix Qhom
X = []
with codecs.open(path + '/outputs/QMatrix.txt', encoding='utf-8-sig') as f:
for line in f:
s = line.split()
X.append([float(s[i]) for i in range(len(s))])
Q = np.array([[X[0][2], X[1][2], X[2][2]],
[X[3][2], X[4][2], X[5][2]],
[X[6][2], X[7][2], X[8][2]] ])
# -- Read Beff (as Vector)
X = []
with codecs.open(path + '/outputs/BMatrix.txt', encoding='utf-8-sig') as f:
for line in f:
s = line.split()
X.append([float(s[i]) for i in range(len(s))])
B = np.array([X[0][2], X[1][2], X[2][2]])
return Q, B
def RunCellProblem(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
with open(InputFilePath, 'r') as file:
filedata = file.read()
filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
f = open(InputFilePath,'w')
f.write(filedata)
f.close()
# --- Run Cell-Problem
# Optional: Check Time
# t = time.time()
subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
capture_output=True, text=True)
# print('elapsed time:', time.time() - t)
# --------------------------------------------------------------------------------------
def Compare_Classification(alpha,beta,theta,gamma,mu1,rho1, InputFilePath = os.path.dirname(os.getcwd()) +"/inputs/computeMuGamma.parset"):
# ---------------------------------------------------------------
# Comparison of the analytical Classification 'ClassifyMin'
# and the symbolic Minimizatio + Classification 'symMinimization'
# ----------------------------------------------------------------
comparison_successful = True
eps = 1e-8
# 1. Substitute Input-Parameters for the Cell-Problem
with open(InputFilePath, 'r') as file:
filedata = file.read()
filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
filedata = re.sub('(?m)^alpha=.*','alpha='+str(alpha),filedata)
filedata = re.sub('(?m)^beta=.*','beta='+str(beta),filedata)
filedata = re.sub('(?m)^theta=.*','theta='+str(theta),filedata)
filedata = re.sub('(?m)^mu1=.*','mu1='+str(mu1),filedata)
filedata = re.sub('(?m)^rho1=.*','rho1='+str(rho1),filedata)
f = open(InputFilePath,'w')
f.write(filedata)
f.close()
# 2. --- Run Cell-Problem
print('Run Cell-Problem ...')
# Optional: Check Time
# t = time.time()
subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
capture_output=True, text=True)
# 3. --- Run Matlab symbolic minimization program: 'symMinimization'
eng = matlab.engine.start_matlab()
# s = eng.genpath(path + '/Matlab-Programs')
s = eng.genpath(path)
eng.addpath(s, nargout=0)
# print('current Matlab folder:', eng.pwd(nargout=1))
eng.cd('Matlab-Programs', nargout=0) #switch to Matlab-Programs folder
# print('current Matlab folder:', eng.pwd(nargout=1))
Inp = False
print('Run symbolic Minimization...')
G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4) #Name of program:symMinimization
# G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4) #Optional: add Path
G = np.asarray(G) #cast Matlab Outout to numpy array
# --- print Output ---
print('Minimizer G:')
print(G)
print('angle:', angle)
print('type:', type )
print('curvature:', kappa)
# 4. --- Read the effective quantities (output values of the Cell-Problem)
# Read Output Matrices (effective quantities)
print('Read effective quantities...')
Q, B = ReadEffectiveQuantities()
print('Q:', Q)
print('B:', B)
q1 = Q[0][0]
q2 = Q[1][1]
q3 = Q[2][2]
q12 = Q[0][1]
b1 = B[0]
b2 = B[1]
b3 = B[2]
# print("q1:", q1)
# print("q2:", q2)
# print("q3:", q3)
# print("q12:", q12)
# print("b1:", b1)
# print("b2:", b2)
# print("b3:", b3)
# --- Check Assumptions:
# Assumption of Classification-Lemma1.6: [b3 == 0] & [Q orthotropic]
# Check if b3 is close to zero
assert (b3 < eps), "ClassifyMin only defined for b3 == 0 !"
# Check if Q is orthotropic i.e. q13 = q31 = q23 = q32 == 0
assert(Q[2][0] < eps and Q[0][2] < eps and Q[1][2] < eps and Q[2][1] < eps), "Q is not orthotropic !"
# 5. --- Get output from the analytical Classification 'ClassifyMin'
G_ana, angle_ana, type_ana, kappa_ana = classifyMin(q1, q2, q3, q12, b1, b2)
# print('Minimizer G_ana:')
# print(G_ana)
# print('angle_ana:', angle_ana)
# print('type_ana:', type_ana )
# print('curvature_ana:', kappa_ana)
# 6. Compare
# print('DifferenceMatrix:', G_ana - G )
# print('MinimizerError (FrobeniusNorm):', np.linalg.norm(G_ana - G , 'fro') )
if np.linalg.norm(G_ana - G , 'fro') > eps :
comparison_successful = False
print('Difference between Minimizers is too large !')
if type != type_ana :
comparison_successful = False
print('Difference in type !')
if abs(angle-angle_ana) > eps :
comparison_successful = False
print('Difference in angle is too large!')
if abs(kappa-kappa_ana) > eps :
comparison_successful = False
print('Difference in curvature is too large!')
if comparison_successful:
print('Comparison successful')
return comparison_successful
# ----------------------------------------------------------------------------------------------------------------------------
# ----- Setup Paths -----
......@@ -188,97 +32,6 @@ print("OutputFilepath: ", OutputFilePath)
print("Path: ", path)
#1. Define Inputs Gamma-Array..
#2. for(i=0; i<length(array)) ..compute Q_hom, B_eff from Cell-Problem
#3
# matrix = np.loadtxt(path + 'Qmatrix.txt', usecols=range(3))
# print(matrix)
# Use Shell Commands:
# subprocess.run('ls', shell=True)
#
# #-------------------- PLOT OPTION -------------------------------------------
#
# Gamma_Values = np.linspace(0.01, 2.5, num=6) # TODO variable Input Parameters...alpha,beta...
# print('(Input) Gamma_Values:', Gamma_Values)
# mu_gamma = []
#
#
#
# # --- Options
# RUN = True
# # RUN = False
# # make_Plot = False
# make_Plot = True # vll besser : Plot_muGamma
#
# if RUN:
# for gamma in Gamma_Values:
# print("Run Cell-Problem for Gamma = ", gamma)
# # print('gamma='+str(gamma))
# with open(InputFilePath, 'r') as file:
# filedata = file.read()
# filedata = re.sub('(?m)^gamma=.*','gamma='+str(gamma),filedata)
# f = open(InputFilePath,'w')
# f.write(filedata)
# f.close()
# # --- Run Cell-Problem
# t = time.time()
# # subprocess.run(['./build-cmake/src/Cell-Problem', './inputs/cellsolver.parset'],
# # capture_output=True, text=True)
# # --- Run Cell-Problem_muGama -> faster
# # subprocess.run(['./build-cmake/src/Cell-Problem_muGamma', './inputs/cellsolver.parset'],
# # capture_output=True, text=True)
# # --- Run Compute_muGamma (2D Problem much much faster)
# subprocess.run(['./build-cmake/src/Compute_MuGamma', './inputs/computeMuGamma.parset'],
# capture_output=True, text=True)
# print('elapsed time:', time.time() - t)
#
# #Extract mu_gamma from Output-File TODO: GENERALIZED THIS FOR QUANTITIES OF INTEREST
# with open(OutputFilePath, 'r') as file:
# output = file.read()
# tmp = re.search(r'(?m)^mu_gamma=.*',output).group() # Not necessary for Intention of Program t output Minimizer etc.....
# s = re.findall(r"[-+]?\d*\.\d+|\d+", tmp)
# mu_gammaValue = float(s[0])
# print("mu_gamma:", mu_gammaValue)
# mu_gamma.append(mu_gammaValue)
# # ------------end of for-loop -----------------
# print("(Output) Values of mu_gamma: ", mu_gamma)
# # ----------------end of if-statement -------------
#
# # mu_gamma=[2.06099, 1.90567, 1.905]
# # mu_gamma=[2.08306, 1.905, 1.90482, 1.90479, 1.90478, 1.90477]
#
# ##Gamma_Values = np.linspace(0.01, 20, num=12) :
# #mu_gamma= [2.08306, 1.91108, 1.90648, 1.90554, 1.90521, 1.90505, 1.90496, 1.90491, 1.90487, 1.90485, 1.90483, 1.90482]
#
# ##Gamma_Values = np.linspace(0.01, 2.5, num=12)
# # mu_gamma=[2.08306, 2.01137, 1.96113, 1.93772, 1.92592, 1.91937, 1.91541, 1.91286, 1.91112, 1.90988, 1.90897, 1.90828]
#
# Gamma_Values = np.linspace(0.01, 2.5, num=6)
# mu_gamma=[2.08306, 1.95497, 1.92287, 1.91375, 1.9101, 1.90828]
#
#
#
# # Make Plot
# if make_Plot:
# plt.figure()
# plt.title(r'$\mu_\gamma(\gamma)$-Plot')
# plt.plot(Gamma_Values, mu_gamma)
# plt.scatter(Gamma_Values, mu_gamma)
# # plt.axis([0, 6, 0, 20])
# plt.axhline(y = 1.90476, color = 'b', linestyle = ':', label='$q_1$')
# plt.axhline(y = 2.08333, color = 'r', linestyle = 'dashed', label='$q_2$')
# plt.xlabel("$\gamma$")
# plt.ylabel("$\mu_\gamma$")
# plt.legend()
# plt.show()
#
# # ------------------------------------------------------------------------
print('---- Input parameters: -----')
mu1 = 10.0
rho1 = 1.0
......@@ -298,23 +51,17 @@ print('----------------------------')
#
# print('RunCellProblem...')
# RunCellProblem(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
# TEST read Matrix file
# Test = loadmat(path + '/outputs/QMatrix.mat')
# TEST:
print('Compare_Classification...')
Compare_Classification(alpha,beta,theta,gamma,mu1,rho1,InputFilePath)
# ------------- RUN Matlab symbolic minimization program
# ------------- RUN Matlab symbolic minimization program 'symMinimization'
eng = matlab.engine.start_matlab()
# s = eng.genpath(path + '/Matlab-Programs')
s = eng.genpath(path)
......@@ -323,11 +70,12 @@ eng.addpath(s, nargout=0)
eng.cd('Matlab-Programs', nargout=0) #switch to Matlab-Programs folder
# print('current Matlab folder:', eng.pwd(nargout=1))
Inp = False
Inp_T = True
print('Run symbolic Minimization...')
G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp, nargout=4) #Name of program:symMinimization
#Arguments: symMinization(print_Input,print_statPoint,print_Output,make_FunctionPlot, InputPath)
G, angle, type, kappa = eng.symMinimization(Inp_T,Inp,Inp,Inp, nargout=4) #Name of program:symMinimization
# G, angle, type, kappa = eng.symMinimization(Inp,Inp,Inp,Inp,path + "/outputs", nargout=4) #Optional: add Path
G = np.asarray(G) #cast Matlab Outout to numpy array
# --- print Output ---
print('Minimizer G:')
print(G)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment