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

tuto shape optimization in Magnetostatics

parent 241e73fc
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Ccore/ccore.geo 0 → 100644
+ 167
− 0
View file @ f0cb6778
/*
All geometrical parameters of the model are used in both
Gmsh (CAO) and GetDP (shape optimisation).
They are therefore defined in the file "ccore_common.pro".
*/
Include "ccore_common.pro";
// The mesh refinement factor is rather a global user parameter,
// therefore defined as a Onelab variable.
R = DefineNumber[1, Name "Model/Parameters/Mesh refinement factor",
Help "R=5 far raw mesh, R=1 normal mesh"];
// Center C-core in air box
CoreX=(L-A)/2; // X-position of C-core's bottom left corner
CoreY=(L-B)/2; // Y-position of C-core's bottom left corner
// external boundary
lc1=L/10.*R;
Point(1) = {0,0,0,lc1};
Point(2) = {L,0,0,lc1};
Point(3) = {L,L,0,lc1};
Point(4) = {0,L,0,lc1};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Curve Loop(1) = {1 ... 4};
// magnetic C-core
lc2=A/10.*R;
lc3=D/2.*R;
Point(5) = {CoreX,CoreY,0,lc2};
Point(6) = {CoreX+A,CoreY,0,lc2};
Point(7) = {CoreX+A,CoreY+(B-D)/2.,0,lc3};
Point(8) = {CoreX+A-E,CoreY+(B-D)/2.,0,lc3};
Point(9) = {CoreX+A-E,CoreY+E,0,lc2};
Point(10) = {CoreX+E,CoreY+E,0,lc2};
Point(11) = {CoreX+E,CoreY+B-E,0,lc2};
Point(12) = {CoreX+A-E,CoreY+B-E,0,lc2};
Point(13) = {CoreX+A-E,CoreY+(B+D)/2.,0,lc3};
Point(14) = {CoreX+A,CoreY+(B+D)/2.,0,lc3};
Point(15) = {CoreX+A,CoreY+B,0,lc2};
Point(16) = {CoreX,CoreY+B,0,lc2};
Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {7,8};
Line(8) = {8,9};
Line(9) = {9,10};
Line(10) = {10,11};
Line(11) = {11,12};
Line(12) = {12,13};
Line(13) = {13,14};
Line(14) = {14,15};
Line(15) = {15,16};
Line(16) = {16,5};
Curve Loop(2) = {5 ... 16};
// inductors
lc4=lc2; //F/2.*R;
Point(17) = {CoreX+E+F,CoreY+E+F,0,lc4};
Point(18) = {CoreX+E+F+G,CoreY+E+F,0,lc4};
Point(19) = {CoreX+E+F+G,CoreY+B-E-F,0,lc4};
Point(20) = {CoreX+E+F,CoreY+B-E-F,0,lc4};
Line(17) = {17,18};
Line(18) = {18,19};
Line(19) = {19,20};
Line(20) = {20,17};
Curve Loop(3) = {17 ... 20};
Point(21) = {CoreX-F-G,CoreY+E+F,0,lc4};
Point(22) = {CoreX-F,CoreY+E+F,0,lc4};
Point(23) = {CoreX-F,CoreY+B-E-F,0,lc4};
Point(24) = {CoreX-F-G,CoreY+B-E-F,0,lc4};
Line(21) = {21,22};
Line(22) = {22,23};
Line(23) = {23,24};
Line(24) = {24,21};
Curve Loop(4) = {21 ... 24};
Plane Surface(1) = {1,2,3,4};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Physical Surface("AIR", 1) = 1;
Physical Surface("CORE", 2) = 2;
Physical Surface("COILP", 3) = 3;
Physical Surface("COILN", 4) = 4;
Physical Line("DIR", 11) = {1 ... 4};
/*
The remainder of the file is a generic block of commands
for the automatic computation of velocity fields.
It should be copied unmodified at the bottom
of all shape optimisation .geo files.
*/
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/ccore.pro 0 → 100644
+ 457
− 0
View file @ f0cb6778
Include "ccore_common.pro";
DefineConstant[
R_ = {"GetPerformances", Name "GetDP/1ResolutionChoices", Visible 0}
C_ = {"-solve -v 5 -v2", Name "GetDP/9ComputeCommand", Visible 0}
P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
VelocityTag = -1
OptiIterNumber = 0
];
NL_tol_abs = 1e-8; // absolute tolerance on residual for noninear iterations
NL_tol_relax = 1.0; // relaxation on residual for noninear iterations
NL_iter_max = 50; // maximum number of noninear iterations
Flag_NL =
DefineNumber[ 1, Name "Model/Parameters/Non linear core material", Choices {0,1}];
murCore =
DefineNumber[300, Name "Model/Parameters/Mur core", Visible 1,//!Flag_NL,
Help "Magnetic relative permeability of Core"];
Flag_Jfixed =
DefineNumber[ 1, Name "Model/Parameters/Fixed current density", Choices {0,1}];
CurrentDensity = 1e6 *
DefineNumber[1, Name "Model/Parameters/Current Density", Visible Flag_Jfixed,
Help "Current density in coil [A/mm2]"];
NbTurns = 300;
Current =
DefineNumber[CurrentDensity*CoilSection_ref/NbTurns,
Name "Model/Parameters/Current [A]", Visible !Flag_Jfixed,
Help "Current in coil [A]"];
CoilSection = G*(B-2*E-2*F);
If(Flag_Jfixed)
Mmf = CurrentDensity * CoilSection;
Current = Mmf / NbTurns;
Else
Mmf = NbTurns * Current;
CurrentDensity = Mmf / CoilSection;
EndIf
By_target = 1.;
Px_target = CoreX + A - E/2 ;
Py_target = CoreY + B/2 ;
Group {
// Physical regions (in capital letters):
AIR = Region[ 1 ];
CORE = Region[ 2 ];
COILP = Region[ 3 ];
COILN = Region[ 4 ];
NOFLUX = Region[ 11 ];
// Abstract regions
Vol_Mag = Region[ {AIR, CORE, COILP, COILN} ];
Vol_S_Mag = Region[ {COILP, COILN} ];
Sur_Dir_Mag = Region[ {NOFLUX} ];
Sur_Neu_Mag = Region[ {} ];
Vol_NL_Mag = Region[ {} ];
If(Flag_NL)
Vol_NL_Mag = Region[ {CORE} ];
EndIf
Vol_L_Mag = Region[ {Vol_Mag,-Vol_NL_Mag} ];
}
Function {
mu0 = 4.e-7 * Pi ;
nu0 = 1. / mu0 ;
nu [ Region[{AIR, COILP, COILN}] ] = nu0;
Mat_h_r = { 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 125, 150, 175, 200, 250,
300, 400, 500, 600, 700, 800, 900, 1000, 1250, 1500, 2000, 2500, 5000,
7500, 10000, 15000, 20000, 59000, 174000, 514000, 1520000, 4470000,
13200000, 38900000, 115000000, 339000000, 1000000000 } ;
Mat_b_r = { 0.0, 0.194880829963, 0.377143018857, 0.537767739762, 0.672888260835,
0.783043000477, 0.871342430831,0.941778611986, 0.998183303557, 1.04378111223,
1.08110469369, 1.14963767549, 1.19607212343, 1.22964695907, 1.25515221835,
1.29162498935, 1.31678879432, 1.35015120537, 1.37220092877, 1.38859114656,
1.4017440574, 1.41287024565, 1.42264180514, 1.43146158921, 1.45082466146,
1.46784549989, 1.49819370601, 1.52578650709, 1.64314027719, 1.73458485332,
1.8039068939,1.89568786291, 1.95213815187, 2.1390774927, 2.45827909293,
3.32303272825, 5.85485500678, 13.2701832298, 35.2114648741, 99.8027446541,
291.062951228, 854.036370229, 2515.3105707 } ;
Mat_b2_r = Mat_b_r()^2;
Mat_nu_r = Mat_h_r()/Mat_b_r();
Mat_nu_r(0) = Mat_nu_r(1);
Mat_nu_b2_r = ListAlt[Mat_b2_r(), Mat_nu_r()];
nu_interp_r[] = InterpolationLinear[ SquNorm[$1] ]{Mat_nu_b2_r()};
dnudb2_interp_r[] = dInterpolationLinear[SquNorm[$1]]{Mat_nu_b2_r()};
h_interp_r[] = nu_interp_r[$1] * $1 ;
If(!Flag_NL)
nu [ CORE ] = nu0/murCore;
dhdb_NL[ CORE ] = 0;
mu_analytic = mu0*murCore;
Else
nu [ Vol_NL_Mag ] = nu_interp_r[$1] ;
dhdb_NL[ Vol_NL_Mag ] = 2*dnudb2_interp_r[$1#2]*SquDyadicProduct[#2];
mu_analytic = 1/Mat_nu_r(0);
EndIf
// shape function of J : js[] = nI * Jshape[]
Jshape[ COILP ] = Vector[0,0,+1]/CoilSection;
Jshape[ COILN ] = Vector[0,0,-1]/CoilSection;
js[] = Mmf * Jshape[];
// Analytical solution
RelGap = D/(mu0*E);
RelCore = (2*(A+B-2*E)-D)/(mu_analytic*E);
Flux = NbTurns*Mmf/(RelGap+RelCore); // analytic approximation
SetNumber("Model/Results/Flux analytic", Flux);
//Flux = 0;
If( ParamIndex == 0 ) // D
dFluxdTau = Flux/E*Flux/(NbTurns*Mmf)*(1/mu_analytic-1/mu0);
ElseIf( ParamIndex == 1 ) // E
dFluxdTau = Flux/E*(1.+4*Flux/(NbTurns*Mmf*mu_analytic));
If( Flag_Jfixed )
dFluxdTau -= 2*Flux*G/CoilSection;
EndIf
ElseIf( ParamIndex == 2 ) // F
dFluxdTau = 0;
EndIf
VV[] = Vector[$1,$2,0*$3] ;
dV[] = TensorV[$1,$2,0*$3] ;
Lie2form[] = TTrace[$2]*$1 - Transpose[$2]*$1; // $1=2-form (vector field), $2=dV (tensor)
If(Flag_Jfixed)
LieOf_js[] = Lie2form[js[], $1]; // $1=dV
Else
LieOf_js[] = 0;
EndIf
// Lie derivative of H(B) where B is a 2-form and H a 1-form ($1:{d a}, $2:dV)
LieOf_HB[] = nu[$1] * (Transpose[$2] * $1 - TTrace[$2] * $1 + $2 * $1) ;
LieOf_HB_NL[] = dhdb_NL[$1] * (Transpose[$2] - TTrace[$2] * TensorDiag[1,1,1]) * $1;
}
Group {
Dom_Hcurl_a_Mag_2D = Region[ {Vol_Mag, Sur_Neu_Mag} ];
}
Function{
l_a = ListFromServer["Optimization/Results/a"];
aFromServer[] = ValueFromIndex[]{l_a()};
For i In {1:3}
l_v~{i} = ListFromServer[Sprintf["Optimization/Results/velocity_%g_%g",VelocityTag,i]];
velocity~{i}[] = ValueFromIndex[]{l_v~{i}()};
EndFor
}
Constraint {
{ Name Dirichlet_a_Mag;
Case {
{ Region NOFLUX; Type Assign; Value 0; }
}
}
{ Name aFromServer;
Case {
{ Region Vol_Mag; Type Assign; Value aFromServer[]; }
{ Region NOFLUX; Type Assign; Value 0; }
}
}
For i In {1:3}
{ Name velocity~{i} ;
Case {
{ Region Vol_Mag ; Value velocity~{i}[]; }
}
}
EndFor
}
FunctionSpace {
{ Name Hcurl_a_2D; Type Form1P;
BasisFunction {
{ Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
Support Dom_Hcurl_a_Mag_2D ; Entity NodesOf[ All ]; }
}
Constraint {
{ NameOfCoef ae; EntityType NodesOf;
NameOfConstraint Dirichlet_a_Mag; }
}
}
{ Name Hcurl_a_2D_fullyfixed; Type Form1P;
BasisFunction{
{ Name se1; NameOfCoef ae1; Function BF_PerpendicularEdge;
Support Vol_Mag; Entity NodesOf[All]; }
}
Constraint{
{ NameOfCoef ae1; EntityType NodesOf; NameOfConstraint aFromServer; }
}
}
For i In {1:3}
{ Name H_v~{i} ; Type Form0;
BasisFunction {
{ Name sn ; NameOfCoef un ; Function BF_Node ;
Support Vol_Mag; Entity NodesOf[ All ] ; }
}
Constraint {
{ NameOfCoef un ; EntityType NodesOf ; NameOfConstraint velocity~{i}; }
}
}
EndFor
}
Jacobian {
{ Name Vol ;
Case {
{ Region All ; Jacobian Vol ; }
}
}
{ Name Sur;
Case {
{ Region All; Jacobian Sur; }
}
}
}
Integration {
{ Name Int ;
Case { { Type Gauss ;
Case {
{ GeoElement Point; NumberOfPoints 1; }
{ GeoElement Line; NumberOfPoints 5; }
{ GeoElement Triangle; NumberOfPoints 4; }
{ GeoElement Quadrangle; NumberOfPoints 4; }
}
}
}
}
}
// -------------------------------------------------------------------------
// This resolution solves the direct problem to compute the sensitivity
// of the induction flux with respect to a given design variable.
Formulation {
{ Name MagSta_a; Type FemEquation;
Quantity {
{ Name a; Type Local; NameOfSpace Hcurl_a_2D; }
}
Equation {
Integral { [ nu[{d a}] * Dof{d a} , {d a} ];
In Vol_Mag ; Jacobian Vol; Integration Int; }
Integral { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Integral { [ -js[] , {a} ];
In Vol_S_Mag; Jacobian Vol; Integration Int; }
}
}
}
Resolution {
{ Name GetPerformances;
System {
{ Name SYS; NameOfFormulation MagSta_a;}
}
Operation {
If( OptiIterNumber == 1 )
CreateDir[ResDir];
DeleteFile[StrCat[ResDir,"w.txt"]];
DeleteFile[StrCat[ResDir,"Lie_w.txt"]];
EndIf
InitSolution[SYS];
IterativeLoop[NL_iter_max, NL_tol_abs, NL_tol_relax]{
GenerateJac[SYS];
SolveJac[SYS];
}
PostOperation[Get_ObjectiveConstraints];
}
}
}
PostProcessing {
{ Name MagSta_a_2D; NameOfFormulation MagSta_a;
Quantity {
{ Name az;
Value {
Term { [ CompZ[{a}] ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
{ Name b;
Value {
Term { [ {d a} ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
{ Name by;
Value {
Term { [ CompY[{d a}] ]; In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
{ Name js;
Value {
Term { [ js[] ]; In Vol_S_Mag; Jacobian Vol; }
}
}
{ Name Flux ; Value { Integral {
[ {a} * NbTurns * Jshape[] ] ;
In Vol_S_Mag ; Jacobian Vol ; Integration Int ; } }
}
{ Name fobj_b_at_point;
Value {
Term { [ SquNorm[ CompY[{d a}] - By_target ] ];
In Dom_Hcurl_a_Mag_2D; Jacobian Vol; }
}
}
}
}
}
PostOperation {
{ Name Get_ObjectiveConstraints; NameOfPostProcessing MagSta_a_2D;
Operation{
CreateDir["res"];
//Print[js, OnElementsOf Vol_S_Mag, File "res/js.pos"];
Print[b, OnElementsOf Vol_Mag, File "res/b.pos"];
Print[az, OnElementsOf Vol_Mag, File "res/az.pos"];
Echo[ Str["k= PostProcessing.NbViews-1;",
"View[k].IntervalsType = 3;",
"View[0].NbIso = 30;"],
File "tmp.geo", LastTimeStepOnly] ;
// In memory transfer of the a field for assembling sensitivity equations
Print[az, OnElementsOf Vol_Mag, Format NodeTable, File "",
SendToServer "Optimization/Results/a", Hidden 1];
Print[ Flux[ Vol_S_Mag ], OnGlobal,
Format Table, File > "",
SendToServer "Model/Results/Flux computed"];
Print[ by, OnPoint { Px_target, Py_target, 0},
Format Table, File "",
SendToServer "Model/Results/b probe"];
Print[ fobj_b_at_point, OnPoint { Px_target, Py_target, 0},
Format Table, File > "res/w.txt",
SendToServer "Optimization/Results/w"];
}
}
}
// -------------------------------------------------------------------------
// This resolution solves the direct problem to compute the sensitivity
// of the induction flux with respect to a given design variable.
Formulation {
{ Name DirectFormulationOf_MagSta_a_2D; Type FemEquation ;
Quantity {
{ Name Lva; Type Local; NameOfSpace Hcurl_a_2D; }
{ Name a; Type Local; NameOfSpace Hcurl_a_2D_fullyfixed; }
For i In {1:3}
{ Name v~{i} ; Type Local ; NameOfSpace H_v~{i}; }
EndFor
}
Equation {
Galerkin { [ nu[{d a}] * Dof{d Lva}, {d Lva} ];
In Vol_Mag; Jacobian Vol; Integration Int; }
Galerkin { [ dhdb_NL[{d a}] * Dof{d Lva}, {d Lva} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Galerkin { [LieOf_HB[{d a},dV[{d v_1},{d v_2},{d v_3}]], {d Lva}];
In Vol_Mag; Jacobian Vol; Integration Int; }
Galerkin { [LieOf_HB_NL[{d a},dV[{d v_1},{d v_2},{d v_3}]], {d Lva}];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Integral { [ -LieOf_js[dV[{d v_1},{d v_2},{d v_3}]], {Lva} ] ;
In Vol_S_Mag ; Jacobian Vol; Integration Int; }
Integral { [ 0*Dof{a}, {a} ] ;
In Vol_Mag; Jacobian Vol ; Integration Int ; }
For i In {1:3}
Integral { [ 0*Dof{v~{i}}, {v~{i}} ] ;
In Vol_Mag; Jacobian Vol ; Integration Int ; }
EndFor
}
}
}
Resolution {
{ Name GetGradient_wrt_dv;
System {
{ Name DIR; NameOfFormulation DirectFormulationOf_MagSta_a_2D; }
}
Operation {
InitSolution[DIR];
Generate[DIR];
Solve[DIR];
PostOperation[Get_GradOf_w];
}
}
}
PostProcessing{
{ Name Direct_MagSta; NameOfFormulation DirectFormulationOf_MagSta_a_2D;
PostQuantity {
{ Name az;
Value{
Term{ [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; }
}
}
{ Name VV;
Value{
Term{ [ Vector[{v_1},{v_2},{v_3}] ]; In Vol_Mag ; Jacobian Vol;}
}
}
{ Name Lie_az;
Value {
Term { [ CompZ[ {Lva}] ] ; In Vol_Mag ; Jacobian Vol; }
}
}
{ Name Lie_b_at_point ; Value {
Term { [ 2 * (CompY[{d a}] - By_target) * CompY[{d Lva}] ] ;
In Vol_Mag; Jacobian Vol; }
}
}
}
}
}
PostOperation{
{ Name Get_GradOf_w; NameOfPostProcessing Direct_MagSta;
Operation {
CreateDir["res"];
Print[ VV, OnElementsOf Vol_Mag,File "res/velocity.pos" ];
Echo[ Str["View[PostProcessing.NbViews-1].GlyphLocation = 2;"],
File "tmp.geo", LastTimeStepOnly] ;
//Print[az, OnElementsOf Vol_Mag, File "res/azstar.pos"];
Echo[ Str["k= PostProcessing.NbViews-1;",
"View[k].IntervalsType = 3;",
"View[0].NbIso = 30;"],
File "tmp.geo", LastTimeStepOnly] ;
Print[ Lie_az, OnElementsOf Vol_Mag, File "res/Lie_az.pos" ];
Echo[ Str["k= PostProcessing.NbViews-1;",
"View[k].IntervalsType = 3;",
"View[0].NbIso = 30;"],
File "tmp.geo", LastTimeStepOnly] ;
Print[ Lie_b_at_point, OnPoint{ Px_target, Py_target, 0},
Format Table, File > "res/Lie_w.txt",
SendToServer Sprintf["Optimization/Results/dwdtau_%g",VelocityTag]];
}
}
}
Ccore/ccore_common.pro 0 → 100644
+ 62
− 0
View file @ f0cb6778
/*
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
*/
ParamIndex =
DefineNumber[0, Name "Optimization/Sweep on parameter",
Choices {0="None",
1="D: air gap length",
2="E: C-core thickness",
3="F: Core-coil gap"} ];
StepsInRange =
DefineNumber[20, Name "Optimization/Steps in range", Visible ParamIndex];
w =
DefineNumber[0, Name "Optimization/Results/w",
Graph "030000000000000000000000000000000000"];
dwdt =
DefineNumber[0, Name Sprintf["Optimization/Results/dwdtau_%g",0],
Graph "000300000000000000000000000000000000"];
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]
];
DefineConstant[
// all dimensions in meter
L=1 // sidelength of the air box square
CoreX=0.3 // X-position of C-core's bottom left corner
CoreY=0.3 // Y-position of C-core's bottom left corner
A=0.4 // C-core length along X-axis
B=0.4 // C-core length along Y-axis
D=0.01 // air gap length
E=0.1 // C-core thickness
F=0.01 // gap between C-core and coil
G=0.05 // coil width along X-axis
CoilSection_ref = G*(B-2*E-2*F)
R=1 // refinement factor: R=5 far raw mesh, R=1 normal mesh
];
Ccore/shp.plt 0 → 100644
+ 53
− 0
View file @ f0cb6778
#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
Ccore/shp.py 0 → 100644
+ 245
− 0
View file @ f0cb6778
# Open this file with Gmsh (interactively with File->Open->shape.py, or on the
# command line with 'gmsh shape.py')
#
import numpy as np
import conveks
import onelab
from scipy.optimize import minimize
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))
modelName = 'ccore'
# 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)
# define optimization parameters as Onelab parameter (editable in the GUI)
maxIter = c.defineNumber('Optimization/00Max iterations', value=100)
maxChange = c.defineNumber('Optimization/01Max change', value=1e-0)
# 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 = [c.getNumber('Optimization/Results/dwdtau_0'), c.getNumber('Optimization/Results/dwdtau_1')]
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
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
# x is the 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:
x[0] = designParameters[index];
test(int(c.getNumber('Optimization/Steps in range')))
#callback=callbackF,
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:
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))]
## print x
## print values, bounds
## print fobj(values)
## print grad(values)
## exit(1)
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
exit(1)
# number of design variables and number of constraints
numVariables = len(x.keys())
# initial value of the variables in the design space,
# the lower bound for design variables,
# as well as the upper bound for design variables.
initialPoint = np.zeros(numVariables)
lowerBound = np.zeros(numVariables)
upperBound = np.zeros(numVariables)
for label, var in x.iteritems():
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-6 :
c.setNumber('Optimization/01Current iteration', value=it, readOnly=1,
attributes={'Highlight':'LightYellow'})
# get (copy of) current point
xFromMMA = np.array(conveks.mma.getCurrentPoint())
print 'xFromMMA = ', xFromMMA
# get the value of the objective function and of the constraints
# as well as their sensitivity with respect to design variables at `x'
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, 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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment