From 713202bd8fd682e90f85d6951622a944376f0ac5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois=20Henrotte?= <francois.henrotte@uclouvain.be> Date: Thu, 25 Oct 2018 11:42:29 +0200 Subject: [PATCH] reshape optimizer as an independent module --- Ccore/OnelabOptimize.geo | 69 +++++++++++ Ccore/OnelabOptimize.py | 261 +++++++++++++++++++++++++++++++++++++++ Ccore/shp.py | 233 ---------------------------------- 3 files changed, 330 insertions(+), 233 deletions(-) create mode 100644 Ccore/OnelabOptimize.geo create mode 100644 Ccore/OnelabOptimize.py delete mode 100644 Ccore/shp.py diff --git a/Ccore/OnelabOptimize.geo b/Ccore/OnelabOptimize.geo new file mode 100644 index 0000000..d4ce25f --- /dev/null +++ b/Ccore/OnelabOptimize.geo @@ -0,0 +1,69 @@ +/* + 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 diff --git a/Ccore/OnelabOptimize.py b/Ccore/OnelabOptimize.py new file mode 100644 index 0000000..63576ba --- /dev/null +++ b/Ccore/OnelabOptimize.py @@ -0,0 +1,261 @@ +import numpy as np +import conveks +import onelab +import scipy.optimize as spo + +c = onelab.client(__file__) +x = {} + +def setNodeTable(var, nodes, val): + data = np.ravel(np.column_stack((nodes, val))) + data = np.insert(data, 0, float(len(nodes))) + c.setNumber(var, values=data, visible=0) + +def getVelocityField(xvar, perturb=1e-7): + for id, var in xvar.iteritems(): + c.runNonBlockingSubClient('myGmsh', mygmsh \ + + ' -setnumber PerturbMesh 1'\ + + ' -setnumber Perturbation '+str(perturb)\ + + ' -setnumber VelocityTag '+str(id)\ + + ' -setnumber '+var[-1]+' '+str(var[1]+perturb)\ + + ' -run') + c.waitOnSubClients() + # send the x,y,z components of each velocity field + for label, var in x.iteritems(): + d = c.getNumbers('Optimization/Results/velocity_'+str(label)) + if label==0: + nodes = np.array(d[1::4]) + setNodeTable('Optimization/Results/velocity_'+str(label)+'_1',nodes,np.array(d[2::4])) + setNodeTable('Optimization/Results/velocity_'+str(label)+'_2',nodes,np.array(d[3::4])) + setNodeTable('Optimization/Results/velocity_'+str(label)+'_3',nodes,np.array(d[4::4])) + c.setNumber(var[0],value=x[label][1]) + +def fobj(xFromOpti): + # send the current point to GetDP model + for label, var in x.iteritems(): + x[label][1] = xFromOpti[label] + c.setNumber(var[0],value=x[label][1]) + + # reload the geometry in GUI to see the geometrical changes + c.openProject(file_geo) + + c.runSubClient('myGmsh', mygmsh + ' -2') + getVelocityField(x) + c.runSubClient('myGetDP', mygetdp + '-solve GetPerformances' \ + + ' - setnumber OuterIteration ' + str(0)) # it !! + return c.getNumber('Optimization/Results/w') + + +def grad(xFromOpti): + # send the current point to GetDP model + for label, var in x.iteritems(): + x[label][1] = xFromOpti[label] + c.setNumber(var[0],value=x[label][1]) + + # as well as their sensitivity with respect to design variables at `x' + for dv, var in x.iteritems(): + c.runNonBlockingSubClient('myGetDP', mygetdp \ + + ' -setnumber VelocityTag '+str(dv)\ + + ' -solve GetGradient_wrt_dv') + c.waitOnSubClients() + grads=[] + for dv, var in x.iteritems(): + grads.append(c.getNumber('Optimization/Results/dwdtau_'+str(dv))) + return np.asarray(grads) + +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]) + 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]) + print "%4d" % step, + for item in tab[step] : print "%3.6f" % item, + print + return + +def callbackF(Xi): + print Xi + if (c.getString(c.name+'/Action') == 'stop'): + raise StopIteration + + +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==1: # 'mma' + numVariables = len(x) + initialPoint = np.zeros(numVariables) + lowerBound = np.zeros(numVariables) + upperBound = np.zeros(numVariables) + 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', 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_mma = 1.0 + while it <= maxIter : + c.setNumber('Optimization/01Current iteration', value=it, readOnly=1, + attributes={'Highlight':'LightYellow'}) + + 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 + + c.sendInfo('Optimization: it. {}, obj. {}, constr. {}'.format(it,objective,constraints[0])) + + # call MMA and update the current point + #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) + + + + + diff --git a/Ccore/shp.py b/Ccore/shp.py deleted file mode 100644 index e9eaaca..0000000 --- a/Ccore/shp.py +++ /dev/null @@ -1,233 +0,0 @@ -# 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' - -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]) - -def setNodeTable(var, nodes, val): - data = np.ravel(np.column_stack((nodes, val))) - data = np.insert(data, 0, float(len(nodes))) - c.setNumber(var, values=data, visible=0) - -def getVelocityField(xvar, perturb=1e-7): - for id, var in xvar.iteritems(): - c.runNonBlockingSubClient('myGmsh', mygmsh \ - + ' -setnumber PerturbMesh 1'\ - + ' -setnumber Perturbation '+str(perturb)\ - + ' -setnumber VelocityTag '+str(id)\ - + ' -setnumber '+var[-1]+' '+str(var[1]+perturb)\ - + ' -run') - c.waitOnSubClients() - # send the x,y,z components of each velocity field - for label, var in x.iteritems(): - d = c.getNumbers('Optimization/Results/velocity_'+str(label)) - if label==0: - nodes = np.array(d[1::4]) - setNodeTable('Optimization/Results/velocity_'+str(label)+'_1',nodes,np.array(d[2::4])) - setNodeTable('Optimization/Results/velocity_'+str(label)+'_2',nodes,np.array(d[3::4])) - setNodeTable('Optimization/Results/velocity_'+str(label)+'_3',nodes,np.array(d[4::4])) - c.setNumber(var[0],value=x[label][1]) - -def fobj(xFromOpti): - # send the current point to GetDP model - for label, var in x.iteritems(): - x[label][1] = xFromOpti[label] - c.setNumber(var[0],value=x[label][1]) - - # reload the geometry in GUI to see the geometrical changes - c.openProject(file_geo) - - c.runSubClient('myGmsh', mygmsh + ' -2') - getVelocityField(x) - c.runSubClient('myGetDP', mygetdp + '-solve GetPerformances' \ - + ' - setnumber OuterIteration ' + str(0)) # it !! - return c.getNumber('Optimization/Results/w') - - -def grad(xFromOpti): - # send the current point to GetDP model - for label, var in x.iteritems(): - x[label][1] = xFromOpti[label] - c.setNumber(var[0],value=x[label][1]) - - # as well as their sensitivity with respect to design variables at `x' - for dv, var in x.iteritems(): - c.runNonBlockingSubClient('myGetDP', mygetdp \ - + ' -setnumber VelocityTag '+str(dv)\ - + ' -solve GetGradient_wrt_dv') - c.waitOnSubClients() - grads=[] - for dv, var in x.iteritems(): - 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): - 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]) - print "%4d" % step, - for item in tab[step] : print "%3.6f" % item, - print - - return tab - - -# 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} ) - print res - - elif Optimizer=='mma': - numVariables = len(designParameters.keys()) - initialPoint = np.zeros(numVariables) - lowerBound = np.zeros(numVariables) - upperBound = np.zeros(numVariables) - for label, var in designParameters.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 - - # 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 : - 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 - - 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() - it = conveks.mma.getOuterIteration() - - # This should be called at the end - conveks.mma.finalize() - - -- GitLab