From 8c78d32ebe5e134c83d56ab7620f0903e61fd8fc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Henrotte?= <francois.henrotte@uclouvain.be>
Date: Mon, 15 Oct 2018 14:55:10 +0200
Subject: [PATCH] added comments and explanations
---
Ccore/ccore.geo | 7 +++++
Ccore/ccore.pro | 68 ++++++++++++++++++++++++++++++++++++++++--
Ccore/ccore_common.pro | 40 +++++++++++++------------
Ccore/shp.py | 14 ++++-----
4 files changed, 100 insertions(+), 29 deletions(-)
diff --git a/Ccore/ccore.geo b/Ccore/ccore.geo
index 348dff3..d5440d1 100644
--- a/Ccore/ccore.geo
+++ b/Ccore/ccore.geo
@@ -106,8 +106,15 @@ Physical Line("DIR", 11) = {1 ... 4};
/*
The remainder of the 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 copied unmodified 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[
diff --git a/Ccore/ccore.pro b/Ccore/ccore.pro
index e32cb7f..d7054f6 100644
--- a/Ccore/ccore.pro
+++ b/Ccore/ccore.pro
@@ -1,3 +1,66 @@
+/* -------------------------------------------------------------------
+ Tutorial: Optimization of ferromagnetic C-core
+
+ Features:
+ - Nonlinear magnetostatic model of a C-core
+ - Imposed current density or imposed current
+ - Air gap and core width as design parameters
+ - Prescribed By in the air gap as objective function
+ - Direct approach of the shape sensitivity analysis
+ - Sensivity analysis with Lie derivative
+
+ To compute the solution in a terminal:
+ getdp shape.pro -solve GetPerformancesAndGradients
+
+ To compute the solution interactively from the Gmsh GUI:
+ File > Open > shape.pro
+ Run (button at the bottom of the left panel)
+ ------------------------------------------------------------------- */
+
+/*
+ This tutorial deals with shape sensitivity analysis in Magnetostatics.
+ The model consists of a ferromagnetic C-core powered
+ by a coil with 'NbTurns' windings suplied with a DC current.
+ The ferromagnetic behavior of the C-core material can be
+ represented by either a linear of a nonlinear saturable law.
+ The power supply is set to maintain either a constant total current
+ or a constant current density as the design parameters vary.
+ The optimisation problem aims at having a prescribed value
+ for the By component of the magnetic flux density
+ in the middle of the air gap.
+ The two design variables selected for achieving this objective are
+ - the air gap width 'D',
+ - the core width 'E'.
+
+ The 'shp.py' is organised so as to allow either
+ - to check the correctness of the sensitivity computation,
+ - effectively solve the optimisation problem.
+
+ The computation of sensitivities is a delicate step of the optimisation,
+ which can however be one hundred percent validated by finite difference.
+ This tutorial is thus dedicated to both the explanation and verification
+ of the direct computation of sensitivities in Magnetostatics.
+
+ Communication between the optimisation loop ('shp.py')
+ and the Onelab model is done via the Onelab database.
+ Design variables are defined as Onelab variables whose name
+ has the reserved prefix 'Optimization/Parameters'.
+ Similarly, velocities and sensibilities are also exchanged as
+ Onelab variables with reserved prefixes 'Optimization/Results/*'.
+
+ Finally, an extra mechanism is implemented
+ to apply a temporary small perturbation of a design parameter
+ in order to compute the associated velocity field
+ without interfering with the value stored in the Onelab database.
+ This affectation is done with a '-setnumber' line argument
+ when calling Gmsh.
+ In consequence, the definition of the design variables
+ as Onelab parameters must be enclosed in a 'DefineConstant[]' statement.
+ Remember that affectations within a 'DefineConstant[]'
+ are ignored if the variable already exists because,
+ and in particular if it has been set by a '-setnumber' line argument.
+ */
+
Include "ccore_common.pro";
DefineConstant[
@@ -8,6 +71,7 @@ DefineConstant[
OptiIterNumber = 0
];
+// Onelab parameters of the C-core model
NL_tol_abs = 1e-8; // absolute tolerance on residual for noninear iterations
NL_tol_relax = 1.0; // relaxation on residual for noninear iterations
@@ -37,9 +101,7 @@ Else
CurrentDensity = Mmf / CoilSection;
EndIf
-By_target = 1.;
-Px_target = CoreX + A - E/2 ;
-Py_target = CoreY + B/2 ;
+
Group {
// Physical regions (in capital letters):
diff --git a/Ccore/ccore_common.pro b/Ccore/ccore_common.pro
index 3898872..6ef1cab 100644
--- a/Ccore/ccore_common.pro
+++ b/Ccore/ccore_common.pro
@@ -1,28 +1,14 @@
-/*
- This tutorial deals with shape sensitivity analysis in Magnetostatics.
- The model consists of a ferromagnetic C-core powered
- by a coil with 'NbTurns' windings suplied with a DC current.
- The ferromagnetic behavior of the C-core material can be
- represented by either a linear of a nonlinear saturable law.
- The power supply is set to maintain either a constant total current
- or a constant current density as the design parameters vary.
-
- Three geometrical parametrs are
-
- */
-
+// Onelab parameters related with the optimisation
ParamIndex =
DefineNumber[0, Name "Optimization/Sweep on parameter",
- Choices {0="None",
+ Choices {0="Optimize on D and E",
1="D: air gap length",
- 2="E: C-core thickness",
- 3="F: Core-coil gap"} ];
+ 2="E: C-core thickness"} ];
StepsInRange =
DefineNumber[20, Name "Optimization/Steps in range", Visible ParamIndex];
-
w =
DefineNumber[0, Name "Optimization/Results/w",
Graph "030000000000000000000000000000000000"];
@@ -32,15 +18,28 @@ dwdt =
+// Design parameters 'D' and 'E' are defined
+// within a 'DefineConstant[]' statement to be contrallable
+// from the calling optimisation loop with a '-setnumber' line argument
+// for the computation of the velocity field.
+
DefineConstant[
D = DefineNumber[0.01, Name "Optimization/Parameters/D",
Min 0.001, Max 0.050]
E = DefineNumber[0.1 , Name "Optimization/Parameters/E",
Min 0.01 , Max 0.15]
- F = DefineNumber[0.01, Name "Optimization/Parameters/F",
- Min 0.001, Max 0.090]
];
+// A number of geometrical parameters of the physical model
+// are also used in the definition of the finite element model
+// and the optimisation problem.
+// Definitions for 'D' and 'E' must be in a 'DefineConstant[]' statement
+// for the same reason as above
+// (they also must be overriden by the '-setnumber' line argument)
+// They are however all gathered in the '_common.pro' file
+// and in the 'DefineConstant[]' statement for the sake of clarity
+// and a tidy organisation of the model files.
+
DefineConstant[
// all dimensions in meter
@@ -60,3 +59,6 @@ DefineConstant[
R=1 // refinement factor: R=5 far raw mesh, R=1 normal mesh
];
+By_target = 1.;
+Px_target = CoreX + A - E/2 ;
+Py_target = CoreY + B/2 ;
diff --git a/Ccore/shp.py b/Ccore/shp.py
index 32aa8d2..bb23d55 100644
--- a/Ccore/shp.py
+++ b/Ccore/shp.py
@@ -1,6 +1,7 @@
-# Open this file with Gmsh (interactively with File->Open->shape.py, or on the
-# command line with 'gmsh shape.py')
-#
+# 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
@@ -139,19 +140,18 @@ def test(n):
return tab
-# x is the list of design parameters with format
+# List of design parameters with format
# [ "Onelab name", value, Min, Max, "variable name" ]
designParameters={
0 : ['Optimization/Parameters/D', 0.01, 0.001, 0.050, 'D'],
1 : ['Optimization/Parameters/E', 0.01 , 0.01, 0.15 , 'E']
-# 2 : ['Optimization/Parameters/F', 0.01, 0.001, 0.090, 'F']
}
x = {}
index = int(c.getNumber('Optimization/Sweep on parameter'))-1
-if index >= 0:
+if index >= 0: # parameter sweep and optimisation over one design parameter
x[0] = designParameters[index];
test(int(c.getNumber('Optimization/Steps in range')))
#callback=callbackF,
@@ -160,7 +160,7 @@ if index >= 0:
options={'ftol': 1e-5, 'gtol': 1e-3, 'disp':True} )
print res
-else:
+else: # optimisation over full parameter space
for label, var in designParameters.iteritems():
x[label] = var
values = [designParameters[i][1] for i in range(len(designParameters))]
--
GitLab