Skip to content
Snippets Groups Projects
Commit 713202bd authored by François Henrotte's avatar François Henrotte
Browse files

reshape optimizer as an independent module

parent ebbec416
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Ccore/OnelabOptimize.geo 0 → 100644
+ 69
− 0
View file @ 713202bd
/*
This file is a generic block of commands
for the automatic computation of velocity fields.
This velocity is used for the computation of sensitivities [1].
It should be included at the bottom
of all shape optimisation .geo files.
[1]:Erin Kuci, François Henrotte, Pierre Duysinx, and Christophe Geuzaine.
``Design sensitivity analysis for shape optimization based on the Lie derivative'',
Computer Methods in Applied Mechanics and Engineering 317 (2017), pp. 702 –722.
issn: 0045-7825.
*/
DefineConstant[
PerturbMesh = 0
VelocityTag = -1
Perturbation = 1e-10
modelpath = CurrentDir
];
If(PerturbMesh == 1)
Printf("Computing velocity field ...");
modelName = StrPrefix(StrRelative(General.FileName));
SyncModel;
// Merge the original mesh
Merge StrCat(modelpath, modelName, ".msh");
// create a view with the original node coordinates as a vector dataset
Plugin(NewView).NumComp = 3;
Plugin(NewView).Run;
Plugin(ModifyComponents).View = 0;
Plugin(ModifyComponents).Expression0 = "x";
Plugin(ModifyComponents).Expression1 = "y";
Plugin(ModifyComponents).Expression2 = "z";
Plugin(ModifyComponents).Run;
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 "*";
// Create a view with the perturbed node coordinates as vector dataset
Plugin(NewView).NumComp = 3;
Plugin(NewView).Run;
Plugin(ModifyComponents).View = 1;
Plugin(ModifyComponents).Expression0 = "x";
Plugin(ModifyComponents).Expression1 = "y";
Plugin(ModifyComponents).Expression2 = "z";
Plugin(ModifyComponents).Run;
Save View[1] StrCat(modelpath, Sprintf("view_%g.msh",1));
// compute the velocity field by subtracting the two vector datasets
Plugin(ModifyComponents).View = 0;
Plugin(ModifyComponents).OtherView = 1;
Plugin(ModifyComponents).Expression0 = Sprintf("(w0 - v0)/(%g)", Perturbation);
Plugin(ModifyComponents).Expression1 = Sprintf("(w1 - v1)/(%g)", Perturbation);
Plugin(ModifyComponents).Expression2 = Sprintf("(w2 - v2)/(%g)", Perturbation);
Plugin(ModifyComponents).Run;
View.Name = "velocity";
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);
EndIf
Ccore/shp.py → Ccore/OnelabOptimize.py
+ 261
− 0
View file @ 713202bd
# Open this file with Gmsh
# - interactively with File->Open->shp.py,
# - in the command line with 'gmsh shape.py'
import numpy as np
import conveks
import onelab
from scipy.optimize import minimize
modelName = 'ccore'
import scipy.optimize as spo
c = onelab.client(__file__)
# get gmsh and getdp locations from Gmsh options
mygmsh = c.getString('General.ExecutableFileName')
mygetdp = ''
for s in range(9):
n = c.getString('Solver.Name' + str(s))
if(n == 'GetDP'):
mygetdp = c.getString('Solver.Executable' + str(s))
break
if(not len(mygetdp)):
c.sendError('This appears to be the first time you are trying to run GetDP')
c.sendError('Please run a GetDP model interactively once with Gmsh to ' +
'initialize the solver location')
exit(0)
c.sendInfo('Will use gmsh={0} and getdp={1}'.format(mygmsh, mygetdp))
# append correct path to model file names
file_geo = c.getPath(modelName+'.geo')
file_msh = c.getPath(modelName+'.msh')
file_pro = c.getPath(modelName+'.pro')
# build command strings with appropriate parameters and options
mygmsh += ' -p 0 -v 2 -parametric ' + file_geo + ' '
mygetdp += ' -p 0 ' + file_pro + ' -msh ' + file_msh + ' '
# load geometry in GUI to have something to look at
c.openProject(file_geo)
c.sendCommand('Solver.AutoMesh = -1;')
# dry getdp run (without -solve or -pos option) to get model parameters in the GUI
c.runSubClient('myGetDP', mygetdp)
# end of check phase. We're done if we do not start the optimization
if c.action == 'check' :
exit(0)
## def readSimpleTable(path):
## with open(path) as FileObj:
## return np.array([float(lines.split()[-1]) for lines in FileObj])
x = {}
def setNodeTable(var, nodes, val):
data = np.ravel(np.column_stack((nodes, val)))
......@@ -110,20 +63,16 @@ def grad(xFromOpti):
grads.append(c.getNumber('Optimization/Results/dwdtau_'+str(dv)))
return np.asarray(grads)
Nfeval = 1
def callbackF(Xi):
global Nfeval
print '{0:4d} {1: 3.6f} {2: 3.6f} {3: 3.6f}'.format(Nfeval, Xi[0], fobj(Xi), grad(Xi)[0])
Nfeval += 1
def test(n):
if n==0: return
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])
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
......@@ -133,101 +82,180 @@ def test(n):
print "%4d" % step,
for item in tab[step] : print "%3.6f" % item,
print
return tab
return
def callbackF(Xi):
print Xi
if (c.getString(c.name+'/Action') == 'stop'):
raise StopIteration
# List of design parameters with format
# [ "Onelab name", value, Min, Max, "variable name" ]
designParameters={
0 : ['Optimization/Parameters/D', 0.02, 0.001, 0.050, 'D'],
1 : ['Optimization/Parameters/E', 0.01 , 0.01, 0.15 , 'E']
}
Optimizers={
0 : 'scipy',
1 : 'mma'
}
x = {}
index = int(c.getNumber('Optimization/Sweep on parameter'))-1
Optimizer = Optimizers[int(c.getNumber('Optimization/Optimizer'))]
# define optimization parameters as Onelab parameter (editable in the GUI)
maxIter = c.defineNumber('Optimization/00Max iterations', value=100, visible=(Optimizer=='mma'))
maxChange = c.defineNumber('Optimization/01Max change', value=1e-1, visible=(Optimizer=='mma'))
if index >= 0: # parameter sweep and optimisation over one selected design parameter
x[0] = designParameters[index];
test(int(c.getNumber('Optimization/Steps in range')))
res = minimize(fobj, [x[0][1]], jac=grad, bounds=[(x[0][2],x[0][3])],\
callback=callbackF, method = 'L-BFGS-B', tol=None,\
options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
print res
else: # optimisation over full parameter space
if Optimizer=='scipy':
for label, var in designParameters.iteritems():
x[label] = var
values = [designParameters[i][1] for i in range(len(designParameters))]
bounds = [(designParameters[i][2], designParameters[i][3]) for i in range(len(designParameters))]
res = minimize(fobj, values, jac=grad, bounds=bounds, callback=callbackF,\
method = 'L-BFGS-B', tol=None, options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
def Optimize(Optimizer):
print Optimizer
print x
if Optimizer==0: # 'scipy'
values = [var[1] for i,var in x.iteritems()]
bounds = [(var[2], var[3]) for i,var in x.iteritems()]
print values, bounds
try:
res = spo.minimize(fobj, values, jac=grad, bounds=bounds,\
callback=callbackF,\
method = 'L-BFGS-B', tol=None,\
options={'ftol': 1e-6, 'gtol': 1e-3, 'disp':True} )
except StopIteration:
print 'Stopped by Gmsh'
return
print res
elif Optimizer=='mma':
numVariables = len(designParameters.keys())
elif Optimizer==1: # 'mma'
numVariables = len(x)
initialPoint = np.zeros(numVariables)
lowerBound = np.zeros(numVariables)
upperBound = np.zeros(numVariables)
for label, var in designParameters.iteritems():
for label, var in x.iteritems():
x[label] = var
initialPoint[label] = var[1]
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', 0) # 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', 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
# Get iteration count (here it will be 1 - could be different in case of restart)
it = conveks.mma.getOuterIteration()
change = 1.0
while it <= maxIter and c.getString(c.name+'/Action') != 'stop' and change > 1e-4 :
change_mma = 1.0
while it <= maxIter :
c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
attributes={'Highlight':'LightYellow'})
xFromMMA = np.array(conveks.mma.getCurrentPoint())
objective = fobj(xFromMMA)
grad_objective = grad(xFromMMA)
constraints = np.array([np.sum(xFromMMA)/100.0-1.0])
grad_constraints = np.ones(numVariables)/100.0
if it == 1: fscale = 1.0 / np.abs(objective)
objective *= fscale
grad_objective *= fscale
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
grad_obj = grad(xFromMMA)
if it == 1:
fscale = 1.0
# fscale = 1.0 / np.abs(obj)
print 'Initial: %4d %3.6f' %(it, obj),
if it > 1:
change_rel = np.abs(obj-obj_prec)*fscale
print '%4d %3.6f %3.6f %3.6f' %(it, obj, change_rel, change_mma),
if change_rel < maxChange: break
print xFromMMA, grad_obj
obj_prec = obj
objective = obj*fscale
grad_objective = grad_obj*fscale
print it, change, xFromMMA, objective, grad_objective
c.sendInfo('Optimization: it. {}, obj. {}, constr. {}'.format(it,objective,constraints[0]))
# call MMA and update the current point
conveks.mma.setMoveLimits(lowerBound, upperBound, 1e-0) # parametre interessant
conveks.mma.updateCurrentPoint(constraints,grad_objective,grad_constraints)
change = conveks.mma.getDesignChange()
#conveks.mma.setMoveLimits(lowerBound, upperBound, 1e-4) # parametre interessant
conveks.mma.updateCurrentPoint(constraints, grad_objective, grad_constraints)
change_mma = conveks.mma.getDesignChange()
if change_mma < maxChange: break
it = conveks.mma.getOuterIteration()
print
print 'Converged: %4d %3.6f %3.6f %3.6f' %(it, obj, change_rel, change_mma),
print xFromMMA, grad_obj
# This should be called at the end
conveks.mma.finalize()
def minimize(name):
global x, modelName, file_geo, file_pro, file_msh, mygmsh, mygetdp, maxIter, maxChange
modelName = name
# get gmsh and getdp locations from Gmsh options
mygmsh = c.getString('General.ExecutableFileName')
mygetdp = ''
for s in range(9):
n = c.getString('Solver.Name' + str(s))
if(n == 'GetDP'):
mygetdp = c.getString('Solver.Executable' + str(s))
break
if(not len(mygetdp)):
c.sendError('This appears to be the first time you are trying to run GetDP')
c.sendError('Please run a GetDP model interactively once with Gmsh to ' +
'initialize the solver location')
exit(0)
c.sendInfo('Will use gmsh={0} and getdp={1}'.format(mygmsh, mygetdp))
# append correct path to model file names
file_geo = c.getPath(modelName+'.geo')
file_msh = c.getPath(modelName+'.msh')
file_pro = c.getPath(modelName+'.pro')
# build command strings with appropriate parameters and options
mygmsh += ' -p 0 -v 2 -parametric ' + file_geo + ' '
mygetdp += ' -p 0 ' + file_pro + ' -msh ' + file_msh + ' '
# load geometry in GUI to have something to look at
c.openProject(file_geo)
c.sendCommand('Solver.AutoMesh = -1;')
# dry getdp run (without -solve or -pos option) to get model parameters in the GUI
c.runSubClient('myGetDP', mygetdp)
# get design space from Onelab
sweepOnParameter = int(c.getNumber('Optimization/Sweep on parameter'))
if sweepOnParameter:
designSpace = [ sweepOnParameter ]
else:
designSpace = c.getNumberChoices('Optimization/Design space')
designParams = c.getStringChoices('Optimization/Design parameters')
designParameters = {}
iParam = 0
for i in designSpace:
label = designParams[int(i)-1]
myrange = c.getNumberRange(label)
param = [label]
param.extend(myrange)
param.extend([label.split('/')[-1]])
designParameters[iParam]=param
iParam +=1
print
print 'Working with design space of dimension', len(designParameters), ' :'
for i, val in designParameters.iteritems():
print i, val
## Optimizer2 = c.defineString('Optimization/Algorithm', value='scipy',
## choices=('mma', 'scipy'))
## print Optimizer2
maxIter = c.defineNumber('Optimization/00Max iterations', value=30, visible=1)
maxChange = c.defineNumber('Optimization/01Max change', value=1e-4, visible=1)
# end of check phase. We're done if we do not start the optimization
if c.action == 'check' :
exit(0)
Optimizer = c.getNumber('Optimization/Algorithm')
x=designParameters # x is global
if sweepOnParameter:
test(int(c.getNumber('Optimization/Steps in range')))
Optimize(Optimizer)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment