From bb2f418ba9da6facaf62c151c557e885712837fd Mon Sep 17 00:00:00 2001
From: ErinKuci <Erin.Kuci@ulg.ac.be>
Date: Wed, 31 Oct 2018 21:44:33 +0100
Subject: [PATCH] check sesnitivity with fd; adapt conveks or the unconstrained
case
---
Ccore/onelab_optimize.geo | 8 +++---
Ccore/onelab_optimize.py | 56 ++++++++++++++++++---------------------
Ccore/shape.py | 2 +-
Ccore/shp.plt | 53 ------------------------------------
Lbracket/topo.py | 2 +-
Team25/shape.pro | 2 ++
Team25/shape.py | 2 +-
7 files changed, 35 insertions(+), 90 deletions(-)
delete mode 100644 Ccore/shp.plt
diff --git a/Ccore/onelab_optimize.geo b/Ccore/onelab_optimize.geo
index d4ce25f..49a767b 100644
--- a/Ccore/onelab_optimize.geo
+++ b/Ccore/onelab_optimize.geo
@@ -36,13 +36,13 @@ If(PerturbMesh == 1)
Plugin(ModifyComponents).Expression1 = "y";
Plugin(ModifyComponents).Expression2 = "z";
Plugin(ModifyComponents).Run;
- Save View[0] StrCat(modelpath, Sprintf("view_%g.msh",0));
+ //Save View[0] StrCat(modelpath, Sprintf("view_%g.msh",0));
// relocate the vertices of the original mesh on the perturbed geometry
-
RelocateMesh Point "*";
RelocateMesh Line "*";
RelocateMesh Surface "*";
+ Save StrCat(modelpath, modelName, Sprintf("Perturb_%g.msh",VelocityTag));
// Create a view with the perturbed node coordinates as vector dataset
Plugin(NewView).NumComp = 3;
@@ -52,7 +52,7 @@ If(PerturbMesh == 1)
Plugin(ModifyComponents).Expression1 = "y";
Plugin(ModifyComponents).Expression2 = "z";
Plugin(ModifyComponents).Run;
- Save View[1] StrCat(modelpath, Sprintf("view_%g.msh",1));
+ //Save View[1] StrCat(modelpath, Sprintf("view_%g.msh",1));
// compute the velocity field by subtracting the two vector datasets
Plugin(ModifyComponents).View = 0;
@@ -62,7 +62,7 @@ If(PerturbMesh == 1)
Plugin(ModifyComponents).Expression2 = Sprintf("(w2 - v2)/(%g)", Perturbation);
Plugin(ModifyComponents).Run;
View.Name = "velocity";
- Save View[0] StrCat(modelpath, Sprintf("view_%g.msh",2));
+ //Save View[0] StrCat(modelpath, Sprintf("view_%g.msh",2));
Delete View[1];
Save View[0] StrCat(modelpath, Sprintf("velocity_%g.msh",VelocityTag));
SendToServer View[0] Sprintf("Optimization/Results/velocity_%g",VelocityTag);
diff --git a/Ccore/onelab_optimize.py b/Ccore/onelab_optimize.py
index e2a5124..0ec1b8f 100644
--- a/Ccore/onelab_optimize.py
+++ b/Ccore/onelab_optimize.py
@@ -2,6 +2,7 @@ import numpy as np
import conveks
import onelab
import scipy.optimize as spo
+from shutil import copyfile
c = onelab.client(__file__)
x = {}
@@ -63,33 +64,38 @@ def grad(xFromOpti):
grads.append(c.getNumber('Optimization/Results/dwdtau_'+str(dv)))
return np.asarray(grads)
+def grad_fd(xFromOpti, obj, basefile, perturb):
+ # generate the perturbed meshes (for each design variable)
+ getVelocityField(x,perturb)
+
+ # run the simulation for the perturbed mesh
+ grads = []
+ for dv, var in x.iteritems():
+ copyfile(basefile+'Perturb_'+str(dv)+'.msh', basefile+'.msh')
+ c.runSubClient('myGetDP', mygetdp + '-solve GetPerformances')
+ obj_perturb = c.getNumber('Optimization/Results/w')
+ grads.append((obj_perturb-obj)/perturb)
+ return np.asarray(grads)
+
def test(n):
if n==0: return
+ perturb = 1e-5
tab = np.zeros((n+1,4))
for step in range(n+1):
xstep = x[0][2] + (x[0][3]-x[0][2])/float(n)*step
- f,g,d = fobj([xstep]), grad([xstep])[0], 0
- tab[step] = [xstep, f, g, d]
- if step >= 2 :
- tab[step-1][3] = (tab[step][1]-tab[step-2][1]) \
- / (tab[step][0]-tab[step-2][0])
- if c.getString(c.name+'/Action') == 'stop':
- exit(1)
- print
- for step in range(n+1):
- ind1, ind2 = max(step-1,0), min(step+1,n)
- tab[step][3] = (tab[ind2][1]-tab[ind1][1])/(tab[ind2][0]-tab[ind1][0])
+ f,g = fobj([xstep]), grad([xstep])[0]
+ g_fd = grad_fd([xstep], f, modelName, perturb)[0]
+ tab[step] = [xstep, f, g, g_fd]
print "%4d" % step,
for item in tab[step] : print "%3.6f" % item,
print
- return
+ if c.getString(c.name+'/Action') == 'stop': exit(1)
def callbackF(Xi):
print Xi
if (c.getString(c.name+'/Action') == 'stop'):
raise StopIteration
-
def Optimize(Optimizer):
print Optimizer
print x
@@ -114,22 +120,15 @@ def Optimize(Optimizer):
for label, var in x.iteritems():
x[label] = var
initialPoint[label] = var[1]
- lowerBound[label] = var[2]
- upperBound[label] = var[3]
+ lowerBound[label] = var[2]
+ upperBound[label] = var[3]
# Initialize the MMA optimizer
conveks.mma.initialize(initialPoint, lowerBound, upperBound)
# Set some options for MMA
- conveks.mma.option.setNumber('General.Verbosity', 4) # def 4
- #conveks.mma.option.setNumber('General.SaveHistory', 0)
- #conveks.mma.option.setNumber('SubProblem.isRobustMoveLimit', 1) # def 1
- #conveks.mma.option.setNumber('SubProblem.isRobustAsymptotes', 1) # def 1
- #conveks.mma.option.setNumber('SubProblem.type', 0)
- #conveks.mma.option.setNumber('SubProblem.addConvexity', 0)
- #conveks.mma.option.setNumber('SubProblem.asymptotesRmax', 100.0)
- #conveks.mma.option.setNumber('SubProblem.asymptotesRmin', 0.001) # def 0.001
-
+ conveks.mma.option.setNumber('General.Verbosity', 0) # def 4
+
# Get iteration count (here it will be 1 - could be different in case of restart)
it = conveks.mma.getOuterIteration()
change_mma = 1.0
@@ -139,9 +138,6 @@ def Optimize(Optimizer):
xFromMMA = conveks.mma.getCurrentPoint()
- constraints = np.array([np.sum(xFromMMA)/10.0-1.0])
- grad_constraints = np.ones(numVariables)/10.0
-
if (c.getString(c.name+'/Action') == 'stop'): break
obj = fobj(xFromMMA)
if (c.getString(c.name+'/Action') == 'stop'): break
@@ -161,11 +157,11 @@ def Optimize(Optimizer):
objective = obj*fscale
grad_objective = grad_obj*fscale
- c.sendInfo('Optimization: it. {}, obj. {}, constr. {}'.format(it,objective,constraints[0]))
+ c.sendInfo('Optimization: it. {}, obj. {}'.format(it,objective))
# call MMA and update the current point
- #conveks.mma.setMoveLimits(lowerBound, upperBound, 1e-4) # parametre interessant
- conveks.mma.updateCurrentPoint(constraints, grad_objective, grad_constraints)
+ #conveks.mma.setMoveLimits(lowerBound, upperBound, 0.01) # parametre interessant
+ conveks.mma.updateCurrentPoint(grad_objective)
change_mma = conveks.mma.getDesignChange()
if change_mma < maxChange: break
it = conveks.mma.getOuterIteration()
diff --git a/Ccore/shape.py b/Ccore/shape.py
index 8c75326..029dc79 100644
--- a/Ccore/shape.py
+++ b/Ccore/shape.py
@@ -212,7 +212,7 @@ elif Optimizer=='mma':
# call MMA and update the current point
#conveks.mma.setMoveLimits(lowerBound, upperBound, 1e-0) # parametre interessant
- conveks.mma.updateCurrentPoint(constraints,grad_objective,grad_constraints)
+ conveks.mma.updateCurrentPoint(grad_objective,constraints,grad_constraints)
change = conveks.mma.getDesignChange()
it = conveks.mma.getOuterIteration()
diff --git a/Ccore/shp.plt b/Ccore/shp.plt
deleted file mode 100644
index 7fee2e8..0000000
--- a/Ccore/shp.plt
+++ /dev/null
@@ -1,53 +0,0 @@
-#set terminal pdf
-
-set grid
-set ytics nomirror
-set y2tics
-set style data lines
-
-set xlabel "tau"
-set y2label "w"
-set ylabel "dw/dtau"
-
-set output "aaa.pdf"
-plot "aaa" u 2:3 axis x1y2 t"w",\
- "aaa" u 2:4 t"dwdt",\
- "aaa" u 2:5 w linesp t"dwdt FD"
-
-pause -1
-
-set output "optishape1.pdf"
-plot "Torque_linTH.txt" u 2:3 axis x1y2 t"w",\
- "Torque_linTH.txt" u 2:4 t"dwdt",\
- "Torque_linTH.txt" u 2:5 w linesp t"dwdt FD"
-
-pause -1
-
-set output "optishape2.pdf"
-plot "LieDevSmoother.txt" u 2:3 axis x1y2 t"w",\
- "LieDevSmoother.txt" u 2:4 t"dwdt",\
- "LieDevSmoother.txt" u 2:5 w linesp t"dwdt FD"
-
-pause -1
-
-set output "optishape3.pdf"
-plot "Torque_linTH2.txt" u 2:3 axis x1y2 t"w",\
- "Torque_linTH2.txt" u 2:4 t"dwdt",\
- "Torque_linTH2.txt" u 2:5 w linesp t"dwdt FD"
-
-pause -1
-
-set output "optishape4.pdf"
-plot "LieDevSmoother2.txt" u 2:3 axis x1y2 t"w",\
- "LieDevSmoother2.txt" u 2:4 t"dwdt",\
- "LieDevSmoother2.txt" u 2:5 w linesp t"dwdt FD"
-
-pause -1
-
-set output "optishape5.pdf"
-plot "Torque_TH.txt" u 2:3 axis x1y2 t"w",\
- "Torque_TH.txt" u 2:4 t"dwdt",\
- "Torque_TH.txt" u 2:5 w linesp t"dwdt FD"
-
-pause -1
-
diff --git a/Lbracket/topo.py b/Lbracket/topo.py
index 40c77b3..92dd5b2 100644
--- a/Lbracket/topo.py
+++ b/Lbracket/topo.py
@@ -155,7 +155,7 @@ while change > maxChange and it <= maxIter and c.getString('topo/Action') != 'st
# call MMA and update the current point
#conveks.mma.setMoveLimits(lowerBound, upperBound, 0.1)
- conveks.mma.updateCurrentPoint(constraints, grad_objective, grad_constraints)
+ conveks.mma.updateCurrentPoint(grad_objective, constraints, grad_constraints)
change = conveks.mma.getDesignChange()
it = conveks.mma.getOuterIteration()
diff --git a/Team25/shape.pro b/Team25/shape.pro
index 8d3983d..d27346f 100644
--- a/Team25/shape.pro
+++ b/Team25/shape.pro
@@ -283,6 +283,8 @@ Formulation {
Equation {
Galerkin { [ nu[{d a}] * Dof{d Lva}, {d Lva} ];
In Domain; Jacobian Vol; Integration Int2; }
+ Galerkin { [ dhdb_NL[{d a}] * Dof{d Lva} , {d Lva} ];
+ In Core; Jacobian Vol ; Integration Int2 ; }
Galerkin { [LieOf_HB[{d a},dV[{d v_1},{d v_2},{d v_3}]],{d Lva}];
In Domain; Jacobian Vol; Integration Int2; }
Galerkin { [LieOf_HB_NL[{d a},dV[{d v_1},{d v_2},{d v_3}]],{d Lva}];
diff --git a/Team25/shape.py b/Team25/shape.py
index 25ca547..35e0834 100644
--- a/Team25/shape.py
+++ b/Team25/shape.py
@@ -158,7 +158,7 @@ while it <= maxIter and c.getString('shape/Action') != 'stop':
# call MMA and update the current point
conveks.mma.setMoveLimits(lowerBound, upperBound, 1.0e-3)
- conveks.mma.updateCurrentPoint(constraints,grad_objective,grad_constraints)
+ conveks.mma.updateCurrentPoint(grad_objective,constraints,grad_constraints)
change = conveks.mma.getDesignChange()
it = conveks.mma.getOuterIteration()
--
GitLab