Skip to content
Snippets Groups Projects
Commit bb2f418b authored by Erin Kuci's avatar Erin Kuci
Browse files

check sesnitivity with fd; adapt conveks or the unconstrained case

parent bbfd67c2
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Ccore/onelab_optimize.geo
+ 4
− 4
View file @ bb2f418b
......@@ -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);
......
Ccore/onelab_optimize.py
+ 26
− 30
View file @ bb2f418b
......@@ -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
......@@ -121,14 +127,7 @@ def Optimize(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()
......@@ -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()
......
Ccore/shape.py
+ 1
− 1
View file @ bb2f418b
......@@ -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()
......
Ccore/shp.plt deleted 100644 → 0
+ 0
− 53
View file @ bbfd67c2
#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
Lbracket/topo.py
+ 1
− 1
View file @ bb2f418b
......@@ -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()
......
Team25/shape.pro
+ 2
− 0
View file @ bb2f418b
......@@ -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}];
......
Team25/shape.py
+ 1
− 1
View file @ bb2f418b
......@@ -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()
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment