Commit e4c8f816 authored by Christophe Geuzaine's avatar Christophe Geuzaine

use explicit nonlinear loop in template

parent d32540fd
Pipeline #1469 passed with stage
in 9 minutes 28 seconds
......@@ -921,6 +921,30 @@ void Cal_GalerkinTermOfFemEquation(struct Element * Element,
}
}
// store unassembled RHS in DofData, per element
#if 0
for (j = 0 ; j < Nbr_Dof ; j++){
if(QuantityStorageDof_P->BasisFunction[j].Dof->Type == DOF_FIXED &&
QuantityStorageDof_P->BasisFunction[j].Dof->Entity == 0){
std::vector<std::pair<int, double> > &vec =
Current.DofData->unassembledRHS[Element->Num];
printf("ele %d: ", Element->Num);
for (i = 0 ; i < Nbr_Equ ; i++) {
if(QuantityStorageEqu_P->BasisFunction[i].Dof->Type == DOF_UNKNOWN){
for(k = 0 ; k < Current.NbrHar ; k++){
int n = QuantityStorageEqu_P->BasisFunction[i].Dof->Case.Unknown.NumDof;
printf("(%d, %g) ", n, Ek[i][j][k]) ;
vec.push_back(std::pair<int, double>(n, Ek[i][j][k]);
}
}
}
printf("\n") ;
}
}
#endif
/* Assembly in global matrix */
/* ------------------------- */
if (!Current.flagAssDiag) /*+++prov*/
......
......@@ -3550,11 +3550,19 @@ void Treatment_Operation(struct Resolution * Resolution_P,
}
}
break ;
case OPERATION_BREAK :
Flag_Break = 1;
break ;
/*
case OPERATION_GETSIZE:
{
Init_OperationOnSystem("GetSize",
Resolution_P, Operation_P, DofData_P0, GeoData_P0,
&DefineSystem_P, &DofData_P, Resolution2_P) ;
Cal_ZeroValue(&Value);
Value.Type = SCALAR;
Value.Val[0] = DofData_P->NbrDof;
Cal_StoreInVariable(&Value, Operation_P->Case.GetSize.VariableName);
}
break;
*/
case OPERATION_SLEEP :
Get_ValueOfExpressionByIndex(Operation_P->Case.Sleep.ExpressionIndex,
NULL, 0., 0., 0., &Value) ;
......
......@@ -25,16 +25,14 @@ Group {
Flag_NL = { 0, Choices{0,1}, Name "Parameters/Materials/1Nonlinear BH-curve"}
];
// These are the generic group names that are used in "Magnetostatics.pro"
Vol_Inf_Mag = Region[ AirInf ] ;
Vol_M_Mag = Region[ Magnet ] ;
Vol_L_Mag = Region[ {Air, AirInf, AirGap} ] ;
Vol_NL_Mag = Region[ {} ] ;
// These are the generic group names that are used in the
// "Lib_MagSta_a_phi.pro" template included below
Vol_Mag = Region[ {Air, AirInf, Core, AirGap, Magnet} ] ;
If(Flag_NL)
Vol_NL_Mag += Region[ {Core} ] ;
Else
Vol_L_Mag += Region[ {Core} ] ;
Vol_NL_Mag = Region[ {Core} ] ;
EndIf
Vol_Inf_Mag = Region[ AirInf ] ;
Vol_M_Mag = Region[ Magnet ] ;
}
Function {
......@@ -44,20 +42,18 @@ Function {
Name "Parameters/Materials/Core relative permeability"} ];
nu [ Region[{Air, AirInf, AirGap, Magnet}] ] = 1. / mu0;
mu [ Region[{Air, AirInf, AirGap, Magnet}] ] = mu0 ;
If(!Flag_NL)
nu [ Core ] = 1. / (murCore * mu0) ;
mu [ Core ] = murCore * mu0;
EndIf
If(Flag_NL)
Else
nu [ Core ] = SteelGeneric_nu[$1] ;
dhdb_NL [ Core ] = SteelGeneric_dhdb_NL[$1];
dhdb [ Core ] = SteelGeneric_dhdb[$1];
mu [ Core ] = SteelGeneric_mu[$1] ;
dbdh_NL [ Core ] = SteelGeneric_dbdh_NL[$1];
dbdh [ Core ] = SteelGeneric_dbdh[$1];
EndIf
mu [ Region[{Air, AirInf, AirGap, Magnet}] ] = mu0 ;
DefineConstant[ Hc = {920000,
Name "Parameters/Materials/hc", Label "Magnet coercive field (A/m)"} ];
hc [ Magnet ] = Rotate[ Vector[Hc, 0, 0.], 0, 0, Pi/2] ;
......
......@@ -19,10 +19,15 @@ DefineConstant[
Group {
DefineGroup[
Vol_Ele, // dielectric domain
// Full dielectric domain:
Vol_Ele,
// The following are subsets of Vol_Ele:
Vol_Q_Ele, // domain with imposed volume charge density
Vol_Inf_Ele, // infinite region
Sur_Neu_Ele, // Non-homogeneous Neumann boundary conditions (n.d)
// Boundaries:
Sur_Neu_Ele, // non-homogeneous Neumann boundary conditions (n.d)
Sur_C_Ele // boundary of conductors
];
Dom_Ele = Region[ {Vol_Ele, Sur_Neu_Ele} ];
......
......@@ -9,30 +9,33 @@
DefineConstant[
modelPath = "", // default path of the model
resPath = StrCat[modelPath, "res/"], // path for post-operation files
mu0 = 4*Pi*1e-7, // magnetic permeability of vacuum
Flag_NewtonRaphson = 1, // Newton-Raphson or Picard method for nonlinear iterations
modelDim = 2, // default model dimension (2D)
Val_Rint = 0, // internal radius of Vol_Inf_Ele annulus
Val_Rext = 0, // external radius of Vol_Inf_Ele annulus
Val_Cx = 0, // x-coordinate of center of Vol_Inf_Ele
Val_Cy = 0, // y-coordinate of center of Vol_Inf_Ele
Val_Cz = 0, // z-coordinate of center of Vol_Inf_Ele
NL_max_iter = 20, // max number of nonlinear iterations
NL_relax = 1, // nonlinear iteration relexation
NL_tol = 1e-5 // stopping criterion for nonlinear iterations
NL_tol_abs = 1e-6, // absolute tolerance on residual for noninear iterations
NL_tol_rel = 1e-6, // relative tolerance on residual for noninear iterations
NL_iter_max = 20 // maximum number of noninear iterations
];
Group {
DefineGroup[
Vol_L_Mag, // linear magnetic materials
// The full magnetic domain:
Vol_Mag,
// Subsets of Vol_Mag:
Vol_NL_Mag, // nonlinear magnetic materials
Vol_M_Mag, // permenent magnets
Vol_S0_Mag, // imposed current density
Vol_Inf_Mag, // infinite domains
// Boundaries:
Sur_Dir_Mag // Dirichlet boundary conditions
Sur_Neu_Mag // Non-homogeneous Neumann boundary conditions
];
Vol_Mag = Region[{Vol_L_Mag, Vol_NL_Mag, Vol_M_Mag, Vol_S0_Mag, Vol_Inf_Mag}];
Dom_Mag = Region[{Vol_Mag, Sur_Neu_Mag}];
}
Function{
......@@ -41,8 +44,8 @@ Function{
nu, // magnetic reluctivity (= 1/nu)
hc, // coercive magnetic field (in magnets)
js, // source current density
dhdb_NL, // nonlinar part of Jacobian for phi formulation
dbdh_NL, // nonlinear part of Jacobian for a formulation
dhdb, // Jacobian for Newton-Raphson method a formulation
dbdh, // Jacobian for Newton-Raphson method phi formulation
bn, // normal magnetic flux density on Sur_Neu_Mag for phi formulation
nxh // tangential magnetic field on Sur_Neu_Mag for a formulation
];
......@@ -50,6 +53,13 @@ Function{
// End of default definitions.
Group {
// linear materials
Vol_L_Mag = Region[{Vol_Mag, -Vol_NL_Mag}];
// all volumes + surfaces on which integrals must be computed
Dom_Mag = Region[{Vol_Mag, Sur_Neu_Mag}];
}
Jacobian {
{ Name Vol;
Case {
......@@ -66,7 +76,7 @@ Jacobian {
}
Integration {
{ Name I1;
{ Name Int;
Case {
{ Type Gauss;
Case {
......@@ -133,14 +143,26 @@ Formulation {
{ Name phi; Type Local; NameOfSpace Hgrad_phi; }
}
Equation {
Integral { [ - mu[-{d phi}] * Dof{d phi} , {d phi} ];
In Vol_Mag; Jacobian Vol; Integration I1; }
Integral { JacNL [ - dbdh_NL[-{d phi}] * Dof{d phi} , {d phi} ];
In Vol_NL_Mag; Jacobian Vol; Integration I1; }
Integral { [ - mu[] * Dof{d phi} , {d phi} ];
In Vol_L_Mag; Jacobian Vol; Integration Int; }
If(Flag_NewtonRaphson)
Integral { [ - mu[-{d phi}] * {d phi} , {d phi} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Integral { [ - dbdh[-{d phi}] * Dof{d phi} , {d phi}];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Integral { [ dbdh[-{d phi}] * {d phi} , {d phi}];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Else
Integral { [ - mu[-{d phi}] * Dof{d phi} , {d phi} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
EndIf
Integral { [ - mu[] * hc[] , {d phi} ];
In Vol_M_Mag; Jacobian Vol; Integration I1; }
In Vol_M_Mag; Jacobian Vol; Integration Int; }
Integral { [ bn[] , {phi} ];
In Sur_Neu_Mag; Jacobian Sur; Integration I1; }
In Sur_Neu_Mag; Jacobian Sur; Integration Int; }
}
}
{ Name MagSta_a; Type FemEquation;
......@@ -148,16 +170,29 @@ Formulation {
{ Name a; Type Local; NameOfSpace Hcurl_a; }
}
Equation {
Integral { [ nu[{d a}] * Dof{d a} , {d a} ];
In Vol_Mag; Jacobian Vol; Integration I1; }
Integral { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ];
In Vol_NL_Mag; Jacobian Vol; Integration I1; }
Integral { [ nu[] * Dof{d a} , {d a} ];
In Vol_L_Mag; Jacobian Vol; Integration Int; }
If(Flag_NewtonRaphson)
Integral { [ nu[{d a}] * {d a} , {d a} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Integral { [ dhdb[{d a}] * Dof{d a} , {d a} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Integral { [ - dhdb[{d a}] * {d a} , {d a} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
Else
Integral { [ nu[{d a}] * Dof{d a}, {d a} ];
In Vol_NL_Mag; Jacobian Vol; Integration Int; }
EndIf
Integral { [ hc[] , {d a} ];
In Vol_M_Mag; Jacobian Vol; Integration I1; }
Integral { [ -js[] , {a} ];
In Vol_S0_Mag; Jacobian Vol; Integration I1; }
In Vol_M_Mag; Jacobian Vol; Integration Int; }
Integral { [ - js[] , {a} ];
In Vol_S0_Mag; Jacobian Vol; Integration Int; }
Integral { [ nxh[] , {a} ];
In Sur_Neu_Mag; Jacobian Sur; Integration I1; }
In Sur_Neu_Mag; Jacobian Sur; Integration Int; }
}
}
}
......@@ -168,11 +203,19 @@ Resolution {
{ Name A; NameOfFormulation MagSta_phi; }
}
Operation {
If(!NbrRegions[Vol_NL_Mag])
Generate[A]; Solve[A];
Else
IterativeLoop[NL_max_iter, NL_tol, NL_relax]{
GenerateJac[A]; SolveJac[A];
InitSolution[A];
Generate[A]; Solve[A];
If(NbrRegions[Vol_NL_Mag])
Generate[A]; GetResidual[A, $res0];
Evaluate[ $res = $res0, $iter = 0 ];
Print[{$iter, $res, $res / $res0},
Format "Residual %03g: abs %14.12e rel %14.12e"];
While[$res > NL_tol_abs && $res / $res0 > NL_tol_rel &&
$res / $res0 <= 1 && $iter < NL_iter_max]{
Solve[A]; Generate[A]; GetResidual[A, $res];
Evaluate[ $iter = $iter + 1 ];
Print[{$iter, $res, $res / $res0},
Format "Residual %03g: abs %14.12e rel %14.12e"];
}
EndIf
SaveSolution[A];
......@@ -183,11 +226,19 @@ Resolution {
{ Name A; NameOfFormulation MagSta_a; }
}
Operation {
If(!NbrRegions[Vol_NL_Mag])
Generate[A]; Solve[A];
Else
IterativeLoop[NL_max_iter, NL_tol, NL_relax]{
GenerateJac[A]; SolveJac[A];
InitSolution[A];
Generate[A]; Solve[A];
If(NbrRegions[Vol_NL_Mag])
Generate[A]; GetResidual[A, $res0];
Evaluate[ $res = $res0, $iter = 0 ];
Print[{$iter, $res, $res / $res0},
Format "Residual %03g: abs %14.12e rel %14.12e"];
While[$res > NL_tol_abs && $res / $res0 > NL_tol_rel &&
$res / $res0 <= 1 && $iter < NL_iter_max]{
Solve[A]; Generate[A]; GetResidual[A, $res];
Evaluate[ $iter = $iter + 1 ];
Print[{$iter, $res, $res / $res0},
Format "Residual %03g: abs %14.12e rel %14.12e"];
}
EndIf
SaveSolution[A];
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment