From 1c4e0f8de6faa2457f79e80645fde61a851b2db4 Mon Sep 17 00:00:00 2001
From: Kim Liegeois <kimliegeois@ymail.com>
Date: Mon, 10 Apr 2017 08:47:02 +0200
Subject: [PATCH] ICOMP second numerical SA
---
.vscode/settings.json | 13 +
Metafor/config/bord01_numericalSA2.py | 282 +++++++
Metafor/models/bord01/createNumericalSA.m | 49 +-
.../bord01/numericalSATrainingPoints2.ascii | 500 +++++++++++++
.../models/bord01/numericalSAWeights2.ascii | 500 +++++++++++++
Metafor/models/bord01/tombeBord.py | 68 +-
Metafor/models/bord01/tombeBord_old.py | 247 +++++++
mrstlnos/matlab/Belos/test_belos8.m | 25 +
mrstlnos/matlab/Multigrid/read_agregate.m | 20 +
.../matlab/Multigrid/testOneDProlongation.m | 86 ++-
.../matlab/Multigrid/testOneDRestriction.m | 82 ++-
mrstlnos/src/IterativeSolver.cpp | 12 +-
mrstlnos/src/Preconditioner.h | 7 +-
mrstlnos/src/wExample2.cpp | 687 +++++++++++++++---
mrstlnos/src/wExample2.h | 10 +-
mrstlnos/tests/bolt_iterative.py | 193 +++++
mrstlnos/tests/cube.geo | 2 +-
mrstlnos/tests/cube_new_iterative.py | 45 +-
mrstlnos/tests/example2.py | 32 +
19 files changed, 2674 insertions(+), 186 deletions(-)
create mode 100644 .vscode/settings.json
create mode 100644 Metafor/config/bord01_numericalSA2.py
create mode 100644 Metafor/models/bord01/numericalSATrainingPoints2.ascii
create mode 100644 Metafor/models/bord01/numericalSAWeights2.ascii
create mode 100644 Metafor/models/bord01/tombeBord_old.py
create mode 100644 mrstlnos/matlab/Belos/test_belos8.m
create mode 100644 mrstlnos/matlab/Multigrid/read_agregate.m
create mode 100644 mrstlnos/tests/bolt_iterative.py
create mode 100644 mrstlnos/tests/example2.py
diff --git a/.vscode/settings.json b/.vscode/settings.json
new file mode 100644
index 00000000..c390280e
--- /dev/null
+++ b/.vscode/settings.json
@@ -0,0 +1,13 @@
+{
+ "files.associations": {
+ "exception": "cpp",
+ "new": "cpp",
+ "cstdint": "cpp",
+ "tuple": "cpp",
+ "type_traits": "cpp",
+ "typeinfo": "cpp",
+ "initializer_list": "cpp",
+ "ostream": "cpp",
+ "cmath": "cpp"
+ }
+}
\ No newline at end of file
diff --git a/Metafor/config/bord01_numericalSA2.py b/Metafor/config/bord01_numericalSA2.py
new file mode 100644
index 00000000..4b6ce77a
--- /dev/null
+++ b/Metafor/config/bord01_numericalSA2.py
@@ -0,0 +1,282 @@
+#!/usr/bin/env python
+# -*- coding: latin-1; -*-
+
+import numpy as np
+import os
+import shutil
+
+import Metafor.Mparams.TaskManager as tm
+import Metafor.Mparams.mparams as mp
+
+
+def main():
+ if tm.rank == 0:
+ #--------------------------------------------------------
+ # Model parameters
+ #--------------------------------------------------------
+ model_parameters = {}
+ isUnix = lambda: os.name == 'posix'
+
+ if isUnix():
+ model_parameters['Metafor_exe'] = '/home/kliegeois/Codes/Metafor/oo_metaB/bin/Metafor'
+ model_parameters['Metafor_model'] = '/home/kliegeois/Codes/waves/Metafor/models/bord01'
+ else:
+ model_parameters['Metafor_exe'] = 'C:/Users/kimli_000/dev/oo_metaB/bin/Release/Metafor'
+ model_parameters['Metafor_model'] = 'C:/Users/kimli_000/dev/waves/Metafor/models/bord01'
+ #--------------------------------------------------------
+ # Surrogate parameters
+ #--------------------------------------------------------
+ surrogate_parameters = {}
+ surrogate_parameters['display'] = False
+ if not isUnix():
+ surrogate_parameters['display'] = False
+ surrogate_parameters['normalisation'] = False
+ #--------------------------------------------------------
+ surrogate_parameters['fname'] = 'mparams.log'
+ d = 6
+ surrogate_parameters['d'] = d
+ surrogate_parameters['P_name'] = np.array([ 'NRTol', 'peno', 'peta', 'PEAS', 'nx', 'ny'])
+
+ M = np.zeros((d,), dtype =np.int)
+ for j in range(0,d):
+ M[j] = 5
+
+ surrogate_parameters['M'] = M
+
+ surrogate_parameters['Dom_min'] = np.array([ 1400., 350.])
+ surrogate_parameters['Dom_max'] = np.array([ 1600., 450.])
+
+ N = np.zeros((d,), dtype =np.int)
+ for j in range(0,d):
+ N[j] = 2
+
+ surrogate_parameters['N'] = N
+ #--------------------------------------------------------
+ #--------------------------------------------------------
+ surrogate_parameters['function_type'] = 'monomial_1norm'
+ surrogate_parameters['training_points'] = 'precomputed'
+ surrogate_parameters['precomputed_training_points'] = '/home/kliegeois/Codes/waves/Metafor/models/bord01/numericalSATrainingPoints2.ascii'
+ surrogate_parameters['precomputed_weights'] = '/home/kliegeois/Codes/waves/Metafor/models/bord01/numericalSAWeights2.ascii'
+ surrogate_parameters['solver'] = 'QR'
+ surrogate_parameters['means'] = np.array([1495.,396.])
+ surrogate_parameters['variances'] = np.array([1390.,660.])
+ rho = -0.233
+ surrogate_parameters['correlations'] = np.array([[1.,rho],[rho,1.]])
+ surrogate_parameters['f_mapping_method'] = 'Cholesky'
+ surrogate_parameters['rho'] = rho
+ surrogate_parameters['nstddeviation'] = 2.6281407
+ #--------------------------------------------------------
+ surrogate_parameters['m1'] = 40
+ surrogate_parameters['m2'] = 40
+ #--------------------------------------------------------
+ # Images parameters
+ #--------------------------------------------------------
+ images_parameters = {}
+ images_parameters['xlabel'] = '$H\;[$MPa$]$'
+ images_parameters['ylabel'] = '$S\;[$MPa$]$'
+ images_parameters['zlabel'] = 'Angle $[$rad$]$'
+ images_parameters['z_min'] = 0.04
+ images_parameters['z_max'] = 0.07
+
+ if mp.AvailableComputedResults():
+ if tm.rank == 0:
+ print "the master/sol.json file will be loaded"
+ print "the master/Z.json file will be loaded"
+ print "the master/W_sr.json file will be loaded"
+ sol, Z, W_sr = mp.LoadComputedResults()
+
+
+ from Metafor.Mparams.sequential import Z_evaluate
+ from Metafor.Mparams.sequential import scatterplot
+ from Metafor.Mparams.surrogate import surrogate_creation
+
+ surrogate_parameters['nstddeviation'] = 1
+ Z1, W1_sr, Xi, Xi2 = Z_evaluate(model_parameters, surrogate_parameters)
+
+ i_max = 1
+ cond = np.zeros((i_max,1))
+ ZWZrank = np.zeros((i_max,1))
+ ZWZdet = np.zeros((i_max,1))
+ residual = np.zeros((i_max,1))
+ max_residual = np.zeros((i_max,1))
+ min_residual = np.zeros((i_max,1))
+ w_residual = np.zeros((i_max,1))
+ max_w_residual = np.zeros((i_max,1))
+ min_w_residual = np.zeros((i_max,1))
+
+ np.set_printoptions(linewidth = 210)
+
+ for j in range(0,d):
+ N[j] = 9
+
+ surrogate_parameters['N'] = N
+
+
+ nstddeviation = np.linspace(0.1,10.,i_max)
+ for i in range(0,i_max):
+ surrogate_parameters['nstddeviation'] = nstddeviation[i]
+ Z, W_sr, Xi, Xi2 = Z_evaluate(model_parameters, surrogate_parameters)
+
+
+ #W_sr = np.eye(W_sr.shape[0])
+
+
+
+ W_srZ = np.dot(W_sr,Z)
+ cond[i] = np.linalg.cond(W_srZ)
+
+ '''
+ Image = scatterplot(sol, surrogate_parameters, images_parameters)
+ s_hat = surrogate_creation(sol, Z, W_sr, surrogate_parameters, Image, images_parameters)
+
+ ZWZ = np.dot(np.transpose(W_srZ),W_srZ)
+
+ ZWZrank[i] = np.linalg.matrix_rank(ZWZ)
+ ZWZdet[i] = np.linalg.det(ZWZ)
+ res = np.dot(Z,s_hat)-sol
+ residual[i] = np.linalg.norm(res)
+ max_residual[i] = np.abs(res).max()
+ min_residual[i] = np.abs(res).min()
+
+ w_res = np.dot(W,res)
+ w_residual[i] = np.linalg.norm(w_res)
+ max_w_residual[i] = np.abs(w_res).max()
+ min_w_residual[i] = np.abs(w_res).min()
+ '''
+ print i
+
+ #import matplotlib.pyplot as plt
+ #plt.semilogy(nstddeviation, cond)
+ #plt.show()
+ print 'Condition number of W_srZ'
+ print cond.T
+ print 'nstddeviation'
+ print nstddeviation.T
+
+ #'''
+ surrogate_parameters['nstddeviation'] = 2.6281407
+ for j in range(0,d):
+ N[j] = 9
+
+ surrogate_parameters['N'] = N
+ surrogate_parameters['display'] = True
+ Z, W_sr, Xi, Xi2 = Z_evaluate(model_parameters, surrogate_parameters)
+
+ Image = scatterplot(sol, surrogate_parameters, images_parameters)
+ s_hat_w = surrogate_creation(sol, Z, W_sr, surrogate_parameters, Image, images_parameters)
+ print s_hat_w.T
+
+ W_sr = np.eye(W_sr.shape[0])
+
+ Image = scatterplot(sol, surrogate_parameters, images_parameters)
+ s_hat_wout = surrogate_creation(sol, Z, W_sr, surrogate_parameters, Image, images_parameters)
+ print s_hat_wout.T
+
+ ## ---------------------------
+ from params.src.monomials_1norm import phi
+ from params.src.monomials_1norm import compute_n_max
+ n_max = compute_n_max(N,d)
+ import params.arnst_ponthot.src.gamma as tp
+
+ means = surrogate_parameters['means']
+ variances = surrogate_parameters['variances']
+ correlations = surrogate_parameters['correlations']
+
+ n_MC = 200
+
+
+ values_w = np.zeros((n_MC,1))
+ values_wout = np.zeros((n_MC,1))
+ xi_s = np.zeros((2,n_MC))
+ for i in range(0,n_MC):
+ if np.mod(i,1000) == 0:
+ print i
+ xi_in = np.random.normal(0., 1., 2)
+
+
+
+ xi_in[0] = 5.273790423018157
+
+
+ xi = tp.f_mapping(xi_in,2,means,variances,correlations,surrogate_parameters['f_mapping_method'])
+ xi_s[0,i] = (xi[0]-means[0])/(surrogate_parameters['nstddeviation']*np.sqrt(variances[0]))
+ xi_s[1,i] = (xi[1]-means[1])/(surrogate_parameters['nstddeviation']*np.sqrt(variances[1]))
+ phi_xi = phi(xi_s[:,i],N,n_max,d)
+ values_w[i] = np.dot(np.transpose(phi_xi),s_hat_w)
+ values_wout[i] = np.dot(np.transpose(phi_xi),s_hat_wout)
+ '''
+ mean = 0.
+ variance = 0.
+ exceedence = 0.
+ for i in range(0,n_MC):
+ mean = mean + values[i]
+
+ if values[i] >= 0.065:
+ exceedence = exceedence + 1.
+
+ mean = mean/n_MC
+
+ xi_s = np.zeros((2,1))
+ for i in range(0,n_MC):
+ variance = variance + (mean - values[i])**2
+
+ variance = variance/(n_MC-1.)
+
+ print 'Number of MC points'
+ print n_MC
+ print 'Mean'
+ print mean
+ print 'Variance'
+ print variance
+ print 'exceedence'
+ print exceedence
+ '''
+ '''
+ import codecs, json
+ tmp = values_w.tolist()
+ file_path = "values_w.json"
+ json.dump(tmp, codecs.open(file_path, 'w', encoding='utf-8'), separators=(',', ':'), sort_keys=True, indent=4)
+
+ tmp = values_wout.tolist()
+ file_path = "values_wout.json"
+ json.dump(tmp, codecs.open(file_path, 'w', encoding='utf-8'), separators=(',', ':'), sort_keys=True, indent=4)
+
+ tmp = xi_s.tolist()
+ file_path = "xi_s.json"
+ json.dump(tmp, codecs.open(file_path, 'w', encoding='utf-8'), separators=(',', ':'), sort_keys=True, indent=4)
+ '''
+ #import matplotlib.pyplot as plt
+ #count, bins, ignored = plt.hist(values, 100, normed=True)
+ #plt.show()
+ #'''
+ else:
+ if tm.rank == 0:
+ ##
+
+ sol, Z, W_sr, s_hat = mp.master(model_parameters,surrogate_parameters,images_parameters)
+
+ ##
+
+ import codecs, json
+
+ tmp = sol.tolist()
+ file_path = "sol.json"
+ json.dump(tmp, codecs.open(file_path, 'w', encoding='utf-8'), separators=(',', ':'), sort_keys=True, indent=4)
+
+ tmp = Z.tolist()
+ file_path = "Z.json"
+ json.dump(tmp, codecs.open(file_path, 'w', encoding='utf-8'), separators=(',', ':'), sort_keys=True, indent=4)
+
+ tmp = W_sr.tolist()
+ file_path = "W_sr.json"
+ json.dump(tmp, codecs.open(file_path, 'w', encoding='utf-8'), separators=(',', ':'), sort_keys=True, indent=4)
+
+ tmp = s_hat.tolist()
+ file_path = "s_hat.json"
+ json.dump(tmp, codecs.open(file_path, 'w', encoding='utf-8'), separators=(',', ':'), sort_keys=True, indent=4)
+
+ else:
+ mp.worker()
+
+if __name__ == "__main__":
+ main()
diff --git a/Metafor/models/bord01/createNumericalSA.m b/Metafor/models/bord01/createNumericalSA.m
index 135fcefe..4e529936 100644
--- a/Metafor/models/bord01/createNumericalSA.m
+++ b/Metafor/models/bord01/createNumericalSA.m
@@ -4,6 +4,7 @@ clc
format long
+%% old
NRTolExpo = -5:-0.5:-8;
NRTol = 10.^(NRTolExpo);
penoExpo = 5: 0.5: 8;
@@ -11,6 +12,8 @@ peno = 10.^(penoExpo);
nx = 20:5:60;
ny = 4:2:16;
+
+
i = 1;
trainingPoints = zeros(length(NRTol)*length(nx)*length(ny)*length(peno),4);
weights = ones(length(NRTol)*length(nx)*length(ny)*length(peno),1);
@@ -33,4 +36,48 @@ for i = 1:length(NRTol)*length(nx)*length(ny)*length(peno)
end
fclose(fid)
%save numericalSATrainingPoints.ascii trainingPoints -ASCII -DOUBLE;
-save numericalSAWeights.ascii weights -ASCII -DOUBLE;
\ No newline at end of file
+save numericalSAWeights.ascii weights -ASCII -DOUBLE;
+
+clear all
+close all
+clc
+
+%% new
+NRTolExpo = linspace(-4,-5,5);
+NRTol = 10.^(NRTolExpo);
+
+penoExpo = linspace(5,6,5);
+peno = 10.^(penoExpo);
+peta = peno/10;
+
+PEASExpo = linspace(-6,-7,5);
+PEAS = 10.^(NRTolExpo);
+
+ny = 6:2:12;
+nx = ny*10;
+
+i = 1;
+trainingPoints = zeros(length(NRTol)*length(peno)*length(PEAS)*length(ny),6);
+weights = ones(length(NRTol)*length(peno)*length(PEAS)*length(ny),1);
+for i_NR = 1:length(NRTol)
+for i_peno= 1:length(peno)
+for i_PEAS = 1:length(PEAS)
+for i_ny = 1:length(ny)
+ trainingPoints(i,:) = [NRTol(i_NR), peno(i_peno), peta(i_peno), PEAS(i_PEAS), nx(i_ny), ny(i_ny)];
+ i = i+1;
+end
+end
+end
+end
+ fid = fopen('numericalSATrainingPoints2.ascii', 'wt');
+ for i = 1:length(NRTol)*length(peno)*length(PEAS)*length(ny)
+ fprintf(fid,'%e ',trainingPoints(i,1));
+ fprintf(fid,'%e ',trainingPoints(i,2));
+ fprintf(fid,'%e ',trainingPoints(i,3));
+ fprintf(fid,'%e ',trainingPoints(i,4));
+ fprintf(fid,'%d ',trainingPoints(i,5));
+ fprintf(fid,'%d \n',trainingPoints(i,6));
+ end
+ fclose(fid)
+%save numericalSATrainingPoints2.ascii trainingPoints -ASCII -DOUBLE;
+save numericalSAWeights2.ascii weights -ASCII -DOUBLE;
\ No newline at end of file
diff --git a/Metafor/models/bord01/numericalSATrainingPoints2.ascii b/Metafor/models/bord01/numericalSATrainingPoints2.ascii
new file mode 100644
index 00000000..c92aa1e2
--- /dev/null
+++ b/Metafor/models/bord01/numericalSATrainingPoints2.ascii
@@ -0,0 +1,500 @@
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-04 60 6
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-04 80 8
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-04 100 10
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-04 120 12
+1.000000e-04 1.000000e+05 1.000000e+04 5.623413e-05 60 6
+1.000000e-04 1.000000e+05 1.000000e+04 5.623413e-05 80 8
+1.000000e-04 1.000000e+05 1.000000e+04 5.623413e-05 100 10
+1.000000e-04 1.000000e+05 1.000000e+04 5.623413e-05 120 12
+1.000000e-04 1.000000e+05 1.000000e+04 3.162278e-05 60 6
+1.000000e-04 1.000000e+05 1.000000e+04 3.162278e-05 80 8
+1.000000e-04 1.000000e+05 1.000000e+04 3.162278e-05 100 10
+1.000000e-04 1.000000e+05 1.000000e+04 3.162278e-05 120 12
+1.000000e-04 1.000000e+05 1.000000e+04 1.778279e-05 60 6
+1.000000e-04 1.000000e+05 1.000000e+04 1.778279e-05 80 8
+1.000000e-04 1.000000e+05 1.000000e+04 1.778279e-05 100 10
+1.000000e-04 1.000000e+05 1.000000e+04 1.778279e-05 120 12
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-05 60 6
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-05 80 8
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-05 100 10
+1.000000e-04 1.000000e+05 1.000000e+04 1.000000e-05 120 12
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-04 60 6
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-04 80 8
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-04 100 10
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-04 120 12
+1.000000e-04 1.778279e+05 1.778279e+04 5.623413e-05 60 6
+1.000000e-04 1.778279e+05 1.778279e+04 5.623413e-05 80 8
+1.000000e-04 1.778279e+05 1.778279e+04 5.623413e-05 100 10
+1.000000e-04 1.778279e+05 1.778279e+04 5.623413e-05 120 12
+1.000000e-04 1.778279e+05 1.778279e+04 3.162278e-05 60 6
+1.000000e-04 1.778279e+05 1.778279e+04 3.162278e-05 80 8
+1.000000e-04 1.778279e+05 1.778279e+04 3.162278e-05 100 10
+1.000000e-04 1.778279e+05 1.778279e+04 3.162278e-05 120 12
+1.000000e-04 1.778279e+05 1.778279e+04 1.778279e-05 60 6
+1.000000e-04 1.778279e+05 1.778279e+04 1.778279e-05 80 8
+1.000000e-04 1.778279e+05 1.778279e+04 1.778279e-05 100 10
+1.000000e-04 1.778279e+05 1.778279e+04 1.778279e-05 120 12
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-05 60 6
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-05 80 8
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-05 100 10
+1.000000e-04 1.778279e+05 1.778279e+04 1.000000e-05 120 12
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-04 60 6
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-04 80 8
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-04 100 10
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-04 120 12
+1.000000e-04 3.162278e+05 3.162278e+04 5.623413e-05 60 6
+1.000000e-04 3.162278e+05 3.162278e+04 5.623413e-05 80 8
+1.000000e-04 3.162278e+05 3.162278e+04 5.623413e-05 100 10
+1.000000e-04 3.162278e+05 3.162278e+04 5.623413e-05 120 12
+1.000000e-04 3.162278e+05 3.162278e+04 3.162278e-05 60 6
+1.000000e-04 3.162278e+05 3.162278e+04 3.162278e-05 80 8
+1.000000e-04 3.162278e+05 3.162278e+04 3.162278e-05 100 10
+1.000000e-04 3.162278e+05 3.162278e+04 3.162278e-05 120 12
+1.000000e-04 3.162278e+05 3.162278e+04 1.778279e-05 60 6
+1.000000e-04 3.162278e+05 3.162278e+04 1.778279e-05 80 8
+1.000000e-04 3.162278e+05 3.162278e+04 1.778279e-05 100 10
+1.000000e-04 3.162278e+05 3.162278e+04 1.778279e-05 120 12
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-05 60 6
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-05 80 8
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-05 100 10
+1.000000e-04 3.162278e+05 3.162278e+04 1.000000e-05 120 12
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-04 60 6
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-04 80 8
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-04 100 10
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-04 120 12
+1.000000e-04 5.623413e+05 5.623413e+04 5.623413e-05 60 6
+1.000000e-04 5.623413e+05 5.623413e+04 5.623413e-05 80 8
+1.000000e-04 5.623413e+05 5.623413e+04 5.623413e-05 100 10
+1.000000e-04 5.623413e+05 5.623413e+04 5.623413e-05 120 12
+1.000000e-04 5.623413e+05 5.623413e+04 3.162278e-05 60 6
+1.000000e-04 5.623413e+05 5.623413e+04 3.162278e-05 80 8
+1.000000e-04 5.623413e+05 5.623413e+04 3.162278e-05 100 10
+1.000000e-04 5.623413e+05 5.623413e+04 3.162278e-05 120 12
+1.000000e-04 5.623413e+05 5.623413e+04 1.778279e-05 60 6
+1.000000e-04 5.623413e+05 5.623413e+04 1.778279e-05 80 8
+1.000000e-04 5.623413e+05 5.623413e+04 1.778279e-05 100 10
+1.000000e-04 5.623413e+05 5.623413e+04 1.778279e-05 120 12
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-05 60 6
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-05 80 8
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-05 100 10
+1.000000e-04 5.623413e+05 5.623413e+04 1.000000e-05 120 12
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-04 60 6
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-04 80 8
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-04 100 10
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-04 120 12
+1.000000e-04 1.000000e+06 1.000000e+05 5.623413e-05 60 6
+1.000000e-04 1.000000e+06 1.000000e+05 5.623413e-05 80 8
+1.000000e-04 1.000000e+06 1.000000e+05 5.623413e-05 100 10
+1.000000e-04 1.000000e+06 1.000000e+05 5.623413e-05 120 12
+1.000000e-04 1.000000e+06 1.000000e+05 3.162278e-05 60 6
+1.000000e-04 1.000000e+06 1.000000e+05 3.162278e-05 80 8
+1.000000e-04 1.000000e+06 1.000000e+05 3.162278e-05 100 10
+1.000000e-04 1.000000e+06 1.000000e+05 3.162278e-05 120 12
+1.000000e-04 1.000000e+06 1.000000e+05 1.778279e-05 60 6
+1.000000e-04 1.000000e+06 1.000000e+05 1.778279e-05 80 8
+1.000000e-04 1.000000e+06 1.000000e+05 1.778279e-05 100 10
+1.000000e-04 1.000000e+06 1.000000e+05 1.778279e-05 120 12
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-05 60 6
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-05 80 8
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-05 100 10
+1.000000e-04 1.000000e+06 1.000000e+05 1.000000e-05 120 12
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-04 60 6
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-04 80 8
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-04 100 10
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-04 120 12
+5.623413e-05 1.000000e+05 1.000000e+04 5.623413e-05 60 6
+5.623413e-05 1.000000e+05 1.000000e+04 5.623413e-05 80 8
+5.623413e-05 1.000000e+05 1.000000e+04 5.623413e-05 100 10
+5.623413e-05 1.000000e+05 1.000000e+04 5.623413e-05 120 12
+5.623413e-05 1.000000e+05 1.000000e+04 3.162278e-05 60 6
+5.623413e-05 1.000000e+05 1.000000e+04 3.162278e-05 80 8
+5.623413e-05 1.000000e+05 1.000000e+04 3.162278e-05 100 10
+5.623413e-05 1.000000e+05 1.000000e+04 3.162278e-05 120 12
+5.623413e-05 1.000000e+05 1.000000e+04 1.778279e-05 60 6
+5.623413e-05 1.000000e+05 1.000000e+04 1.778279e-05 80 8
+5.623413e-05 1.000000e+05 1.000000e+04 1.778279e-05 100 10
+5.623413e-05 1.000000e+05 1.000000e+04 1.778279e-05 120 12
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-05 60 6
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-05 80 8
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-05 100 10
+5.623413e-05 1.000000e+05 1.000000e+04 1.000000e-05 120 12
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-04 60 6
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-04 80 8
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-04 100 10
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-04 120 12
+5.623413e-05 1.778279e+05 1.778279e+04 5.623413e-05 60 6
+5.623413e-05 1.778279e+05 1.778279e+04 5.623413e-05 80 8
+5.623413e-05 1.778279e+05 1.778279e+04 5.623413e-05 100 10
+5.623413e-05 1.778279e+05 1.778279e+04 5.623413e-05 120 12
+5.623413e-05 1.778279e+05 1.778279e+04 3.162278e-05 60 6
+5.623413e-05 1.778279e+05 1.778279e+04 3.162278e-05 80 8
+5.623413e-05 1.778279e+05 1.778279e+04 3.162278e-05 100 10
+5.623413e-05 1.778279e+05 1.778279e+04 3.162278e-05 120 12
+5.623413e-05 1.778279e+05 1.778279e+04 1.778279e-05 60 6
+5.623413e-05 1.778279e+05 1.778279e+04 1.778279e-05 80 8
+5.623413e-05 1.778279e+05 1.778279e+04 1.778279e-05 100 10
+5.623413e-05 1.778279e+05 1.778279e+04 1.778279e-05 120 12
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-05 60 6
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-05 80 8
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-05 100 10
+5.623413e-05 1.778279e+05 1.778279e+04 1.000000e-05 120 12
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-04 60 6
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-04 80 8
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-04 100 10
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-04 120 12
+5.623413e-05 3.162278e+05 3.162278e+04 5.623413e-05 60 6
+5.623413e-05 3.162278e+05 3.162278e+04 5.623413e-05 80 8
+5.623413e-05 3.162278e+05 3.162278e+04 5.623413e-05 100 10
+5.623413e-05 3.162278e+05 3.162278e+04 5.623413e-05 120 12
+5.623413e-05 3.162278e+05 3.162278e+04 3.162278e-05 60 6
+5.623413e-05 3.162278e+05 3.162278e+04 3.162278e-05 80 8
+5.623413e-05 3.162278e+05 3.162278e+04 3.162278e-05 100 10
+5.623413e-05 3.162278e+05 3.162278e+04 3.162278e-05 120 12
+5.623413e-05 3.162278e+05 3.162278e+04 1.778279e-05 60 6
+5.623413e-05 3.162278e+05 3.162278e+04 1.778279e-05 80 8
+5.623413e-05 3.162278e+05 3.162278e+04 1.778279e-05 100 10
+5.623413e-05 3.162278e+05 3.162278e+04 1.778279e-05 120 12
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-05 60 6
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-05 80 8
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-05 100 10
+5.623413e-05 3.162278e+05 3.162278e+04 1.000000e-05 120 12
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-04 60 6
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-04 80 8
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-04 100 10
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-04 120 12
+5.623413e-05 5.623413e+05 5.623413e+04 5.623413e-05 60 6
+5.623413e-05 5.623413e+05 5.623413e+04 5.623413e-05 80 8
+5.623413e-05 5.623413e+05 5.623413e+04 5.623413e-05 100 10
+5.623413e-05 5.623413e+05 5.623413e+04 5.623413e-05 120 12
+5.623413e-05 5.623413e+05 5.623413e+04 3.162278e-05 60 6
+5.623413e-05 5.623413e+05 5.623413e+04 3.162278e-05 80 8
+5.623413e-05 5.623413e+05 5.623413e+04 3.162278e-05 100 10
+5.623413e-05 5.623413e+05 5.623413e+04 3.162278e-05 120 12
+5.623413e-05 5.623413e+05 5.623413e+04 1.778279e-05 60 6
+5.623413e-05 5.623413e+05 5.623413e+04 1.778279e-05 80 8
+5.623413e-05 5.623413e+05 5.623413e+04 1.778279e-05 100 10
+5.623413e-05 5.623413e+05 5.623413e+04 1.778279e-05 120 12
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-05 60 6
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-05 80 8
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-05 100 10
+5.623413e-05 5.623413e+05 5.623413e+04 1.000000e-05 120 12
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-04 60 6
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-04 80 8
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-04 100 10
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-04 120 12
+5.623413e-05 1.000000e+06 1.000000e+05 5.623413e-05 60 6
+5.623413e-05 1.000000e+06 1.000000e+05 5.623413e-05 80 8
+5.623413e-05 1.000000e+06 1.000000e+05 5.623413e-05 100 10
+5.623413e-05 1.000000e+06 1.000000e+05 5.623413e-05 120 12
+5.623413e-05 1.000000e+06 1.000000e+05 3.162278e-05 60 6
+5.623413e-05 1.000000e+06 1.000000e+05 3.162278e-05 80 8
+5.623413e-05 1.000000e+06 1.000000e+05 3.162278e-05 100 10
+5.623413e-05 1.000000e+06 1.000000e+05 3.162278e-05 120 12
+5.623413e-05 1.000000e+06 1.000000e+05 1.778279e-05 60 6
+5.623413e-05 1.000000e+06 1.000000e+05 1.778279e-05 80 8
+5.623413e-05 1.000000e+06 1.000000e+05 1.778279e-05 100 10
+5.623413e-05 1.000000e+06 1.000000e+05 1.778279e-05 120 12
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-05 60 6
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-05 80 8
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-05 100 10
+5.623413e-05 1.000000e+06 1.000000e+05 1.000000e-05 120 12
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-04 60 6
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-04 80 8
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-04 100 10
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-04 120 12
+3.162278e-05 1.000000e+05 1.000000e+04 5.623413e-05 60 6
+3.162278e-05 1.000000e+05 1.000000e+04 5.623413e-05 80 8
+3.162278e-05 1.000000e+05 1.000000e+04 5.623413e-05 100 10
+3.162278e-05 1.000000e+05 1.000000e+04 5.623413e-05 120 12
+3.162278e-05 1.000000e+05 1.000000e+04 3.162278e-05 60 6
+3.162278e-05 1.000000e+05 1.000000e+04 3.162278e-05 80 8
+3.162278e-05 1.000000e+05 1.000000e+04 3.162278e-05 100 10
+3.162278e-05 1.000000e+05 1.000000e+04 3.162278e-05 120 12
+3.162278e-05 1.000000e+05 1.000000e+04 1.778279e-05 60 6
+3.162278e-05 1.000000e+05 1.000000e+04 1.778279e-05 80 8
+3.162278e-05 1.000000e+05 1.000000e+04 1.778279e-05 100 10
+3.162278e-05 1.000000e+05 1.000000e+04 1.778279e-05 120 12
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-05 60 6
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-05 80 8
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-05 100 10
+3.162278e-05 1.000000e+05 1.000000e+04 1.000000e-05 120 12
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-04 60 6
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-04 80 8
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-04 100 10
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-04 120 12
+3.162278e-05 1.778279e+05 1.778279e+04 5.623413e-05 60 6
+3.162278e-05 1.778279e+05 1.778279e+04 5.623413e-05 80 8
+3.162278e-05 1.778279e+05 1.778279e+04 5.623413e-05 100 10
+3.162278e-05 1.778279e+05 1.778279e+04 5.623413e-05 120 12
+3.162278e-05 1.778279e+05 1.778279e+04 3.162278e-05 60 6
+3.162278e-05 1.778279e+05 1.778279e+04 3.162278e-05 80 8
+3.162278e-05 1.778279e+05 1.778279e+04 3.162278e-05 100 10
+3.162278e-05 1.778279e+05 1.778279e+04 3.162278e-05 120 12
+3.162278e-05 1.778279e+05 1.778279e+04 1.778279e-05 60 6
+3.162278e-05 1.778279e+05 1.778279e+04 1.778279e-05 80 8
+3.162278e-05 1.778279e+05 1.778279e+04 1.778279e-05 100 10
+3.162278e-05 1.778279e+05 1.778279e+04 1.778279e-05 120 12
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-05 60 6
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-05 80 8
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-05 100 10
+3.162278e-05 1.778279e+05 1.778279e+04 1.000000e-05 120 12
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-04 60 6
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-04 80 8
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-04 100 10
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-04 120 12
+3.162278e-05 3.162278e+05 3.162278e+04 5.623413e-05 60 6
+3.162278e-05 3.162278e+05 3.162278e+04 5.623413e-05 80 8
+3.162278e-05 3.162278e+05 3.162278e+04 5.623413e-05 100 10
+3.162278e-05 3.162278e+05 3.162278e+04 5.623413e-05 120 12
+3.162278e-05 3.162278e+05 3.162278e+04 3.162278e-05 60 6
+3.162278e-05 3.162278e+05 3.162278e+04 3.162278e-05 80 8
+3.162278e-05 3.162278e+05 3.162278e+04 3.162278e-05 100 10
+3.162278e-05 3.162278e+05 3.162278e+04 3.162278e-05 120 12
+3.162278e-05 3.162278e+05 3.162278e+04 1.778279e-05 60 6
+3.162278e-05 3.162278e+05 3.162278e+04 1.778279e-05 80 8
+3.162278e-05 3.162278e+05 3.162278e+04 1.778279e-05 100 10
+3.162278e-05 3.162278e+05 3.162278e+04 1.778279e-05 120 12
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-05 60 6
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-05 80 8
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-05 100 10
+3.162278e-05 3.162278e+05 3.162278e+04 1.000000e-05 120 12
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-04 60 6
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-04 80 8
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-04 100 10
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-04 120 12
+3.162278e-05 5.623413e+05 5.623413e+04 5.623413e-05 60 6
+3.162278e-05 5.623413e+05 5.623413e+04 5.623413e-05 80 8
+3.162278e-05 5.623413e+05 5.623413e+04 5.623413e-05 100 10
+3.162278e-05 5.623413e+05 5.623413e+04 5.623413e-05 120 12
+3.162278e-05 5.623413e+05 5.623413e+04 3.162278e-05 60 6
+3.162278e-05 5.623413e+05 5.623413e+04 3.162278e-05 80 8
+3.162278e-05 5.623413e+05 5.623413e+04 3.162278e-05 100 10
+3.162278e-05 5.623413e+05 5.623413e+04 3.162278e-05 120 12
+3.162278e-05 5.623413e+05 5.623413e+04 1.778279e-05 60 6
+3.162278e-05 5.623413e+05 5.623413e+04 1.778279e-05 80 8
+3.162278e-05 5.623413e+05 5.623413e+04 1.778279e-05 100 10
+3.162278e-05 5.623413e+05 5.623413e+04 1.778279e-05 120 12
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-05 60 6
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-05 80 8
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-05 100 10
+3.162278e-05 5.623413e+05 5.623413e+04 1.000000e-05 120 12
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-04 60 6
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-04 80 8
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-04 100 10
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-04 120 12
+3.162278e-05 1.000000e+06 1.000000e+05 5.623413e-05 60 6
+3.162278e-05 1.000000e+06 1.000000e+05 5.623413e-05 80 8
+3.162278e-05 1.000000e+06 1.000000e+05 5.623413e-05 100 10
+3.162278e-05 1.000000e+06 1.000000e+05 5.623413e-05 120 12
+3.162278e-05 1.000000e+06 1.000000e+05 3.162278e-05 60 6
+3.162278e-05 1.000000e+06 1.000000e+05 3.162278e-05 80 8
+3.162278e-05 1.000000e+06 1.000000e+05 3.162278e-05 100 10
+3.162278e-05 1.000000e+06 1.000000e+05 3.162278e-05 120 12
+3.162278e-05 1.000000e+06 1.000000e+05 1.778279e-05 60 6
+3.162278e-05 1.000000e+06 1.000000e+05 1.778279e-05 80 8
+3.162278e-05 1.000000e+06 1.000000e+05 1.778279e-05 100 10
+3.162278e-05 1.000000e+06 1.000000e+05 1.778279e-05 120 12
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-05 60 6
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-05 80 8
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-05 100 10
+3.162278e-05 1.000000e+06 1.000000e+05 1.000000e-05 120 12
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-04 60 6
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-04 80 8
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-04 100 10
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-04 120 12
+1.778279e-05 1.000000e+05 1.000000e+04 5.623413e-05 60 6
+1.778279e-05 1.000000e+05 1.000000e+04 5.623413e-05 80 8
+1.778279e-05 1.000000e+05 1.000000e+04 5.623413e-05 100 10
+1.778279e-05 1.000000e+05 1.000000e+04 5.623413e-05 120 12
+1.778279e-05 1.000000e+05 1.000000e+04 3.162278e-05 60 6
+1.778279e-05 1.000000e+05 1.000000e+04 3.162278e-05 80 8
+1.778279e-05 1.000000e+05 1.000000e+04 3.162278e-05 100 10
+1.778279e-05 1.000000e+05 1.000000e+04 3.162278e-05 120 12
+1.778279e-05 1.000000e+05 1.000000e+04 1.778279e-05 60 6
+1.778279e-05 1.000000e+05 1.000000e+04 1.778279e-05 80 8
+1.778279e-05 1.000000e+05 1.000000e+04 1.778279e-05 100 10
+1.778279e-05 1.000000e+05 1.000000e+04 1.778279e-05 120 12
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-05 60 6
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-05 80 8
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-05 100 10
+1.778279e-05 1.000000e+05 1.000000e+04 1.000000e-05 120 12
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-04 60 6
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-04 80 8
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-04 100 10
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-04 120 12
+1.778279e-05 1.778279e+05 1.778279e+04 5.623413e-05 60 6
+1.778279e-05 1.778279e+05 1.778279e+04 5.623413e-05 80 8
+1.778279e-05 1.778279e+05 1.778279e+04 5.623413e-05 100 10
+1.778279e-05 1.778279e+05 1.778279e+04 5.623413e-05 120 12
+1.778279e-05 1.778279e+05 1.778279e+04 3.162278e-05 60 6
+1.778279e-05 1.778279e+05 1.778279e+04 3.162278e-05 80 8
+1.778279e-05 1.778279e+05 1.778279e+04 3.162278e-05 100 10
+1.778279e-05 1.778279e+05 1.778279e+04 3.162278e-05 120 12
+1.778279e-05 1.778279e+05 1.778279e+04 1.778279e-05 60 6
+1.778279e-05 1.778279e+05 1.778279e+04 1.778279e-05 80 8
+1.778279e-05 1.778279e+05 1.778279e+04 1.778279e-05 100 10
+1.778279e-05 1.778279e+05 1.778279e+04 1.778279e-05 120 12
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-05 60 6
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-05 80 8
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-05 100 10
+1.778279e-05 1.778279e+05 1.778279e+04 1.000000e-05 120 12
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-04 60 6
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-04 80 8
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-04 100 10
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-04 120 12
+1.778279e-05 3.162278e+05 3.162278e+04 5.623413e-05 60 6
+1.778279e-05 3.162278e+05 3.162278e+04 5.623413e-05 80 8
+1.778279e-05 3.162278e+05 3.162278e+04 5.623413e-05 100 10
+1.778279e-05 3.162278e+05 3.162278e+04 5.623413e-05 120 12
+1.778279e-05 3.162278e+05 3.162278e+04 3.162278e-05 60 6
+1.778279e-05 3.162278e+05 3.162278e+04 3.162278e-05 80 8
+1.778279e-05 3.162278e+05 3.162278e+04 3.162278e-05 100 10
+1.778279e-05 3.162278e+05 3.162278e+04 3.162278e-05 120 12
+1.778279e-05 3.162278e+05 3.162278e+04 1.778279e-05 60 6
+1.778279e-05 3.162278e+05 3.162278e+04 1.778279e-05 80 8
+1.778279e-05 3.162278e+05 3.162278e+04 1.778279e-05 100 10
+1.778279e-05 3.162278e+05 3.162278e+04 1.778279e-05 120 12
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-05 60 6
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-05 80 8
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-05 100 10
+1.778279e-05 3.162278e+05 3.162278e+04 1.000000e-05 120 12
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-04 60 6
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-04 80 8
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-04 100 10
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-04 120 12
+1.778279e-05 5.623413e+05 5.623413e+04 5.623413e-05 60 6
+1.778279e-05 5.623413e+05 5.623413e+04 5.623413e-05 80 8
+1.778279e-05 5.623413e+05 5.623413e+04 5.623413e-05 100 10
+1.778279e-05 5.623413e+05 5.623413e+04 5.623413e-05 120 12
+1.778279e-05 5.623413e+05 5.623413e+04 3.162278e-05 60 6
+1.778279e-05 5.623413e+05 5.623413e+04 3.162278e-05 80 8
+1.778279e-05 5.623413e+05 5.623413e+04 3.162278e-05 100 10
+1.778279e-05 5.623413e+05 5.623413e+04 3.162278e-05 120 12
+1.778279e-05 5.623413e+05 5.623413e+04 1.778279e-05 60 6
+1.778279e-05 5.623413e+05 5.623413e+04 1.778279e-05 80 8
+1.778279e-05 5.623413e+05 5.623413e+04 1.778279e-05 100 10
+1.778279e-05 5.623413e+05 5.623413e+04 1.778279e-05 120 12
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-05 60 6
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-05 80 8
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-05 100 10
+1.778279e-05 5.623413e+05 5.623413e+04 1.000000e-05 120 12
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-04 60 6
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-04 80 8
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-04 100 10
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-04 120 12
+1.778279e-05 1.000000e+06 1.000000e+05 5.623413e-05 60 6
+1.778279e-05 1.000000e+06 1.000000e+05 5.623413e-05 80 8
+1.778279e-05 1.000000e+06 1.000000e+05 5.623413e-05 100 10
+1.778279e-05 1.000000e+06 1.000000e+05 5.623413e-05 120 12
+1.778279e-05 1.000000e+06 1.000000e+05 3.162278e-05 60 6
+1.778279e-05 1.000000e+06 1.000000e+05 3.162278e-05 80 8
+1.778279e-05 1.000000e+06 1.000000e+05 3.162278e-05 100 10
+1.778279e-05 1.000000e+06 1.000000e+05 3.162278e-05 120 12
+1.778279e-05 1.000000e+06 1.000000e+05 1.778279e-05 60 6
+1.778279e-05 1.000000e+06 1.000000e+05 1.778279e-05 80 8
+1.778279e-05 1.000000e+06 1.000000e+05 1.778279e-05 100 10
+1.778279e-05 1.000000e+06 1.000000e+05 1.778279e-05 120 12
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-05 60 6
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-05 80 8
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-05 100 10
+1.778279e-05 1.000000e+06 1.000000e+05 1.000000e-05 120 12
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-04 60 6
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-04 80 8
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-04 100 10
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-04 120 12
+1.000000e-05 1.000000e+05 1.000000e+04 5.623413e-05 60 6
+1.000000e-05 1.000000e+05 1.000000e+04 5.623413e-05 80 8
+1.000000e-05 1.000000e+05 1.000000e+04 5.623413e-05 100 10
+1.000000e-05 1.000000e+05 1.000000e+04 5.623413e-05 120 12
+1.000000e-05 1.000000e+05 1.000000e+04 3.162278e-05 60 6
+1.000000e-05 1.000000e+05 1.000000e+04 3.162278e-05 80 8
+1.000000e-05 1.000000e+05 1.000000e+04 3.162278e-05 100 10
+1.000000e-05 1.000000e+05 1.000000e+04 3.162278e-05 120 12
+1.000000e-05 1.000000e+05 1.000000e+04 1.778279e-05 60 6
+1.000000e-05 1.000000e+05 1.000000e+04 1.778279e-05 80 8
+1.000000e-05 1.000000e+05 1.000000e+04 1.778279e-05 100 10
+1.000000e-05 1.000000e+05 1.000000e+04 1.778279e-05 120 12
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-05 60 6
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-05 80 8
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-05 100 10
+1.000000e-05 1.000000e+05 1.000000e+04 1.000000e-05 120 12
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-04 60 6
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-04 80 8
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-04 100 10
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-04 120 12
+1.000000e-05 1.778279e+05 1.778279e+04 5.623413e-05 60 6
+1.000000e-05 1.778279e+05 1.778279e+04 5.623413e-05 80 8
+1.000000e-05 1.778279e+05 1.778279e+04 5.623413e-05 100 10
+1.000000e-05 1.778279e+05 1.778279e+04 5.623413e-05 120 12
+1.000000e-05 1.778279e+05 1.778279e+04 3.162278e-05 60 6
+1.000000e-05 1.778279e+05 1.778279e+04 3.162278e-05 80 8
+1.000000e-05 1.778279e+05 1.778279e+04 3.162278e-05 100 10
+1.000000e-05 1.778279e+05 1.778279e+04 3.162278e-05 120 12
+1.000000e-05 1.778279e+05 1.778279e+04 1.778279e-05 60 6
+1.000000e-05 1.778279e+05 1.778279e+04 1.778279e-05 80 8
+1.000000e-05 1.778279e+05 1.778279e+04 1.778279e-05 100 10
+1.000000e-05 1.778279e+05 1.778279e+04 1.778279e-05 120 12
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-05 60 6
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-05 80 8
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-05 100 10
+1.000000e-05 1.778279e+05 1.778279e+04 1.000000e-05 120 12
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-04 60 6
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-04 80 8
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-04 100 10
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-04 120 12
+1.000000e-05 3.162278e+05 3.162278e+04 5.623413e-05 60 6
+1.000000e-05 3.162278e+05 3.162278e+04 5.623413e-05 80 8
+1.000000e-05 3.162278e+05 3.162278e+04 5.623413e-05 100 10
+1.000000e-05 3.162278e+05 3.162278e+04 5.623413e-05 120 12
+1.000000e-05 3.162278e+05 3.162278e+04 3.162278e-05 60 6
+1.000000e-05 3.162278e+05 3.162278e+04 3.162278e-05 80 8
+1.000000e-05 3.162278e+05 3.162278e+04 3.162278e-05 100 10
+1.000000e-05 3.162278e+05 3.162278e+04 3.162278e-05 120 12
+1.000000e-05 3.162278e+05 3.162278e+04 1.778279e-05 60 6
+1.000000e-05 3.162278e+05 3.162278e+04 1.778279e-05 80 8
+1.000000e-05 3.162278e+05 3.162278e+04 1.778279e-05 100 10
+1.000000e-05 3.162278e+05 3.162278e+04 1.778279e-05 120 12
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-05 60 6
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-05 80 8
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-05 100 10
+1.000000e-05 3.162278e+05 3.162278e+04 1.000000e-05 120 12
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-04 60 6
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-04 80 8
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-04 100 10
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-04 120 12
+1.000000e-05 5.623413e+05 5.623413e+04 5.623413e-05 60 6
+1.000000e-05 5.623413e+05 5.623413e+04 5.623413e-05 80 8
+1.000000e-05 5.623413e+05 5.623413e+04 5.623413e-05 100 10
+1.000000e-05 5.623413e+05 5.623413e+04 5.623413e-05 120 12
+1.000000e-05 5.623413e+05 5.623413e+04 3.162278e-05 60 6
+1.000000e-05 5.623413e+05 5.623413e+04 3.162278e-05 80 8
+1.000000e-05 5.623413e+05 5.623413e+04 3.162278e-05 100 10
+1.000000e-05 5.623413e+05 5.623413e+04 3.162278e-05 120 12
+1.000000e-05 5.623413e+05 5.623413e+04 1.778279e-05 60 6
+1.000000e-05 5.623413e+05 5.623413e+04 1.778279e-05 80 8
+1.000000e-05 5.623413e+05 5.623413e+04 1.778279e-05 100 10
+1.000000e-05 5.623413e+05 5.623413e+04 1.778279e-05 120 12
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-05 60 6
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-05 80 8
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-05 100 10
+1.000000e-05 5.623413e+05 5.623413e+04 1.000000e-05 120 12
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-04 60 6
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-04 80 8
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-04 100 10
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-04 120 12
+1.000000e-05 1.000000e+06 1.000000e+05 5.623413e-05 60 6
+1.000000e-05 1.000000e+06 1.000000e+05 5.623413e-05 80 8
+1.000000e-05 1.000000e+06 1.000000e+05 5.623413e-05 100 10
+1.000000e-05 1.000000e+06 1.000000e+05 5.623413e-05 120 12
+1.000000e-05 1.000000e+06 1.000000e+05 3.162278e-05 60 6
+1.000000e-05 1.000000e+06 1.000000e+05 3.162278e-05 80 8
+1.000000e-05 1.000000e+06 1.000000e+05 3.162278e-05 100 10
+1.000000e-05 1.000000e+06 1.000000e+05 3.162278e-05 120 12
+1.000000e-05 1.000000e+06 1.000000e+05 1.778279e-05 60 6
+1.000000e-05 1.000000e+06 1.000000e+05 1.778279e-05 80 8
+1.000000e-05 1.000000e+06 1.000000e+05 1.778279e-05 100 10
+1.000000e-05 1.000000e+06 1.000000e+05 1.778279e-05 120 12
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-05 60 6
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-05 80 8
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-05 100 10
+1.000000e-05 1.000000e+06 1.000000e+05 1.000000e-05 120 12
diff --git a/Metafor/models/bord01/numericalSAWeights2.ascii b/Metafor/models/bord01/numericalSAWeights2.ascii
new file mode 100644
index 00000000..5e6de452
--- /dev/null
+++ b/Metafor/models/bord01/numericalSAWeights2.ascii
@@ -0,0 +1,500 @@
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
+ 1.0000000000000000e+00
diff --git a/Metafor/models/bord01/tombeBord.py b/Metafor/models/bord01/tombeBord.py
index d316e3cf..a685f9b1 100644
--- a/Metafor/models/bord01/tombeBord.py
+++ b/Metafor/models/bord01/tombeBord.py
@@ -1,4 +1,4 @@
-# -*- coding: latin-1; -*-
+# -*- coding: latin-1; -*-
# $Id: tombeBordEas2D.py 1422 2011-03-24 08:38:18Z papeleux $
#################################################
# Tombe de bord #
@@ -22,13 +22,26 @@ def getMetafor(_parameters={}):
def buildMetafor(p={}):
#-- parametres par defaut :
- parameters = {'NIPDETA': 6, 'Young':210000.0,'Poisson':0.3,
- 'sigEl0':420.0,'ih_H':1500.0,
- 'npgdy': 16,'nx':35,'ny':8,
- 'Thickness':1.0,'Rm':3.0,'Rt':1.0,
- 'Gap':1.0,'nbArch':1,
- 'peno':1e+08, 'mu_sta':0.6, 'mu_dyn':0.6,
- 'NRTol':1e-07}
+ parameters = {
+ 'NIPDETA' : 2, # nbre de pts de Gauss selon l'epaisseur
+ 'Young' : 210000.0,
+ 'Poisson' : 0.3,
+ 'sigEl0' : 420.0,
+ 'ih_H' : 1500.0,
+ 'nx' : 60,
+ 'ny' : 6,
+ 'Thickness' : 1.0,
+ 'Rm' : 3.0, # rayon matrice
+ 'Rt' : 1.0, # rayon outil
+ 'Gap' : 1.0,
+ 'nbArch' : 1,
+ 'peno' : 5e+05,
+ 'peta' : 5e+04,
+ 'mu_sta' : 0.1,
+ 'mu_dyn' : 0.1,
+ 'NRTol' : 1e-04,
+ 'PEAS' : 1e-06
+ }
print p
parameters.update(p)
@@ -47,13 +60,13 @@ def buildMetafor(p={}):
nx2 = parameters['nx']
ny = parameters['ny']
# Tool
- Rm = parameters['Rm'] #Rm : rayon Matrice
- Rt = parameters['Rt'] #Rt : rayon Tool
- Gap = parameters['Gap'] #Jeu matrice - poincon
- Lp = 3.0 #longueur du poincon
+ Rm = parameters['Rm'] # Rm : rayon Matrice
+ Rt = parameters['Rt'] # Rt : rayon Tool
+ Gap = parameters['Gap'] # Jeu matrice - poincon
+ Lp = 3.0 # longueur du poincon
#-- maillage --
- Lx1 = -2*X0 + (math.pi/2.0) * Rm
+ Lx1 = Lx/2.0 #-2*X0 + (math.pi/2.0) * Rm
Lx2 = Lx - Lx1
from toolbox.meshedGeometry2D import createRectangleSameId
createRectangleSameId(domain, nx1, ny, Lx1, Th, X0, Y0, 0)
@@ -80,24 +93,28 @@ def buildMetafor(p={}):
wireset.add( Wire(201, [curveset(1), curveset(101)]) )
wireset.add( Wire(202, [curveset(3), curveset(103)]) )
+
#Outil 1 : Matrice
ci1 = RdContactInteraction(1)
ci1.setTool(wireset(1001))
ci1.push(wireset(201))
ci1.addProperty(prp_c1)
interactionset.add(ci1)
+
#Outil 2 : Poincon
ci2 = RdContactInteraction(2)
ci2.setTool(wireset(1101))
ci2.push(wireset(202))
ci2.addProperty(prp_c1)
interactionset.add(ci2)
+
#Outil 3 : Serre-Flanc
ci3 = RdContactInteraction(3)
ci3.setTool(wireset(1201))
ci3.push(curveset(3))
ci3.addProperty(prp_c1)
interactionset.add(ci3)
+
#-- FieldApplicator --
app001 = FieldApplicator(99)
app001.push(sideset( 1))
@@ -106,6 +123,7 @@ def buildMetafor(p={}):
app101 = FieldApplicator(101)
app101.push(sideset(101))
interactionset.add(app101)
+
#-- Fnn - Loadings --
loadingset = domain.getLoadingSet()
loadingset.define(curveset(4), Field1D(TX,RE))
@@ -133,7 +151,7 @@ def buildMetafor(p={}):
materset.define(1001, CoulombContactMaterial)
mat1001 = materset(1001)
mat1001.put(PEN_NORMALE, parameters['peno'])
- mat1001.put(PEN_TANGENT, 1e+06)
+ mat1001.put(PEN_TANGENT, parameters['peta'])
mat1001.put(PROF_CONT, 0.9)
mat1001.put(COEF_FROT_STA, parameters['mu_sta']) # usually between 0.3 and 0.6 for dry materials
mat1001.put(COEF_FROT_DYN, parameters['mu_dyn'])
@@ -144,7 +162,8 @@ def buildMetafor(p={}):
prp001.put(MATERIAL, 1)
prp001.put(NIPDETA, parameters['NIPDETA'])
prp001.put(EASS, 2)
- prp001.put(EASV, 4)
+ prp001.put(EASV, 2)
+ prp001.put(PEAS, parameters['PEAS'])
prp001.put(VERBOSE, True)
interactionset(99).addProperty(prp001)
@@ -153,28 +172,28 @@ def buildMetafor(p={}):
prp101.put(MATERIAL, 1)
prp101.put(NIPDETA, parameters['NIPDETA'])
prp101.put(EASS, 2)
- prp101.put(EASV, 4)
+ prp101.put(EASV, 2)
+ prp101.put(PEAS, parameters['PEAS'])
prp101.put(VERBOSE, True)
interactionset(101).addProperty(prp101)
#
mim = metafor.getMechanicalIterationManager()
+ mim.setResidualComputationMethod(Method4ResidualComputation())
mim.setResidualTolerance(parameters['NRTol'])
#-- tsm --
tsm = metafor.getTimeStepManager()
- tsm.setInitialTime(0.0, 0.1)
+ tsm.setInitialTime(0.0, 1e-4)
tsm.setNextTime(10.0, parameters['nbArch'], 1.0)
tsm.setNextTime(20.0, parameters['nbArch'], 1.0)
- #Time Integration
+ # Time Integration
ti = AlphaGeneralizedTimeIntegration(metafor)
- ti.setAlphaM(-0.97)
- ti.setAlphaF(0.01)
metafor.setTimeIntegration(ti)
- #Courbes
+ # Courbes
valuesmanager = metafor.getValuesManager()
pointset = geometry.getPointSet()
curveset = geometry.getCurveSet()
@@ -182,8 +201,8 @@ def buildMetafor(p={}):
ave1 = AngleValueExtractor(Axe(curveset(1003)),Axe(pointset(2),pointset(102)))
valuesmanager.add(11,ave1,'Angle1')
- # Objective Function :
- #---------------------
+ # Objective Function
+ # ------------------
objectivefunctionset = metafor.getObjectiveFunctionSet()
Index = 2000
@@ -200,7 +219,7 @@ def buildMetafor(p={}):
sinAlphaDemi = math.sin(alphaRad / 2.0)
cosAlpha = math.cos(alphaRad)
sinAlpha = math.sin(alphaRad)
- cos90moinsAlpha = math.cos(math.pi / 2.0 - alphaRad)
+ cos90moinsAlpha = math.cos(math.pi / 2.0 - alphaRad)
sinAlpha = math.sin(alphaRad)
from toolbox.domainTools import getGeoReferences
@@ -208,6 +227,7 @@ def buildMetafor(p={}):
piDemi = math.asin(1)
C45 = math.cos(piDemi/2.0)
Z0 = 0.0
+
# Outil 1 : matrice
pointset.define(Index + 1, X0 - dX0, Y0, Z0)
pointset.define(Index + 2, X0, Y0, Z0)
diff --git a/Metafor/models/bord01/tombeBord_old.py b/Metafor/models/bord01/tombeBord_old.py
new file mode 100644
index 00000000..d316e3cf
--- /dev/null
+++ b/Metafor/models/bord01/tombeBord_old.py
@@ -0,0 +1,247 @@
+# -*- coding: latin-1; -*-
+# $Id: tombeBordEas2D.py 1422 2011-03-24 08:38:18Z papeleux $
+#################################################
+# Tombe de bord #
+#===============================================#
+#################################################
+from wrap import *
+import math
+
+metafor = None
+
+def getMetafor(_parameters={}):
+ global metafor
+ if not metafor :
+ '''
+ params = {'npgdy':50}
+ params.update(_parameters)
+ metafor = buildMetafor(params)
+ '''
+ metafor = buildMetafor(_parameters)
+ return metafor
+
+def buildMetafor(p={}):
+ #-- parametres par defaut :
+ parameters = {'NIPDETA': 6, 'Young':210000.0,'Poisson':0.3,
+ 'sigEl0':420.0,'ih_H':1500.0,
+ 'npgdy': 16,'nx':35,'ny':8,
+ 'Thickness':1.0,'Rm':3.0,'Rt':1.0,
+ 'Gap':1.0,'nbArch':1,
+ 'peno':1e+08, 'mu_sta':0.6, 'mu_dyn':0.6,
+ 'NRTol':1e-07}
+
+ print p
+ parameters.update(p)
+ metafor = Metafor()
+ domain = metafor.getDomain()
+
+ geometry = domain.getGeometry()
+ geometry.setDimPlaneStrain(1.0)
+ #-- dimension du maillage
+ X0 = -3.0
+ Y0 = 0.0
+ Z0 = 0.0
+ Lx = 20.0
+ Th = parameters['Thickness']
+ nx1 = parameters['nx']
+ nx2 = parameters['nx']
+ ny = parameters['ny']
+ # Tool
+ Rm = parameters['Rm'] #Rm : rayon Matrice
+ Rt = parameters['Rt'] #Rt : rayon Tool
+ Gap = parameters['Gap'] #Jeu matrice - poincon
+ Lp = 3.0 #longueur du poincon
+
+ #-- maillage --
+ Lx1 = -2*X0 + (math.pi/2.0) * Rm
+ Lx2 = Lx - Lx1
+ from toolbox.meshedGeometry2D import createRectangleSameId
+ createRectangleSameId(domain, nx1, ny, Lx1, Th, X0, Y0, 0)
+ from toolbox.meshedGeometry2D import addQuadrangle2X
+ addQuadrangle2X(domain, 100, nx2, ny, Lx2, Th, X0+Lx1, Y0,
+ 2, 3, 2)
+ #-- Fonctions d'evolution --
+ FCT_CHA = PieceWiseLinearFunction()
+ FCT_CHA.setData(0.0, 0.0)
+ FCT_CHA.setData(10.0, 1.0)
+ FCT_CHA.setData(20.0 , 0.0)
+ #Outils de mise a forme
+ from apps.toolbox.formingTools import createFlangingTools2D
+ createFlangingTools2D(domain, 1000, (X0-1), Y0, Rm, Rt, Gap, Th, Lp, Lx,-Lx, FCT_CHA)
+
+ prp_c1 = ElementProperties(Contact2DElement)
+ prp_c1.put(MATERIAL, 1001)
+ prp_c1.put(AREAINCONTACT,AIC_ONCE)
+
+ interactionset = domain.getInteractionSet()
+ curveset = geometry.getCurveSet()
+ wireset = geometry.getWireSet()
+ sideset = geometry.getSideSet()
+
+ wireset.add( Wire(201, [curveset(1), curveset(101)]) )
+ wireset.add( Wire(202, [curveset(3), curveset(103)]) )
+ #Outil 1 : Matrice
+ ci1 = RdContactInteraction(1)
+ ci1.setTool(wireset(1001))
+ ci1.push(wireset(201))
+ ci1.addProperty(prp_c1)
+ interactionset.add(ci1)
+ #Outil 2 : Poincon
+ ci2 = RdContactInteraction(2)
+ ci2.setTool(wireset(1101))
+ ci2.push(wireset(202))
+ ci2.addProperty(prp_c1)
+ interactionset.add(ci2)
+ #Outil 3 : Serre-Flanc
+ ci3 = RdContactInteraction(3)
+ ci3.setTool(wireset(1201))
+ ci3.push(curveset(3))
+ ci3.addProperty(prp_c1)
+ interactionset.add(ci3)
+ #-- FieldApplicator --
+ app001 = FieldApplicator(99)
+ app001.push(sideset( 1))
+ interactionset.add(app001)
+
+ app101 = FieldApplicator(101)
+ app101.push(sideset(101))
+ interactionset.add(app101)
+ #-- Fnn - Loadings --
+ loadingset = domain.getLoadingSet()
+ loadingset.define(curveset(4), Field1D(TX,RE))
+ loadingset.define(curveset(4), Field1D(TY,RE))
+
+ #-- Materiaux --
+ materset = domain.getMaterialSet()
+ materlawset = domain.getMaterialLawSet()
+ materset.define(1, EvpIsoHHypoMaterial)
+ materset(1).put(MASS_DENSITY, 8900.0E-12)
+ materset(1).put(ELASTIC_MODULUS, parameters['Young'])
+ materset(1).put(POISSON_RATIO, parameters['Poisson'])
+ materset(1).put(YIELD_NUM, 2)
+ '''
+ materlawset.define(1, RambergOsgoodIsotropicHardening)
+ materlawset(1).put(IH_SIGEL, 330.3)
+ materlawset(1).put(IH_A, 591.72)
+ materlawset(1).put(IH_N, 5.35)
+ '''
+ materlawset.define(2, LinearIsotropicHardening)
+ materlawset(2).put(IH_SIGEL, parameters['sigEl0'])
+ materlawset(2).put(IH_H, parameters['ih_H'])
+
+
+ materset.define(1001, CoulombContactMaterial)
+ mat1001 = materset(1001)
+ mat1001.put(PEN_NORMALE, parameters['peno'])
+ mat1001.put(PEN_TANGENT, 1e+06)
+ mat1001.put(PROF_CONT, 0.9)
+ mat1001.put(COEF_FROT_STA, parameters['mu_sta']) # usually between 0.3 and 0.6 for dry materials
+ mat1001.put(COEF_FROT_DYN, parameters['mu_dyn'])
+
+ #-- Propelem --
+ prp001 = ElementProperties(Volume2DElement)
+ prp001.put(CAUCHYMECHVOLINTMETH, VES_CMVIM_EAS)
+ prp001.put(MATERIAL, 1)
+ prp001.put(NIPDETA, parameters['NIPDETA'])
+ prp001.put(EASS, 2)
+ prp001.put(EASV, 4)
+ prp001.put(VERBOSE, True)
+ interactionset(99).addProperty(prp001)
+
+ prp101 = ElementProperties(Volume2DElement)
+ prp101.put(CAUCHYMECHVOLINTMETH, VES_CMVIM_EAS)
+ prp101.put(MATERIAL, 1)
+ prp101.put(NIPDETA, parameters['NIPDETA'])
+ prp101.put(EASS, 2)
+ prp101.put(EASV, 4)
+ prp101.put(VERBOSE, True)
+ interactionset(101).addProperty(prp101)
+
+ #
+
+ mim = metafor.getMechanicalIterationManager()
+ mim.setResidualTolerance(parameters['NRTol'])
+
+ #-- tsm --
+ tsm = metafor.getTimeStepManager()
+ tsm.setInitialTime(0.0, 0.1)
+ tsm.setNextTime(10.0, parameters['nbArch'], 1.0)
+ tsm.setNextTime(20.0, parameters['nbArch'], 1.0)
+
+ #Time Integration
+ ti = AlphaGeneralizedTimeIntegration(metafor)
+ ti.setAlphaM(-0.97)
+ ti.setAlphaF(0.01)
+ metafor.setTimeIntegration(ti)
+
+ #Courbes
+ valuesmanager = metafor.getValuesManager()
+ pointset = geometry.getPointSet()
+ curveset = geometry.getCurveSet()
+ valuesmanager.add(1, MiscValueExtractor(metafor,EXT_T),'time')
+ ave1 = AngleValueExtractor(Axe(curveset(1003)),Axe(pointset(2),pointset(102)))
+ valuesmanager.add(11,ave1,'Angle1')
+
+ # Objective Function :
+ #---------------------
+ objectivefunctionset = metafor.getObjectiveFunctionSet()
+
+ Index = 2000
+ X0 = -10.0
+ dX0 = 3.0
+ Y0 = -10.0
+ Z0 = 0.0
+ Rm = 4.0
+ Ll = 50.0
+ alphaDeg = 80.0
+
+ alphaRad = alphaDeg * math.pi / 180.0
+ cosAlphaDemi = math.cos(alphaRad / 2.0)
+ sinAlphaDemi = math.sin(alphaRad / 2.0)
+ cosAlpha = math.cos(alphaRad)
+ sinAlpha = math.sin(alphaRad)
+ cos90moinsAlpha = math.cos(math.pi / 2.0 - alphaRad)
+ sinAlpha = math.sin(alphaRad)
+
+ from toolbox.domainTools import getGeoReferences
+ [pointset, curveset, wireset, surfaceset, sideset, skinset, volumeset] = getGeoReferences(domain)
+ piDemi = math.asin(1)
+ C45 = math.cos(piDemi/2.0)
+ Z0 = 0.0
+ # Outil 1 : matrice
+ pointset.define(Index + 1, X0 - dX0, Y0, Z0)
+ pointset.define(Index + 2, X0, Y0, Z0)
+ pointset.define(Index + 3, X0 + (Rm * sinAlphaDemi), (Y0 - Rm * (1.0-cosAlphaDemi)), Z0)
+ pointset.define(Index + 4, X0 + (Rm * sinAlpha), (Y0 - Rm * (1.0-cosAlpha)), Z0)
+ pointset.define(Index + 5, X0 + (Rm * sinAlpha + Ll * cosAlpha), (Y0 - Rm * (1.0 - cosAlpha) - Ll * sinAlpha), Z0)
+ pointset.define(Index + 6, X0 - dX0 , (Y0 + 2.0), Z0)
+
+ curveset.add( Line(Index + 1, pointset(Index + 1), pointset(Index + 2) ))
+ curveset.add( Arc(Index + 2, pointset(Index + 2), pointset(Index + 3), pointset(Index + 4) ))
+ curveset.add( Line(Index + 3, pointset(Index + 4), pointset(Index + 5) ))
+ curveset.add( Line(Index + 4, pointset(Index + 6), pointset(Index + 1) ))
+
+ wireset.add(MultiProjWire(Index + 1, [curveset(Index+i) for i in (1, 2, 3)]))
+ wireset.add(Wire(Index+2, [curveset(1), curveset(101)]))
+
+ sve1 = ShapeValueExtractor(wireset(Index+1), wireset(Index+2))
+ sve1.setTriedreRef(Triedre(pointset(Index+1),pointset(Index+2),pointset(Index+6)))
+ sve1.setTriedreMesh(Triedre(pointset(1),pointset(1002),pointset(4)))
+
+ valuesmanager.add(12,sve1,'shape1')
+ valuesmanager.add(13,sve1,NormOperator(),'shape2')
+
+ sof1 = ValueExtractorObjectiveFunction(1,ave1)
+ objectivefunctionset.add(sof1)
+ sof2 = ValueExtractorObjectiveFunction(2,sve1)
+ objectivefunctionset.add(sof2)
+
+ # Verification des Objective Function in .res File
+ testSuite = metafor.getTestSuiteChecker()
+ testSuite.checkObjectiveFunction(2)
+ testSuite.checkExtractor(11)
+ testSuite.checkExtractor(12,0)
+ testSuite.checkExtractor(12,5)
+ testSuite.checkExtractor(13)
+
+ return metafor
diff --git a/mrstlnos/matlab/Belos/test_belos8.m b/mrstlnos/matlab/Belos/test_belos8.m
new file mode 100644
index 00000000..a157e5f1
--- /dev/null
+++ b/mrstlnos/matlab/Belos/test_belos8.m
@@ -0,0 +1,25 @@
+%close all
+clc
+clear all
+
+[it_1g, res_1g] = read_belos('multiFGuess_belos.text');
+[it_2g, res_2g] = read_belos('multiFGuess_belos_preconditioned.text');
+
+%%
+close all
+figure
+semilogy(it_1g, res_1g,'o','LineWidth',1,...
+ 'MarkerSize',8,...
+ 'MarkerEdgeColor','k',...
+ 'MarkerFaceColor','r')
+hold on
+semilogy(it_2g, res_2g,'o','LineWidth',1,...
+ 'MarkerSize',8,...
+ 'MarkerEdgeColor','k',...
+ 'MarkerFaceColor','g')
+grid on
+box on
+set(gca,'fontsize',18)
+xlabel('Iteration')
+ylabel('(2-Norm Res Vec) / (2-Norm Prec Res0)')
+hold off
diff --git a/mrstlnos/matlab/Multigrid/read_agregate.m b/mrstlnos/matlab/Multigrid/read_agregate.m
new file mode 100644
index 00000000..051fa1e6
--- /dev/null
+++ b/mrstlnos/matlab/Multigrid/read_agregate.m
@@ -0,0 +1,20 @@
+function [indices] = read_agregate(filename,agg_id,proc)
+
+ indices = [];
+
+ fid = fopen(filename);
+
+ line_start = ['Agg ',int2str(agg_id), ' Proc ', int2str(proc),':'];
+ tline = fgetl(fid);
+ while ischar(tline)
+ current_length = length(tline);
+
+ if tline(1:length(line_start)) == line_start
+ indices = str2num(tline((length(line_start)+1):current_length));
+
+ break
+ end
+ tline = fgetl(fid);
+ end
+ fclose(fid);
+end
\ No newline at end of file
diff --git a/mrstlnos/matlab/Multigrid/testOneDProlongation.m b/mrstlnos/matlab/Multigrid/testOneDProlongation.m
index b9ce6db2..dfbfc579 100644
--- a/mrstlnos/matlab/Multigrid/testOneDProlongation.m
+++ b/mrstlnos/matlab/Multigrid/testOneDProlongation.m
@@ -1,30 +1,96 @@
-
clear all
close all
clc
format long
L = 1;
-n = 7;
-[ A,Ah,H ] = oneDLaplacianMatrix( n, L );
+nH = 7;
+n = 2*nH+1;
+mode = 4;
+[ Ah,Ahh,h ] = oneDLaplacianMatrix( n, L );
+
+
+
+H = 2 *h;
+
+[I_h_H,n2] = oneDRestriction(n);
+[I_H_h,n2] = oneDProlongation(n);
+AH = I_h_H * Ah * I_H_h;
+
+
-[V,D] = eig(A);
+[Vh,D] = eig(Ah);
d = diag(D);
+[VH,DH] = eig(AH);
-[I_H_h,n2] = oneDProlongation(n*2+1);
-xH = V(:,1);
+V = zeros(n,n);
+for k =1:n
+ for i =1:n
+ V(i,k) = sin(i*k*pi/(n+1));
+ end
+end
+
+
+[Vh,Dh] = eig(Ah);
+
+tmp = diag(V'*Vh);
+for i=1:n
+ if tmp(i) < 0
+ Vh(:,i) = -Vh(:,i);
+ end
+end
+tmp = diag(VH'*I_h_H*Vh);
+for i=1:((n-1)/2)
+ VH(:,i) = VH(:,i)/tmp(i); % change the norm
+end
+
+xH = VH(:,mode);
maxx = max(abs(xH));
xH = xH / maxx;
xh = I_H_h * xH;
figure
-subplot(1,2,1)
+subplot(2,3,1)
+plot([0:H:L],[0;xH;0],'o-','linewidth',2)
+
+box on
+grid on
+subplot(2,3,2)
+plot([0:h:L],[0;xh;0],'o-','linewidth',2)
+box on
+grid on
+subplot(2,3,3)
+plot([0:h:L],[0;xh;0],'o-','linewidth',2)
+hold on
plot([0:H:L],[0;xH;0],'o-','linewidth',2)
+plot([0:H:L],[0;VH(:,mode);0]*norm([0;xH;0])/norm([0;VH(:,mode);0]),'o--','linewidth',2)
+box on
+grid on
+subplot(2,3,4)
+
+plot(0:((n-1)/2-1),VH'*xH,'o-','linewidth',2)
box on
grid on
-subplot(1,2,2)
-plot([0:H/2:L],[0;xh;0],'o-','linewidth',2)
+subplot(2,3,5)
+plot(0:(n-1),Vh'*xh,'o-','linewidth',2)
box on
-grid on
\ No newline at end of file
+grid on
+subplot(2,3,6)
+plot(0:(n-1),Vh'*xh,'o-','linewidth',2)
+hold on
+plot(0:((n-1)/2-1),VH'*xH,'o-','linewidth',2)
+plot(0:((n-1)/2-1),VH'*VH(:,mode)*norm([0;xH;0])/norm([0;VH(:,mode);0]),'o--','linewidth',2)
+
+box on
+grid on
+
+figure
+tmp = Vh'*xh;
+subplot(2,2,1)
+plot(0:(nH-1),VH'*xH,'o-','linewidth',2)
+subplot(2,2,2)
+plot(0:(nH-1),tmp(1:nH),'o-','linewidth',2)
+subplot(2,2,4)
+plot(0:(nH-1),tmp(n:-1:n-nH+1),'o-','linewidth',2)
diff --git a/mrstlnos/matlab/Multigrid/testOneDRestriction.m b/mrstlnos/matlab/Multigrid/testOneDRestriction.m
index ff1bd783..4680a027 100644
--- a/mrstlnos/matlab/Multigrid/testOneDRestriction.m
+++ b/mrstlnos/matlab/Multigrid/testOneDRestriction.m
@@ -1,32 +1,98 @@
clear all
-%close all
+close all
clc
format long
L = 1;
n = 15;
-[ A,Ah,h ] = oneDLaplacianMatrix( n, L );
+nH = (n-1)/2;
+mode = 12;
+[ Ah,Ahh,h ] = oneDLaplacianMatrix( n, L );
+
-[V,D] = eig(A);
-d = diag(D);
H = 2 *h;
[I_h_H,n2] = oneDRestriction(n);
+[I_H_h,n2] = oneDProlongation(n);
+AH = I_h_H * Ah * I_H_h;
+
+
+
+[Vh,D] = eig(Ah);
+d = diag(D);
+[VH,DH] = eig(AH);
-xh = V(:,1);
+
+V = zeros(n,n);
+for k =1:n
+ for i =1:n
+ V(i,k) = sin(i*k*pi/(n+1));
+ end
+end
+
+
+[Vh,Dh] = eig(Ah);
+
+tmp = diag(V'*Vh);
+for i=1:n
+ if tmp(i) < 0
+ Vh(:,i) = -Vh(:,i);
+ end
+end
+tmp = diag(VH'*I_h_H*Vh);
+for i=1:((n-1)/2)
+ VH(:,i) = VH(:,i)/tmp(i); % change the norm
+end
+
+xh = Vh(:,mode);
maxx = max(abs(xh));
xh = xh / maxx;
xH = I_h_H * xh;
figure
-subplot(1,2,1)
+subplot(2,3,1)
plot([0:h:L],[0;xh;0],'o-','linewidth',2)
box on
grid on
-subplot(1,2,2)
+subplot(2,3,2)
+plot([0:H:L],[0;xH;0],'o-','linewidth',2)
+box on
+grid on
+subplot(2,3,3)
+plot([0:h:L],[0;xh;0],'o-','linewidth',2)
+hold on
plot([0:H:L],[0;xH;0],'o-','linewidth',2)
+if mode <= size(VH,2)
+ plot([0:H:L],[0;VH(:,mode);0]*norm([0;xH;0])/norm([0;VH(:,mode);0]),'o--','linewidth',2)
+end
+box on
+grid on
+subplot(2,3,4)
+plot(0:(n-1),Vh'*xh,'o-','linewidth',2)
+box on
+grid on
+subplot(2,3,5)
+plot(0:((n-1)/2-1),VH'*xH,'o-','linewidth',2)
+box on
+grid on
+subplot(2,3,6)
+plot(0:(n-1),Vh'*xh,'o-','linewidth',2)
+hold on
+plot(0:((n-1)/2-1),VH'*xH,'o-','linewidth',2)
+if mode <= size(VH,2)
+ plot(0:((n-1)/2-1),VH'*VH(:,mode)*norm([0;xH;0])/norm([0;VH(:,mode);0]),'o--','linewidth',2)
+end
box on
-grid on
\ No newline at end of file
+grid on
+
+figure
+tmp = Vh'*xh;
+subplot(2,2,1)
+plot(0:(nH-1),tmp(1:nH),'o-','linewidth',2)
+subplot(2,2,3)
+plot(0:(nH-1),tmp(n:-1:n-nH+1),'o-','linewidth',2)
+subplot(2,2,2)
+plot(0:(nH-1),VH'*xH,'o-','linewidth',2)
\ No newline at end of file
diff --git a/mrstlnos/src/IterativeSolver.cpp b/mrstlnos/src/IterativeSolver.cpp
index 2ee9cd57..b2bd245a 100644
--- a/mrstlnos/src/IterativeSolver.cpp
+++ b/mrstlnos/src/IterativeSolver.cpp
@@ -142,7 +142,7 @@ void IterativeSolver::start()
typedef Sacado::MP::Vector<Stokhos::StaticFixedStorage<int,double,4,Kokkos::Threads>> ensemble_scalar;
IterativeSolver::tstart<ensemble_scalar>();
}
- /*
+ /*
else if (minimum_ensemble_size <= 8)
{
typedef Sacado::MP::Vector<Stokhos::StaticFixedStorage<int,double,8,Kokkos::Threads>> ensemble_scalar;
@@ -321,9 +321,17 @@ template<typename scalar> void IterativeSolver::tstart()
LinearSolver::getTimers()["prec"].start();
std::string muelu_output_file = solverList.get("Output Stream file name muelu","muelu_out.txt");
Teuchos::ParameterList * mueluParams = new Teuchos::ParameterList(solverList.get("mueluParams",Teuchos::ParameterList()));
+ RCP<Tpetra::MultiVector<double > > inCoords (new Tpetra::MultiVector<double >(discrete_problem.algebraic.map.mapDofs, 3, true));
+ /*for (auto i = 0; i < discrete_problem.algebraic.map.mapNodes->getNodeNumElements(); ++i)
+ for( auto j=0; j<3; ++j)
+ {
+ auto global_index = discrete_problem.algebraic.map.mapNodes->getGlobalElement(i);
+ inCoords->replaceGlobalValue(global_index ,j,(double) discrete_problem.domain.nodes_list.nodes(discrete_problem.algebraic.map.mapNodesWO->getLocalElement(global_index),j));
+ }
+ */
RCP<Tpetra::MultiVector<scalar >> nullspace = Teuchos::null;
RCP<Tpetra::Operator<scalar >> M;
- M = mueluPreconditioner<scalar>(*mueluParams,discrete_problem.algebraic.matrices.A,nullspace, muelu_output_file);
+ M = mueluPreconditioner<scalar>(*mueluParams,discrete_problem.algebraic.matrices.A, inCoords, nullspace, muelu_output_file);
LinearSolver::getTimers()["prec"].stop();
LinearSolver::getTimers()["belos"].start();
blinproblem->setRightPrec(M);
diff --git a/mrstlnos/src/Preconditioner.h b/mrstlnos/src/Preconditioner.h
index 17cd81e3..3b70e5f8 100644
--- a/mrstlnos/src/Preconditioner.h
+++ b/mrstlnos/src/Preconditioner.h
@@ -51,7 +51,7 @@
#include "Stokhos_MueLu_MP_Vector.hpp"
#include "MueLu_CreateTpetraPreconditioner.hpp"
-template<typename scalar> Teuchos::RCP<Tpetra::Operator<scalar> > mueluPreconditioner(Teuchos::ParameterList mueluParams, Teuchos::RCP<const Tpetra::Operator<scalar>> C, Teuchos::RCP<Tpetra::MultiVector<scalar>> nullspace, std::string muelu_output_file)
+template<typename scalar> Teuchos::RCP<Tpetra::Operator<scalar> > mueluPreconditioner(Teuchos::ParameterList mueluParams, Teuchos::RCP<const Tpetra::Operator<scalar>> C, Teuchos::RCP<Tpetra::MultiVector<double >> inCoords, Teuchos::RCP<Tpetra::MultiVector<scalar>> nullspace, std::string muelu_output_file)
{
typedef Tpetra::Vector<scalar> TP_Vec;
typedef typename TP_Vec::local_ordinal_type LO;
@@ -67,9 +67,8 @@ template<typename scalar> Teuchos::RCP<Tpetra::Operator<scalar> > mueluPrecondit
//Teuchos::RCP<Teuchos::FancyOStream> test = MueLu::BaseClass::GetDefaultOStream();
Teuchos::RCP<Tpetra::Operator<scalar >> ncC = Teuchos::rcp_const_cast<Tpetra::Operator<scalar >>(C);
- Teuchos::RCP<Tpetra::MultiVector<double >> inCoords = Teuchos::null;
- //Teuchos::RCP<muelu_tpetra_operator_type> M = MueLu::CreateTpetraPreconditioner<scalar,LO,GO,NT>(ncC, mueluParams, inCoords, nullspace);
- Teuchos::RCP<muelu_tpetra_operator_type> M = MueLu::CreateTpetraPreconditioner<scalar,LO,GO,NT>(ncC, mueluParams);
+ Teuchos::RCP<muelu_tpetra_operator_type> M = MueLu::CreateTpetraPreconditioner<scalar,LO,GO,NT>(ncC, mueluParams, inCoords, nullspace);
+ //Teuchos::RCP<muelu_tpetra_operator_type> M = MueLu::CreateTpetraPreconditioner<scalar,LO,GO,NT>(ncC, mueluParams);
return M;
}
diff --git a/mrstlnos/src/wExample2.cpp b/mrstlnos/src/wExample2.cpp
index bcdcbcb0..34cc954e 100644
--- a/mrstlnos/src/wExample2.cpp
+++ b/mrstlnos/src/wExample2.cpp
@@ -1,143 +1,608 @@
#include "wExample2.h"
-#include <Kokkos_Core.hpp>
-#include <cstdio>
-
#include "wTimer.h"
#include "wTimers.h"
-
-#include <Kokkos_Core.hpp>
-#include <Kokkos_Threads.hpp>
-#include <Threads/Kokkos_Threads_TaskPolicy.hpp>
-#include <cstdio>
-#include <cstdlib>
-#include <cmath>
+// Tpetra
+#include "Stokhos_Tpetra_MP_Vector.hpp"
+#include "Stokhos_Tpetra_Utilities_MP_Vector.hpp"
+#include "Tpetra_ConfigDefs.hpp"
+#include "Tpetra_DefaultPlatform.hpp"
+#include "Tpetra_Map.hpp"
+#include "Tpetra_MultiVector.hpp"
+#include "Tpetra_Vector.hpp"
+#include "Tpetra_CrsGraph.hpp"
+#include "Tpetra_CrsMatrix.hpp"
+#include "Stokhos_Tpetra_CG.hpp"
+
+
+#include "TpetraExt_MatrixMatrix.hpp"
+#include "Tpetra_MatrixIO.hpp"
+#include "Tpetra_DefaultPlatform.hpp"
+#include "Tpetra_Vector.hpp"
+#include "Tpetra_CrsMatrixMultiplyOp.hpp"
+#include "Tpetra_Import.hpp"
+#include "Tpetra_Export.hpp"
+#include "Teuchos_DefaultComm.hpp"
+#include "Teuchos_CommandLineProcessor.hpp"
+#include "Teuchos_XMLParameterListHelpers.hpp"
+#include "Teuchos_UnitTestHarness.hpp"
+#include "MatrixMarket_Tpetra.hpp"
+#include "Tpetra_RowMatrixTransposer.hpp"
+
+
+#include <Teuchos_Array.hpp>
+#include <Teuchos_GlobalMPISession.hpp>
+#include <Teuchos_oblackholestream.hpp>
+#include <Teuchos_ScalarTraits.hpp>
+#include <Teuchos_RCP.hpp>
+#include "Teuchos_ParameterList.hpp"
+
+#include "Stokhos_Amesos2_MP_Vector.hpp"
+#include "Amesos2_Factory.hpp"
+
+#include "Stokhos_Ifpack2_MP_Vector.hpp"
+#include "Ifpack2_Factory.hpp"
+
+#include "Belos_Tpetra_MP_Vector.hpp"
+#include "BelosLinearProblem.hpp"
+#include "BelosPseudoBlockGmresSolMgr.hpp"
+#include "BelosPseudoBlockCGSolMgr.hpp"
+#include "BelosBlockGmresSolMgr.hpp"
+
+#include "EnsembleTraits.h"
+
+#include <math.h>
using namespace fwk;
-// Type of a one-dimensional length-N array of int.
-typedef Kokkos::View<int*> view_type;
-typedef view_type::HostMirror host_view_type;
-// This is a "zero-dimensional" View, that is, a View of a single
-// value (an int, in this case). Access the value using operator()
-// with no arguments: e.g., 'count()'.
-//
-// Zero-dimensional Views are useful for reduction results that stay
-// resident in device memory, as well as for irregularly updated
-// shared state. We use it for the latter in this example.
-typedef Kokkos::View<int> count_type;
-typedef count_type::HostMirror host_count_type;
-
-
-// Functor for finding a list of primes in a given set of numbers. If
-// run in parallel, the order of results is nondeterministic, because
-// hardware atomic updates do not guarantee an order of execution.
-struct findprimes {
- view_type data;
- view_type result;
- count_type count;
-
- findprimes (view_type data_, view_type result_, count_type count_) :
- data (data_), result (result_), count (count_)
- {}
-
- // Test if data(i) is prime. If it is, increment the count of
- // primes (stored in the zero-dimensional View 'count') and add the
- // value to the current list of primes 'result'.
- KOKKOS_INLINE_FUNCTION
- void operator() (const int i) const {
- const int number = data(i); // the current number
-
- // Test all numbers from 3 to ceiling(sqrt(data(i))), to see if
- // they are factors of data(i). It's not the most efficient prime
- // test, but it works.
- const int upper_bound = sqrt(1.0*number)+1;
- bool is_prime = !(number%2 == 0);
- int k = 3;
- while (k < upper_bound && is_prime) {
- is_prime = !(number%k == 0);
- k += 2; // don't have to test even numbers
- }
+using namespace mrstlnos;
+
+typedef Tpetra::Vector<>::local_ordinal_type local_ordinal_type;
+typedef Tpetra::Vector<>::global_ordinal_type global_ordinal_type;
+
+template<typename scalar> void writeLocalCrsMatrixToFile(std::string &name, size_t myRank, Teuchos::RCP<Tpetra::CrsMatrix<scalar,local_ordinal_type,global_ordinal_type> > A)
+{
+ std::ofstream outputfile;
+ outputfile.open (name + "_" + std::to_string(myRank) + ".txt");
+
+ typedef std::numeric_limits< double > dbl;
+ outputfile.precision(dbl::max_digits10+2);
+
+ auto line_size = A->getRowMap()->getNodeNumElements ();
+
+ for (local_ordinal_type local_index = 0; local_index < line_size; ++local_index)
+ {
+ global_ordinal_type global_index = A->getRowMap()->getGlobalElement(local_index);
+
+ Teuchos::ArrayView<const local_ordinal_type> local_indices ;
+ Teuchos::ArrayView<const scalar> local_values ;
+ A->getLocalRowView (local_index, local_indices, local_values);
- if (is_prime) {
- // Use an atomic update both to update the current count of
- // primes, and to find a place in the current list of primes for
- // the new result.
- //
- // atomic_fetch_add results the _current_ count, but increments
- // it (by 1 in this case). The current count of primes indexes
- // into the first unoccupied position of the 'result' array.
- const int idx = Kokkos::atomic_fetch_add (&count(), 1);
- result(idx) = number;
+ for ( auto k = 0; k < local_indices.size() ; ++k)
+ outputfile << global_index << " " << A->getColMap()->getGlobalElement(local_indices[k]) << " " << (scalar) local_values[k] << std::endl;
}
- }
+ outputfile.close();
-};
+}
-using namespace mrstlnos;
+template<typename scalar> void writeVectorToFile(std::string &name, size_t myRank, Teuchos::RCP<Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> > b)
+{
+ using Teuchos::Array;
+ using Teuchos::ArrayView;
+ using Teuchos::ArrayRCP;
+ using Teuchos::arcp;
+ using Teuchos::RCP;
+ using Teuchos::rcp;
+ using Teuchos::tuple;
+
+ std::ofstream outputfile;
+ if (myRank ==0)
+ outputfile.open (name);
+
+ typedef std::numeric_limits< double > dbl;
+ outputfile.precision(dbl::max_digits10+2);
+
+ if (myRank ==0)
+ {
+ b->template sync<Kokkos::HostSpace> ();
+ auto b_2d = b->template getLocalView<Kokkos::HostSpace> ();
+ auto B = Kokkos::subview (b_2d, Kokkos::ALL (), 0);
+
+ const int ensemble_size = Kokkos::dimension_scalar(B) == 0 ? 1 : Kokkos::dimension_scalar(B);
+ typedef decltype(B) B_View;
+ typename Kokkos::FlatArrayType<B_View>::type B_flat = B;
+
+ auto line_size = b->getMap()->getGlobalNumElements ();
+ for (global_ordinal_type global_index = 0; global_index < line_size; ++global_index)
+ {
+ for(auto j=0; j<ensemble_size; ++j)
+ outputfile << B_flat[ensemble_size*global_index+j] << " ";
+ outputfile << std::endl;
+ }
+ outputfile.close();
+ }
+}
Example2::Example2()
{
}
-int
-Example2::execute(int numThreads)
+template<typename scalar> void twoGridCycle(int nu1, int nu2, Teuchos::RCP<Tpetra::CrsMatrix<scalar,local_ordinal_type,global_ordinal_type> > A, Teuchos::RCP<Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> > b, Teuchos::RCP<Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> > x, std::string &name, size_t myRank, std::string &name2, int nx,Teuchos::RCP<const Teuchos::Comm<int> > comm)
{
- // You must call initialize() before you may call Kokkos.
- Timers tms;
-
- std::ostream& out = std::cout;
-
- std::cout << "Num of threads: " << numThreads << std::endl;
-
- Kokkos::InitArguments args;
- args.num_threads = numThreads;
- Kokkos::initialize (args);
+ using Teuchos::Array;
+ using Teuchos::ArrayView;
+ using Teuchos::ArrayRCP;
+ using Teuchos::arcp;
+ using Teuchos::RCP;
+ using Teuchos::rcp;
+ using Teuchos::tuple;
+ using std::cerr;
+ using std::cout;
+ using std::endl;
+
+ using Kokkos::subview;
+ using Kokkos::ALL;
+
+ typedef Tpetra::Map<> map_type;
+ typedef Tpetra::Vector<>::scalar_type scalar_type;
+ typedef Tpetra::Vector<>::mag_type magnitude_type;
+ typedef Tpetra::MultiVector<scalar,local_ordinal_type,global_ordinal_type> multivector_type;
+ typedef Tpetra::CrsMatrix<scalar,local_ordinal_type,global_ordinal_type> crs_matrix_type;
+ typedef Tpetra::CrsGraph<local_ordinal_type,global_ordinal_type> crs_graph_type;
+
+ typedef Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> vector_type;
+ typedef Tpetra::Export<> export_type;
+ typedef typename vector_type::mag_type magnitude_type;
+ using Teuchos::RCP;
+ using Teuchos::rcp;
+ std::ofstream outputfile;
+ if (myRank ==0)
+ outputfile.open (name);
+
+ int indexBase = 0;
+ int nx_coarsed = (nx-1)/2;
+ int N_coarsed = nx_coarsed*nx_coarsed;
+
+ RCP<const map_type> map_coarsed_y_lines = rcp (new map_type (nx_coarsed, indexBase, comm));
+
+ const size_t num_my_coarsed_y_lines = map_coarsed_y_lines->getNodeNumElements ();
+ Teuchos::ArrayView<const global_ordinal_type> my_coarsed_y_lines = map_coarsed_y_lines->getNodeElementList ();
+ Array<global_ordinal_type> my_dofs (num_my_coarsed_y_lines*nx_coarsed);
+ for (auto i = 0; i < num_my_coarsed_y_lines; ++i)
+ for (auto j = 0; j < nx_coarsed; ++j)
+ my_dofs[nx_coarsed*i+j] = nx_coarsed*my_coarsed_y_lines[i]+j;
+
+ RCP<const map_type> map_coarsed_dofs = rcp (new map_type (N_coarsed, my_dofs, indexBase, comm));
+ RCP<const map_type> map_dofs = A->getRowMap();
+ // Prolongation graph
+ RCP<crs_graph_type> graph_P (new crs_graph_type (map_dofs,map_coarsed_dofs, 0));
+
+ for (auto i = 0; i < nx; ++i)
+ for (auto j = 0; j < nx_coarsed; ++j)
+ for (auto k = 0; k < nx; ++k)
+ for (auto l = 0; l < nx_coarsed; ++l)
+ {
+ global_ordinal_type global_index_1 = i*nx+k;
+ global_ordinal_type global_index_2 = j*nx_coarsed +l;
+ if(i==j*2 || i==j*2+1 || i==j*2+2)
+ if(k==l*2 || k==l*2+1 || k==l*2+2)
+ graph_P->insertGlobalIndices(global_index_1, tuple<global_ordinal_type> (global_index_2));
+ }
- //Kokkos::Threads::print_configuration(out);
+ graph_P->fillComplete ();
+ // Prolongation matrix
+ RCP<crs_matrix_type> P (new crs_matrix_type (graph_P));
+ for (auto i = 0; i < nx; ++i)
+ for (auto j = 0; j < nx_coarsed; ++j)
+ for (auto k = 0; k < nx; ++k)
+ for (auto l = 0; l < nx_coarsed; ++l)
+ {
+ global_ordinal_type global_index_1 = i*nx+k;
+ global_ordinal_type global_index_2 = j*nx_coarsed +l;
+ if(i==j*2 || i==j*2+2)
+ {
+ if(k==l*2 || k==l*2+2)
+ P->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 0.25));
+ else if (k==l*2+1)
+ P->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 0.5));
+ }
+ else if(i==j*2+1)
+ {
+ if(k==l*2 || k==l*2+2)
+ P->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 0.5));
+ else if (k==l*2+1)
+ P->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 1.));
+ }
+ }
+
+ P->fillComplete ();
+
+ std::string name_P = "prolongation_matrix.txt";
+ writeLocalCrsMatrixToFile<scalar>(name_P, myRank, P);
+
+ // Restriction graph
+ RCP<crs_graph_type> graph_R (new crs_graph_type (map_coarsed_dofs,map_dofs, 0));
+
+ for (auto i = 0; i < nx_coarsed; ++i)
+ for (auto j = 0; j < nx; ++j)
+ for (auto k = 0; k < nx_coarsed; ++k)
+ for (auto l = 0; l < nx; ++l)
+ {
+ global_ordinal_type global_index_1 = i*nx_coarsed +k;
+ global_ordinal_type global_index_2 = j*nx +l;
+ if(j==i*2 || j==i*2+1 || j==i*2+2)
+ if(l==k*2 || l==k*2+1 || l==k*2+2)
+ graph_R->insertGlobalIndices(global_index_1, tuple<global_ordinal_type> (global_index_2));
+ }
- //std::string kokkos_nbthreads_input = "--kokkos-threads="+ std::to_string(numThreads);
+ graph_R->fillComplete ();
+ // Restrictionmatrix
+ RCP<crs_matrix_type> R (new crs_matrix_type (graph_R));
+ for (auto i = 0; i < nx_coarsed; ++i)
+ for (auto j = 0; j < nx; ++j)
+ for (auto k = 0; k < nx_coarsed; ++k)
+ for (auto l = 0; l < nx; ++l)
+ {
+ global_ordinal_type global_index_1 = i*nx_coarsed+k;
+ global_ordinal_type global_index_2 = j*nx +l;
+ if(j==i*2 || j==i*2+2)
+ {
+ if(l==k*2 || l==k*2+2)
+ R->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 0.25/4.));
+ else if (l==k*2+1)
+ R->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 0.5/4.));
+ }
+ else if(j==i*2+1)
+ {
+ if(l==k*2 || l==k*2+2)
+ R->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 0.5/4.));
+ else if (l==k*2+1)
+ R->replaceGlobalValues(global_index_1, tuple<global_ordinal_type> (global_index_2), tuple<scalar> ((scalar) 1./4.));
+ }
+ }
+
+ R->fillComplete ();
- //std::cout<< kokkos_nbthreads_input << std::endl;
+ std::string name_R = "restriction_matrix.txt";
+ writeLocalCrsMatrixToFile<scalar>(name_R, myRank, R);
- //static int argc = 1;
- //char **argv = new char*[2];
- //argv[0] = new char[kokkos_nbthreads_input.length() + 1]; strcpy(argv[0], kokkos_nbthreads_input.c_str());
- //argv[1] = new char[1]; strcpy(argv[1], "");
+ //Coarsed matrix
+ RCP<crs_matrix_type> AP (new crs_matrix_type (map_dofs,map_coarsed_dofs,1));
- //Kokkos::initialize (argc, argv);
+ Tpetra::MatrixMatrix::Multiply(*A, false, *P, false, *AP, false); // AP = A*P
+ AP->fillComplete ();
+ std::cout << *AP << std::endl;
+ RCP<crs_matrix_type> AH (new crs_matrix_type (map_coarsed_dofs,map_coarsed_dofs,1));
+ Tpetra::MatrixMatrix::Multiply(*R, false, *AP, false, *AH, false); // AH = R*AP
+ AH->fillComplete ();
+ std::cout << *AH << std::endl;
+ std::string name_AH = "coarsed_matrix.txt";
+ writeLocalCrsMatrixToFile<scalar>(name_AH, myRank, AH);
- tms["total"].start();
+ const scalar one = (scalar) 1;
+ const scalar zero = (scalar) 0;
- srand (61391); // Set the random seed
+ outputfile.close();
+}
- int nnumbers = 10000000;
- view_type data ("RND", nnumbers);
- view_type result ("Prime", nnumbers);
- count_type count ("Count");
+template<typename scalar> int conjugateGradient(int niters, Teuchos::RCP<Tpetra::CrsMatrix<scalar,local_ordinal_type,global_ordinal_type> > A, Teuchos::RCP<Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> > b, Teuchos::RCP<Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> > x, std::string &name, size_t myRank, std::string &name2)
+{
+ typedef Tpetra::Map<> map_type;
+ typedef Tpetra::Vector<>::scalar_type scalar_type;
+ typedef Tpetra::Vector<>::mag_type magnitude_type;
+ typedef Tpetra::MultiVector<scalar,local_ordinal_type,global_ordinal_type> multivector_type;
+ typedef Tpetra::CrsMatrix<scalar,local_ordinal_type,global_ordinal_type> crs_matrix_type;
+ typedef Tpetra::CrsGraph<local_ordinal_type,global_ordinal_type> crs_graph_type;
+
+ typedef Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> vector_type;
+ typedef Tpetra::Export<> export_type;
+ typedef typename vector_type::mag_type magnitude_type;
+ using Teuchos::RCP;
+ using Teuchos::rcp;
+ std::ofstream outputfile;
+ if (myRank ==0)
+ outputfile.open (name);
+
+
+
+ const scalar one = (scalar) 1;
+ const scalar zero = (scalar) 0;
+ vector_type q (A->getDomainMap ());
+ vector_type z (A->getRangeMap ());
+ vector_type resid_k (A->getRangeMap ());
+ vector_type p_k (A->getRangeMap ());
+
+ magnitude_type tolerance = (magnitude_type) 0.00000000001; //0.00000001;
+
+ magnitude_type alpha_k;
+ magnitude_type beta_k;
+ magnitude_type normz;
+ magnitude_type residual_0;
+ magnitude_type residual_k;
+ magnitude_type residual_kp1;
+
+ const int reportFrequency = 1;
+ const int outputFrequency = 1;
+
+ x->scale (zero, z); //x->randomize();
+ A->apply(*x, z); // z := A * x
+ resid_k.update (one, *b, -one, z, zero); // resid_k := b - z
+ p_k.scale (one, resid_k); // p_k := resid_k
+ residual_k = resid_k.norm2();
+ residual_0 = residual_k;
+
+ std::string name3 = name2 + "_0.txt";
+ writeVectorToFile<scalar>(name3, myRank, x);
+ // Do the power method, until the method has converged or the
+ // maximum iteration count has been reached.
+ for (int iter = 1; iter <= niters; ++iter)
+ {
+ A->apply(p_k, z); // z := A * p_k
+ alpha_k = residual_k*residual_k/(p_k.dot (z)); // alpha_k := resid_k^T * resid_k /(p_k^T*A*p_k)
+ x->update(alpha_k ,p_k,one); // x := x + alpha_k*p_k
+ resid_k.update(-alpha_k,z,one);
+ residual_kp1 = resid_k.norm2();
+ // Compute and report the residual norm every reportFrequency
+ // iterations, or if we've reached the maximum iteration count.
+ if (iter % reportFrequency == 0 || iter + 1 == niters)
+ {
+ if (myRank ==0)
+ outputfile << "Iteration " << iter << ": " << residual_k/residual_0 << " " << residual_kp1/residual_0 << endl;
+ }
+ if (iter % outputFrequency == 0 || iter + 1 == niters)
+ {
+ if (myRank ==0)
+ {
+ std::string name3 = name2 + "_" + std::to_string(iter) + ".txt";
+ writeVectorToFile<scalar>(name3, myRank, x);
+ }
+
+ }
+ if (residual_kp1/residual_0 < tolerance)
+ {
+ if (myRank ==0)
+ outputfile << "Converged after " << iter << " iterations" << endl;
+ niters = iter;
+ break;
+ }
+ else if (iter+1 == niters)
+ {
+ if (myRank ==0)
+ outputfile << "Failed to converge after " << niters << " iterations" << endl;
+ break;
+ }
+ beta_k = residual_kp1*residual_kp1/(residual_k*residual_k);
+ p_k.update(one,resid_k,beta_k);
+
+ residual_k = residual_kp1;
+ }
+ outputfile.close();
+ return niters;
+}
- host_view_type h_data = Kokkos::create_mirror_view (data);
- host_view_type h_result = Kokkos::create_mirror_view (result);
- host_count_type h_count = Kokkos::create_mirror_view (count);
+template<typename scalar> int tstart(int nx, int ny, int ensemble_size, MPI_Comm _comm)
+{
+ using Teuchos::Array;
+ using Teuchos::ArrayView;
+ using Teuchos::ArrayRCP;
+ using Teuchos::arcp;
+ using Teuchos::RCP;
+ using Teuchos::rcp;
+ using Teuchos::tuple;
+ using std::cerr;
+ using std::cout;
+ using std::endl;
+
+ using Kokkos::subview;
+ using Kokkos::ALL;
+
+ typedef EnsembleTraits<scalar> ET;
+ const int real_ensemble_size = ET::size;
+
+ typedef Tpetra::Map<> map_type;
+ typedef Tpetra::Vector<>::scalar_type scalar_type;
+ typedef Tpetra::Vector<>::mag_type magnitude_type;
+ typedef Tpetra::MultiVector<scalar,local_ordinal_type,global_ordinal_type> multivector_type;
+ typedef Tpetra::CrsMatrix<scalar,local_ordinal_type,global_ordinal_type> crs_matrix_type;
+ typedef Tpetra::CrsGraph<local_ordinal_type,global_ordinal_type> crs_graph_type;
+
+ typedef Tpetra::Vector<scalar,local_ordinal_type,global_ordinal_type> vector_type;
+ typedef Tpetra::Export<> export_type;
+ typedef Amesos2::Solver<crs_matrix_type,multivector_type> solver_type;
+ typedef Kokkos::CrsMatrix<scalar,local_ordinal_type,Tpetra::CrsMatrix<>::execution_space> local_matrix_type;
+
+ // MPI comm
+ RCP<const Teuchos::Comm<int> > comm = rcp (new Teuchos::MpiComm<int> (_comm));
+ const size_t myRank = comm->getRank();
+
+ // Make an output stream (for verbose output) that only prints on
+ // Proc 0 of the communicator.
+ Teuchos::oblackholestream blackHole;
+ std::ostream& out = (myRank == 0) ? std::cout : blackHole;
+
+ // Initiate Kokkos with the given number of threads
+ Kokkos::InitArguments args;
+ args.num_threads = 1;
+ Kokkos::initialize (args);
+ {
+ int numiters = ny;
+ //
+ ny = nx;
+ double L = 1;
+ double h = L/(nx+1);
+ //
+ int N = nx*ny;
+
+ Timers tms;
+ tms["total"].start();
+
+ int indexBase = 0;
+ tms["assembling"].start();
+ // Map
+ RCP<const map_type> map_y_lines = rcp (new map_type (ny, indexBase, comm));
+ const size_t num_my_y_lines = map_y_lines->getNodeNumElements ();
+ Teuchos::ArrayView<const global_ordinal_type> my_y_lines = map_y_lines->getNodeElementList ();
+
+ std::vector<global_ordinal_type> my_nodes ={};
+ Array<global_ordinal_type> my_dofs (num_my_y_lines*nx);
+ for (auto i = 0; i < num_my_y_lines; ++i)
+ for (auto j = 0; j < nx; ++j)
+ my_dofs[nx*i+j] = nx*my_y_lines[i]+j;
+
+ RCP<const map_type> map_dofs = rcp (new map_type (N, my_dofs, indexBase, comm));
+
+ const size_t num_lcl_indices = (myRank == 0) ? N : 0;
+ RCP<const Tpetra::Map<>> map_out = rcp (new Tpetra::Map<> (N, num_lcl_indices, indexBase, comm));
+
+ // Graph
+ RCP<crs_graph_type> graph_A (new crs_graph_type (map_dofs, 0));
+
+ for (auto i = 0; i < num_my_y_lines; ++i)
+ for (auto j = 0; j < nx; ++j)
+ {
+ global_ordinal_type global_index = map_dofs->getGlobalElement(nx*i+j);
+ graph_A->insertGlobalIndices(global_index, tuple<global_ordinal_type> (global_index));
+ if(global_index-nx >= 0 )
+ graph_A->insertGlobalIndices(global_index, tuple<global_ordinal_type> (global_index-nx));
+ if(global_index+nx < N )
+ graph_A->insertGlobalIndices(global_index, tuple<global_ordinal_type> (global_index+nx));
+ if(global_index-1 >= 0 && j != 0)
+ graph_A->insertGlobalIndices(global_index, tuple<global_ordinal_type> (global_index-1));
+ if(global_index+1 < N && j != nx-1)
+ graph_A->insertGlobalIndices(global_index, tuple<global_ordinal_type> (global_index+1));
+ }
- typedef view_type::size_type size_type;
- // Fill the 'data' array on the host with random numbers. We assume
- // that they come from some process which is only implemented on the
- // host, via some library. (That's true in this case.)
- for (size_type i = 0; i < data.dimension_0 (); ++i) {
- h_data(i) = rand () % nnumbers;
+ graph_A->fillComplete ();
+
+ // Matrix
+ RCP<crs_matrix_type> A (new crs_matrix_type (graph_A));
+ for (auto i = 0; i < num_my_y_lines; ++i)
+ for (auto j = 0; j < nx; ++j)
+ {
+ global_ordinal_type global_index = map_dofs->getGlobalElement(nx*i+j);
+ A->replaceGlobalValues(global_index, tuple<global_ordinal_type> (global_index), tuple<scalar> ((scalar) 4));
+ if(global_index-nx >= 0 )
+ A->replaceGlobalValues(global_index, tuple<global_ordinal_type> (global_index-nx), tuple<scalar> ((scalar) -1));
+ if(global_index+nx < N )
+ A->replaceGlobalValues(global_index, tuple<global_ordinal_type> (global_index+nx), tuple<scalar> ((scalar) -1));
+ if(global_index-1 >= 0 && j != 0)
+ A->replaceGlobalValues(global_index, tuple<global_ordinal_type> (global_index-1), tuple<scalar> ((scalar) -1));
+ if(global_index+1 < N && j != nx-1)
+ A->replaceGlobalValues(global_index, tuple<global_ordinal_type> (global_index+1), tuple<scalar> ((scalar) -1));
+ }
+ A->fillComplete ();
+
+ // b
+ int x_mode = 3;
+ int y_mode = 2;
+
+ RCP<vector_type> b (new vector_type(map_dofs, true));
+ RCP<vector_type> b_out (new vector_type(map_out, true));
+ for (auto i = 0; i < num_my_y_lines; ++i)
+ for (auto j = 0; j < nx; ++j)
+ {
+ scalar value;
+ if (real_ensemble_size >= 1)
+ for (int sample=0; sample<ensemble_size; ++sample)
+ {
+ ET::coeff(value, sample) = sin((i+1)*(y_mode+sample)*M_PI/(ny+1))* sin((j+1)*(x_mode+sample)*M_PI/(nx+1))/(h*h);
+ }
+ else
+ value = sin((i+1)*y_mode*M_PI/(ny+1))* sin((j+1)*x_mode*M_PI/(nx+1))/(h*h);
+ global_ordinal_type global_index = map_dofs->getGlobalElement(nx*i+j);
+ b->replaceGlobalValue(global_index, value);
+ }
+ RCP<vector_type> x (new vector_type(map_dofs, true));
+ RCP<vector_type> x_out (new vector_type(map_out, true));
+ tms["assembling"].stop();
+
+ tms["KLU2"].start();
+ RCP<solver_type> solverA;
+ solverA = Amesos2::create<crs_matrix_type,multivector_type>("KLU2", A, x, b);
+ solverA->solve();
+ tms["KLU2"].stop();
+
+
+ const export_type export_out(map_dofs,map_out);
+ b_out->doExport(*b, export_out, Tpetra::INSERT);
+ x_out->doExport(*x, export_out, Tpetra::INSERT);
+
+ std::string name_A = "A_E" + std::to_string(ensemble_size);
+ std::string name_b = "b_E" + std::to_string(ensemble_size) + ".txt";
+ std::string name_x = "x_KLU2_E" + std::to_string(ensemble_size) + ".txt";
+ std::string name_x2 = "x_CG_E" + std::to_string(ensemble_size) + ".txt";
+ std::string name_CG = "CG_E" + std::to_string(ensemble_size) + ".txt";
+ std::string name_CG2 = "CG_E" + std::to_string(ensemble_size) + "_iter";
+ std::string name_timer = "timers_E" + std::to_string(ensemble_size) + ".txt";
+ writeLocalCrsMatrixToFile<scalar>(name_A, myRank, A);
+
+ writeVectorToFile<scalar>(name_b, myRank, b);
+ writeVectorToFile<scalar>(name_x, myRank, x);
+ tms["CG"].start();
+ conjugateGradient<scalar> (numiters,A,b,x,name_CG,myRank,name_CG2);
+ tms["CG"].stop();
+ writeVectorToFile<scalar>(name_x2, myRank, x);
+ tms["total"].stop();
+
+ scalar one = Teuchos::ScalarTraits<scalar>::one();
+ twoGridCycle(5,5,A,b,x,name_CG,myRank,name_CG2,nx,comm);
+ //Tpetra::MatrixMatrix::Add(*A,false,one, *A, one);
+
+ std::ofstream outputfile;
+ if (myRank ==0)
+ outputfile.open (name_timer);
+
+ outputfile << "\n---CPU statistics of proc 0 for ensemble of size: "<< ensemble_size<< " ---\n";
+ outputfile << tms;
+ outputfile << "---------------------------\n";
+ outputfile.close();
}
- Kokkos::deep_copy (data, h_data); // copy from host to device
-
- Kokkos::parallel_for (data.dimension_0 (), findprimes (data, result, count));
- Kokkos::deep_copy (h_count, count); // copy from device to host
-
- printf ("Found %i prime numbers in %i random numbers\n", h_count(), nnumbers);
Kokkos::finalize ();
+ comm = Teuchos::null;
+ return 0;
+
+}
- tms["total"].stop();
- std::cout << "\n---CPU statistics of proc ---\n";
- std::cout << tms;
- std::cout << "---------------------------\n";
+int Example2::execute(int nx, int ny, int ensemble_size, MPI_Comm _comm)
+{
+ if (ensemble_size > 0)
+ {
+ if (ensemble_size <= 1)
+ {
+ typedef Sacado::MP::Vector<Stokhos::StaticFixedStorage<int,double,1,Kokkos::Threads>> ensemble_scalar;
+ int out = tstart<ensemble_scalar>(nx,ny,ensemble_size,_comm);
+ return out;
+ }
+ else if (ensemble_size <= 2)
+ {
+ typedef Sacado::MP::Vector<Stokhos::StaticFixedStorage<int,double,2,Kokkos::Threads>> ensemble_scalar;
+ int out = tstart<ensemble_scalar>(nx,ny,ensemble_size,_comm);
+ return out;
+ }
+ else if (ensemble_size <= 4)
+ {
+ typedef Sacado::MP::Vector<Stokhos::StaticFixedStorage<int,double,4,Kokkos::Threads>> ensemble_scalar;
+ int out = tstart<ensemble_scalar>(nx,ny,ensemble_size,_comm);
+ return out;
+ }
+ else if (ensemble_size <= 6)
+ {
+ typedef Sacado::MP::Vector<Stokhos::StaticFixedStorage<int,double,6,Kokkos::Threads>> ensemble_scalar;
+ int out = tstart<ensemble_scalar>(nx,ny,ensemble_size,_comm);
+ return out;
+ }
+ else if (ensemble_size <= 8)
+ {
+ typedef Sacado::MP::Vector<Stokhos::StaticFixedStorage<int,double,8,Kokkos::Threads>> ensemble_scalar;
+ int out = tstart<ensemble_scalar>(nx,ny,ensemble_size,_comm);
+ return out;
+ }
+ }
+ else
+ {
+ typedef double scalar;
+ int out = tstart<scalar>(nx,ny,ensemble_size,_comm);
+ return out;
+ }
}
diff --git a/mrstlnos/src/wExample2.h b/mrstlnos/src/wExample2.h
index a4e0ad2c..675afc15 100644
--- a/mrstlnos/src/wExample2.h
+++ b/mrstlnos/src/wExample2.h
@@ -6,10 +6,8 @@
#include <iostream>
#include <vector>
#include <string>
-#ifndef SWIG
-#include "Teuchos_Comm.hpp"
-#include "Teuchos_XMLParameterListHelpers.hpp"
-#endif
+
+#include <mpi.h>
namespace mrstlnos {
@@ -21,8 +19,8 @@ class MRSTLNOS_API Example2 : public wObject
{
public:
Example2();
- int execute(int numThreads);
-
+ int execute(int nx, int ny, int ensemble_size, MPI_Comm _comm);
+
};
}
diff --git a/mrstlnos/tests/bolt_iterative.py b/mrstlnos/tests/bolt_iterative.py
new file mode 100644
index 00000000..d923ef62
--- /dev/null
+++ b/mrstlnos/tests/bolt_iterative.py
@@ -0,0 +1,193 @@
+#!/usr/bin/env python
+# -*- coding: latin-1; -*-
+
+import mrstlnos as m
+import tbox
+import tbox.gmsh as gmsh
+
+def main():
+ try:
+ import mpi4py.MPI as mpi
+ comm = mpi.COMM_WORLD
+ rank = comm.rank
+ siz = comm.size
+ name = mpi.Get_processor_name()
+ status = mpi.Status()
+ print "info: MPI found"
+ except:
+ comm = None
+ rank = 0
+ siz = 1
+ name = "noname"
+ print "info: MPI not found => MPI disabled"
+
+
+
+ from fwk.wutils import parseargs
+ args=parseargs()
+
+ from PyTrilinos import Teuchos
+
+ geo_name = 'bolt.geo'
+ mesh_name = 'bolt.msh'
+
+
+ import shlex, subprocess
+ import os
+ import numpy as np
+
+ file_dir = os.path.dirname(__file__)
+ work_dir = os.getcwd()
+
+ if rank == 0:
+ command_line = 'gmsh -3 ' + file_dir + '/' + geo_name + ' -o ' + work_dir + '/' + mesh_name +' -part ' + str(siz)
+ print command_line
+ tmp = shlex.split(command_line)
+ fileout = open('gmsh.log','w')
+ p = subprocess.Popen(tmp, stdin=subprocess.PIPE, stdout=fileout, stderr=fileout, env=os.environ, shell=False,close_fds=True)
+ retcode = p.wait()
+ for i in range(1,siz):
+ comm.send(1, dest=i, tag=11)
+ else:
+ data = comm.recv(source=0, tag=11)
+
+
+ if siz == 1:
+ msh = gmsh.MeshLoader(mesh_name,work_dir).execute()
+ else:
+ msh = gmsh.MeshLoader(mesh_name,work_dir).myrank_execute(rank)
+
+ pbl = m.Problem(msh,comm)
+
+ ensemble_size = 1
+ E = 330000 * np.ones(ensemble_size)
+ nu = 0.36 * np.ones(ensemble_size)
+ k = 0.114
+ d = 0.5E-05
+
+ m.Medium(pbl, "Holder 1", "TZM",E,nu,k,d)
+
+ E = 200000 * np.ones(ensemble_size)
+ nu = 0.36 * np.ones(ensemble_size)
+ k = 0.01
+ d = 1.75E-05
+ m.Medium(pbl, "Holder 2", "SS316LN",E,nu,k,d)
+ m.Medium(pbl, "Bolt", "Inconel 718",E,nu,k,d)
+ E = 125000 * np.ones(ensemble_size)
+ nu = 0.36 * np.ones(ensemble_size)
+ k = 0.325
+ d = 1.75E-05
+ m.Medium(pbl, "Spacer", "CuCrZr",E,nu,k,d)
+
+ m.uDirichlet(pbl, "Holder 1 end", "Heated",1,0.,0,0.,0,0.,1,300.,1)
+ m.uDirichlet(pbl, "Holder 2 end", "Clamped",1,0.,1,0.,1,0.,1,70.,1)
+ m.uDirichlet(pbl, "Sym y","y symetry bc",0,0.,1,0.,0,0.,0,0.,1)
+ m.uDirichlet(pbl, "Bolt center","sym",1,0.,1,0.,0,0.,0,0.,1)
+
+ norm1 = tbox.Pt(0,0,-1)
+ norm2 = tbox.Pt(0,0,1)
+
+ cont1 = m.Contact(pbl,"Bolt to Holder 1","contact",norm1,2)
+ cont1.setMaster(pbl,"Holder 1 to Bolt",norm1)
+
+ cont2 = m.Contact(pbl,"Bolt to Holder 2","contact",norm2,2)
+ cont2.setMaster(pbl,"Holder 2 to Bolt",norm2)
+
+ cont3 = m.Contact(pbl,"Spacer to Holder 1","contact",norm2)
+ cont3.setMaster(pbl,"Holder 1 to Spacer",norm2)
+
+ cont4 = m.Contact(pbl,"Spacer to Holder 2","contact",norm1)
+ cont4.setMaster(pbl,"Holder 2 to Spacer",norm1)
+
+ precondList = Teuchos.ParameterList()
+ precondList = Teuchos.ParameterList()
+ precondSubList = Teuchos.ParameterList()
+ precondSubSubList = Teuchos.ParameterList()
+ precondSubSubList['fact: iluk level-of-fill'] = 2
+ precondSubList['Prec Type'] = "RILUK"
+ precondSubList['Overlap'] = 1
+ precondSubList['Ifpack2 Settings'] = precondSubSubList
+ precondList['Ifpack2'] = precondSubList
+
+ mueluParams = Teuchos.ParameterList()
+ #mueluParamsSub1 = Teuchos.ParameterList()
+ #mueluParamsSub1['PDE equations'] = 4
+ #mueluParams['Matrix'] = mueluParamsSub1
+ mueluParams['number of equations'] = 3
+ mueluParams['max levels'] = 5
+ mueluParams['multigrid algorithm'] = "sa"
+ mueluParams['problem: symmetric'] = True
+ mueluParams['coarse: max size'] = 10
+ #mueluParams['aggregation: type'] = "uncoupled"
+ #mueluParams['aggregation: min agg size'] = 3
+ #mueluParams['aggregation: max agg size'] = 10
+ mueluParams['aggregation: export visualization data'] = True
+ mueluParams['smoother: type'] = "RELAXATION"
+
+ mueluParamsSub2 = Teuchos.ParameterList()
+
+ mueluParamsSub2['relaxation: type'] = "Symmetric Gauss-Seidel"
+ mueluParamsSub2['relaxation: sweeps'] = 5
+ mueluParamsSub2['relaxation: damping factor'] = 0.9
+
+ mueluParams['coarse: type'] = "Klu"
+ mueluParams['coarse: max size'] = 18
+ mueluParams['verbosity'] = "low"
+
+ mueluParamsSub3 = Teuchos.ParameterList()
+
+ mueluParamsSub4 = Teuchos.ParameterList()
+ mueluParamsSub4['A'] = '{0,1,2}'
+ mueluParamsSub4['P'] = '{0,1,2}'
+ mueluParamsSub4['R'] = '{0,1}'
+
+ mueluParams['smoother: params'] = mueluParamsSub2
+ mueluParams['export data'] = mueluParamsSub4
+
+
+ solverList = Teuchos.ParameterList()
+ solverList['type']= "BlockGmres"
+ solverList['Convergence Tolerance'] = 10**(-8)
+ solverList['Maximum Iterations'] = 100
+ solverList['Verbosity'] = 33
+ solverList['Flexible Gmres'] = False
+ solverList['Output Style'] = 1
+ solverList['Output Frequency'] = 1
+ solverList['Output Stream file name belos'] = 'multiFGuess_belos_preconditioned.text'
+ #solverList['Output Stream file name muelu'] = 'muelu_prec.text'
+ solverList['Use preconditioner'] = True
+ solverList['ifpackParams'] = precondList
+ solverList['mueluParams'] = mueluParams
+ solverList['Initial guess'] = 'multiFGuess.text'
+ slv = m.IterativeSolver(pbl,args.k,solverList,3, False, ensemble_size)
+
+ if not args.nogui:
+ if rank == 0:
+ all_msh = gmsh.MeshLoader(mesh_name,work_dir).execute()
+ import mrstlnos.viewer as v
+ gui = v.MeshViewer(pbl,slv,all_msh)
+ gui.vmin = 0
+ gui.vmax = 1
+ gui.start()
+ else:
+ slv.start()
+ else:
+ slv.start()
+
+ timers = slv.tms;
+ if rank == 0:
+ print 'Map real CPU cost'
+ print timers['map'].read().getReal()
+ print 'Graph real CPU cost'
+ print timers['graph'].read().getReal()
+ print 'Matrix computation on each mpi process real CPU cost'
+ print timers['computeMatrices'].read().getReal()
+ print 'MPI sharing CPU cost'
+ print timers['assembly-MPI'].read().getReal()
+ print 'prec'
+ print timers['prec'].read().getReal()
+ print 'belos'
+ print timers['belos'].read().getReal()
+
+if __name__ == "__main__":
+ main()
diff --git a/mrstlnos/tests/cube.geo b/mrstlnos/tests/cube.geo
index 408f84cb..62728c32 100644
--- a/mrstlnos/tests/cube.geo
+++ b/mrstlnos/tests/cube.geo
@@ -2,7 +2,7 @@
// fichier de donnees gmsh
C = 10;
-nC = 8;
+nC = 16;
lc = 1;
diff --git a/mrstlnos/tests/cube_new_iterative.py b/mrstlnos/tests/cube_new_iterative.py
index 40517474..479efb46 100644
--- a/mrstlnos/tests/cube_new_iterative.py
+++ b/mrstlnos/tests/cube_new_iterative.py
@@ -59,20 +59,20 @@ def main():
pbl = m.Problem(msh,comm)
- ensemble_size = 2
+ ensemble_size = 1
E = 1. * np.ones(ensemble_size)
nu = 0.2 * np.ones(ensemble_size)
m.Medium(pbl, "Body", "test",E,nu)
- m.uDirichlet(pbl, "Surf 1", "Clamped",1,0.,1,0.,1,0.,0,0.,ensemble_size)
+ m.uDirichlet(pbl, "Surf 1", "Clamped",0,0.,1,0.,1,0.,0,0.,ensemble_size)
dx = np.zeros(ensemble_size)
dy = 0. * np.ones(ensemble_size)
dz = np.zeros(ensemble_size)
dT = np.zeros(ensemble_size)
- m.uDirichlet(pbl, "Surf 2", "Moved",0,dx,1,dy,0,dz,0,dT)
+ m.uDirichlet(pbl, "Surf 2", "Moved",1,dx,0,dy,1,dz,1,dT)
precondList = Teuchos.ParameterList()
precondList = Teuchos.ParameterList()
@@ -89,45 +89,52 @@ def main():
#mueluParamsSub1['PDE equations'] = 4
#mueluParams['Matrix'] = mueluParamsSub1
mueluParams['number of equations'] = 3
- mueluParams['max levels'] = 2
- mueluParams['multigrid algorithm'] = "unsmoothed"
+ mueluParams['max levels'] = 5
+ mueluParams['multigrid algorithm'] = "sa"
+ mueluParams['problem: symmetric'] = True
+ mueluParams['coarse: max size'] = 10
+ #mueluParams['aggregation: type'] = "uncoupled"
+ #mueluParams['aggregation: min agg size'] = 3
+ #mueluParams['aggregation: max agg size'] = 10
+ mueluParams['aggregation: export visualization data'] = True
mueluParams['smoother: type'] = "RELAXATION"
mueluParamsSub2 = Teuchos.ParameterList()
mueluParamsSub2['relaxation: type'] = "Symmetric Gauss-Seidel"
- mueluParamsSub2['relaxation: sweeps'] = 4
- mueluParamsSub2['relaxation: damping factor'] = 1.25
-
- mueluParams['coarse: type'] = "RELAXATION"
+ mueluParamsSub2['relaxation: sweeps'] = 5
+ mueluParamsSub2['relaxation: damping factor'] = 0.9
+ mueluParams['coarse: type'] = "Klu"
+ mueluParams['coarse: max size'] = 18
mueluParams['verbosity'] = "low"
mueluParamsSub3 = Teuchos.ParameterList()
- mueluParamsSub3['relaxation: type'] = "Symmetric Gauss-Seidel"
- mueluParamsSub3['relaxation: sweeps'] = 4
- mueluParamsSub3['relaxation: damping factor'] = 1.25
-
+ mueluParamsSub4 = Teuchos.ParameterList()
+ mueluParamsSub4['A'] = '{0,1,2}'
+ mueluParamsSub4['P'] = '{0,1,2}'
+ mueluParamsSub4['R'] = '{0,1}'
mueluParams['smoother: params'] = mueluParamsSub2
- mueluParams['coarse: params'] = mueluParamsSub3
+ mueluParams['export data'] = mueluParamsSub4
+
solverList = Teuchos.ParameterList()
solverList['type']= "BlockGmres"
solverList['Convergence Tolerance'] = 10**(-8)
- solverList['Maximum Iterations'] = 1
+ solverList['Maximum Iterations'] = 100
solverList['Verbosity'] = 33
solverList['Flexible Gmres'] = False
solverList['Output Style'] = 1
solverList['Output Frequency'] = 1
- solverList['Output Stream file name belos'] = 'multi_frequence_ensemble1_belos_1it.text'
+ solverList['Output Stream file name belos'] = 'multiFGuess_belos_preconditioned.text'
#solverList['Output Stream file name muelu'] = 'muelu_prec.text'
- solverList['Use preconditioner'] = False
+ solverList['Use preconditioner'] = True
solverList['ifpackParams'] = precondList
solverList['mueluParams'] = mueluParams
- solverList['Initial guess'] = 'multi_frequence_ensemble1_guess.text'
- slv = m.IterativeSolver(pbl,args.k,solverList,3, True, ensemble_size)
+ solverList['Initial guess'] = 'multiFGuess.text'
+ slv = m.IterativeSolver(pbl,args.k,solverList,3, False, ensemble_size)
if not args.nogui:
if rank == 0:
diff --git a/mrstlnos/tests/example2.py b/mrstlnos/tests/example2.py
new file mode 100644
index 00000000..a21e3dee
--- /dev/null
+++ b/mrstlnos/tests/example2.py
@@ -0,0 +1,32 @@
+#!/usr/bin/env python
+# -*- coding: latin-1; -*-
+# runs the basic demo of Trilinos
+# trilinos-12.6.1-Source/demos/buildAgainstTrilinos
+
+import fwk.wutils as wu
+import mrstlnos
+
+def main(**d):
+ try:
+ import mpi4py.MPI as mpi
+ comm = mpi.COMM_WORLD
+ rank = comm.rank
+ siz = comm.size
+ name = mpi.Get_processor_name()
+ status = mpi.Status()
+ print "info: MPI found"
+ except:
+ comm = None
+ rank = 0
+ siz = 1
+ name = "noname"
+ print "info: MPI not found => MPI disabled"
+
+ nx = 99 # = ny
+ num_CG_iters = 200
+ expl = mrstlnos.Example2()
+ expl.execute(nx,num_CG_iters,1,comm)
+
+if __name__ == "__main__":
+ main()
+
--
GitLab