diff --git a/Team32/BH_team32materialdata.pro b/Team32/BH_team32materialdata.pro
new file mode 100755
index 0000000000000000000000000000000000000000..714c4c36e328f115c368a62ee0878c501df8cab0
--- /dev/null
+++ b/Team32/BH_team32materialdata.pro
@@ -0,0 +1,57 @@
+Function{
+
+ /*
+ // The version NON-any below as a form like this (harder to deal with):
+ // (form like this: ___ )
+ // ( __| )
+ // ( __| )
+ */
+
+/*List_H={0,42.9013,51.4518,73.9502,108.958,202.746,251.765,348.401,492.594,727.547,906.266,1232.82,5441.59,1e+09};
+List_H2={0,1840.52,2647.29,5468.63,11871.8,41105.8,63385.5,121383,242649,529325,821319,1.51984e+06,2.9611e+07,1e+18};
+List_Janhy={0,0.706605,0.750429,0.834449,0.927388,1.09949,1.16511,1.26115,1.34964,1.4243,1.45561,1.48941,1.56197,1.58323};
+List_Banhy={0,0.706659,0.750493,0.834542,0.927525,1.09974,1.16542,1.26158,1.35025,1.42521,1.45675,1.49096,1.56881,1258.22};
+List_Banhy2={0,0.499366,0.56324,0.696461,0.860303,1.20943,1.35821,1.59159,1.82319,2.03123,2.12212,2.22296,2.46117,1.58312e+06};
+List_Muanhy={0.0479842,0.0164717,0.0145863,0.0112852,0.00851268,0.00542425,0.00462902,0.00362107,0.00274111,0.00195893,0.00160742,0.00120939,0.0002883,1.25822e-06};
+List_Nuanhy={20.8402,60.7101,68.5573,88.6116,117.472,184.357,216.029,276.161,364.816,510.484,622.115,826.861,3468.61,794773};
+List_J={0,0.150626,0.249049,0.499732,0.750528,1.0153,1.095,1.2004,1.29992,1.39605,1.44988,1.4999,1.56197,1.58323};
+List_B={0,0.15068,0.249113,0.499825,0.750665,1.01556,1.09531,1.20084,1.30054,1.39697,1.45102,1.50145,1.56881,1258.22};
+List_B2={0,0.0227043,0.0620575,0.249825,0.563498,1.03136,1.19971,1.44202,1.6914,1.95152,2.10545,2.25434,2.46117,1.58312e+06};
+List_Mu={0.000572824,0.00351224,0.00484168,0.00675895,0.00688949,0.00500902,0.00435054,0.00344672,0.00264018,0.00192011,0.00160109,0.0012179,0.0002883,1.25822e-06};
+List_Nu={1745.74,284.719,206.54,147.952,145.149,199.64,229.857,290.131,378.762,520.804,624.574,821.086,3468.61,794773};
+*/
+
+// More points
+List_H={0,1,5.06122,9.12245,13.1837,17.2449,21.3061,25.3673,29.4286,33.4898,37.551,41.6122,42.9013,45.6735,49.7347,51.4518,53.7959,57.8571,61.9184,65.9796,70.0408,73.9502,74.102,78.1633,82.2245,86.2857,90.3469,94.4082,98.4694,102.531,106.592,108.958,110.653,114.714,118.776,122.837,126.898,130.959,135.02,139.082,143.143,147.204,151.265,155.327,159.388,163.449,167.51,171.571,175.633,179.694,183.755,187.816,191.878,195.939,200,202.746,251.765,348.401,492.594,727.547,906.266,1232.82,5441.59,1e+09};
+List_H2={0,1,25.616,83.2191,173.809,297.387,453.951,643.502,866.041,1121.57,1410.08,1731.58,1840.52,2086.07,2473.54,2647.29,2894,3347.45,3833.88,4353.31,4905.72,5468.63,5491.11,6109.5,6760.87,7445.22,8162.57,8912.9,9696.22,10512.5,11361.8,11871.8,12244.1,13159.4,14107.6,15088.9,16103.1,17150.3,18230.5,19343.7,20489.9,21669,22881.2,24126.3,25404.5,26715.6,28059.7,29436.8,30846.8,32289.9,33765.9,35275,36817,38392,40000,41105.8,63385.5,121383,242649,529325,821319,1.51984e+06,2.9611e+07,1e+18};
+List_Janhy={0,0.0309761,0.153936,0.266504,0.363752,0.444473,0.510097,0.563225,0.606585,0.642523,0.672877,0.699018,0.706605,0.721948,0.742395,0.750429,0.760892,0.777831,0.793504,0.808137,0.821899,0.834449,0.834925,0.84732,0.859169,0.87054,0.88149,0.892064,0.902301,0.912232,0.921884,0.927388,0.93128,0.94044,0.94938,0.958114,0.966655,0.975014,0.9832,0.991222,0.999086,1.0068,1.01437,1.0218,1.02909,1.03626,1.0433,1.05021,1.05701,1.06368,1.07025,1.0767,1.08304,1.08927,1.0954,1.09949,1.16511,1.26115,1.34964,1.4243,1.45561,1.48941,1.56197,1.58323};
+List_Banhy={0,0.0309774,0.153942,0.266515,0.363768,0.444494,0.510123,0.563257,0.606622,0.642565,0.672924,0.69907,0.706659,0.722005,0.742458,0.750493,0.76096,0.777903,0.793582,0.80822,0.821987,0.834542,0.835018,0.847418,0.859272,0.870649,0.881604,0.892183,0.902425,0.912361,0.922018,0.927525,0.931419,0.940584,0.949529,0.958268,0.966814,0.975178,0.98337,0.991397,0.999266,1.00699,1.01456,1.02199,1.02929,1.03646,1.04351,1.05043,1.05723,1.06391,1.07048,1.07693,1.08328,1.08952,1.09565,1.09974,1.16542,1.26158,1.35025,1.42521,1.45675,1.49096,1.56881,1258.22};
+List_Banhy2={0,0.000959599,0.0236982,0.0710303,0.132327,0.197575,0.260226,0.317259,0.367991,0.41289,0.452827,0.488699,0.499366,0.521292,0.551244,0.56324,0.57906,0.605133,0.629773,0.653219,0.675662,0.696461,0.697254,0.718117,0.738349,0.758029,0.777225,0.795991,0.81437,0.832402,0.850117,0.860303,0.867542,0.884698,0.901605,0.918278,0.93473,0.950973,0.967016,0.982867,0.998533,1.01402,1.02933,1.04447,1.05945,1.07426,1.08891,1.1034,1.11773,1.1319,1.14592,1.15978,1.17349,1.18705,1.20046,1.20943,1.35821,1.59159,1.82319,2.03123,2.12212,2.22296,2.46117,1.58312e+06};
+List_Muanhy={0.0309975,0.0309774,0.030416,0.0292153,0.0275923,0.0257754,0.0239426,0.022204,0.0206134,0.0191869,0.0179203,0.0167996,0.0164717,0.015808,0.0149284,0.0145863,0.0141453,0.0134452,0.0128166,0.0122495,0.0117358,0.0112852,0.0112685,0.0108416,0.0104503,0.0100903,0.00975798,0.00945027,0.00916452,0.00889842,0.00864999,0.00851268,0.00841747,0.00819936,0.00799431,0.00780115,0.00761883,0.00744643,0.00728312,0.00712816,0.0069809,0.00684074,0.00670716,0.00657965,0.0064578,0.00634121,0.00622951,0.00612239,0.00601953,0.00592067,0.00582556,0.00573396,0.00564568,0.00556051,0.00547827,0.00542425,0.00462902,0.00362107,0.00274111,0.00195893,0.00160742,0.00120939,0.0002883,1.25822e-06};
+List_Nuanhy={32.2607,32.2816,32.8775,34.2286,36.242,38.7967,41.7666,45.0369,48.5122,52.1189,55.8027,59.5251,60.7101,63.2592,66.9865,68.5573,70.6948,74.3758,78.0239,81.6357,85.2092,88.6116,88.7431,92.237,95.6909,99.1051,102.48,105.817,109.116,112.379,115.607,117.472,118.801,121.961,125.089,128.186,131.254,134.293,137.304,140.289,143.248,146.183,149.094,151.984,154.851,157.699,160.526,163.335,166.126,168.9,171.657,174.399,177.127,179.84,182.539,184.357,216.029,276.161,364.816,510.484,622.115,826.861,3468.61,794773};
+List_J={0,0.000102733,0.00257973,0.00821228,0.0168,0.0281424,0.0420391,0.0582897,0.0766939,0.0970511,0.119161,0.142823,0.150626,0.175788,0.227968,0.249049,0.27606,0.32408,0.372237,0.418886,0.46238,0.499732,0.501078,0.536706,0.571417,0.60484,0.636603,0.666336,0.693668,0.718228,0.739646,0.750528,0.757812,0.774727,0.790909,0.80639,0.821202,0.835378,0.848948,0.861945,0.874401,0.886347,0.897815,0.908838,0.919446,0.929673,0.93955,0.949108,0.95838,0.967397,0.976192,0.984796,0.993241,1.00156,1.00978,1.0153,1.095,1.2004,1.29992,1.39605,1.44988,1.4999,1.56197,1.58323};
+List_B={0,0.00010399,0.00258609,0.00822374,0.0168165,0.028164,0.0420659,0.0583216,0.0767308,0.0970931,0.119208,0.142875,0.15068,0.175845,0.228031,0.249113,0.276127,0.324153,0.372315,0.418969,0.462468,0.499825,0.501171,0.536805,0.571521,0.604948,0.636716,0.666454,0.693791,0.718357,0.73978,0.750665,0.757951,0.774871,0.791058,0.806545,0.821362,0.835542,0.849118,0.86212,0.87458,0.886532,0.898005,0.909033,0.919647,0.929879,0.93976,0.949324,0.958601,0.967623,0.976423,0.985032,0.993482,1.00181,1.01003,1.01556,1.09531,1.20084,1.30054,1.39697,1.45102,1.50145,1.56881,1258.22};
+List_B2={0,1.08138e-08,6.68786e-06,6.76299e-05,0.000282796,0.000793213,0.00176954,0.00340141,0.00588762,0.00942708,0.0142106,0.0204134,0.0227043,0.0309215,0.0519981,0.0620575,0.0762462,0.105075,0.138619,0.175535,0.213877,0.249825,0.251172,0.288159,0.326636,0.365962,0.405408,0.444161,0.481346,0.516036,0.547274,0.563498,0.57449,0.600425,0.625773,0.650514,0.674635,0.698131,0.721001,0.743251,0.764891,0.785938,0.806413,0.826341,0.84575,0.864674,0.883149,0.901215,0.918915,0.936295,0.953402,0.970288,0.987007,1.00361,1.02017,1.03136,1.19971,1.44202,1.6914,1.95152,2.10545,2.25434,2.46117,1.58312e+06};
+List_Mu={7.43501e-05,0.00010399,0.000510961,0.000901484,0.00127556,0.00163318,0.00197436,0.00229908,0.00260736,0.00289919,0.00317456,0.00343349,0.00351224,0.00385005,0.00458495,0.00484168,0.00513286,0.00560264,0.006013,0.00634998,0.00660284,0.00675895,0.00676325,0.00686774,0.00695074,0.00701099,0.00704746,0.00705929,0.00704576,0.00700627,0.0069403,0.00688949,0.0068498,0.00675479,0.00666011,0.00656599,0.00647262,0.00638017,0.00628881,0.00619866,0.00610984,0.00602247,0.00593662,0.0058524,0.00576987,0.00568911,0.00561017,0.00553311,0.00545799,0.00538484,0.00531372,0.00524466,0.00517769,0.00511285,0.00505017,0.00500902,0.00435054,0.00344672,0.00264018,0.00192011,0.00160109,0.0012179,0.0002883,1.25822e-06};
+List_Nu={13449.9,9616.35,1957.1,1109.28,783.971,612.302,506.494,434.956,383.53,344.924,315.004,291.249,284.719,259.737,218.105,206.54,194.823,178.487,166.306,157.481,151.45,147.952,147.858,145.608,143.87,142.633,141.895,141.657,141.929,142.729,144.086,145.149,145.99,148.043,150.148,152.3,154.497,156.736,159.013,161.325,163.67,166.045,168.446,170.87,173.314,175.775,178.248,180.73,183.218,185.706,188.192,190.67,193.136,195.586,198.013,199.64,229.857,290.131,378.762,520.804,624.574,821.086,3468.61,794773};
+
+
+any_nu_b2 = ListAlt[List_B2, List_Nu] ;
+//any_nu_b2 = ListAlt[List_Banhy2, List_Nuanhy] ; // take the any version (see header)
+
+// Akima Interpolation problematic with Nu ...
+nu_any[] = InterpolationAkima[$1]{List[any_nu_b2]} ;
+dnudb2_any[] = dInterpolationAkima[$1]{List[any_nu_b2]} ;
+
+// Try Linear Interpolation instead with Nu ...
+//nu_any[] = InterpolationLinear[$1]{List[any_nu_b2]} ;
+//dnudb2_any[] = dInterpolationLinear[$1]{List[any_nu_b2]} ;
+
+h_any[] = nu_any[(SquNorm[$1])] * $1 ;
+dhdb_any[] = TensorDiag[1,1,1] * nu_any[SquNorm[$1]] + 2*dnudb2_any[SquNorm[$1]] * SquDyadicProduct[$1] ;
+dhdb_any_NL[] = 2*dnudb2_any[SquNorm[$1]] * SquDyadicProduct[$1] ;
+
+//nu[DomainNL]=nu_any[(SquNorm[$1])]; // already in magstadyn_a.pro
+//dnudbxb[] = 2*dnudb2_any[(SquNorm[$1])] * SquDyadicProduct[$1]; // defined here abobe as dhdb_any_NL
+
+}
diff --git a/Team32/MyCircuit.pro b/Team32/MyCircuit.pro
index 269876a9180351426e47be5d95d7691675fe0b35..94d52a196c0e35f8077cda4984beadb067a46aa4 100644
--- a/Team32/MyCircuit.pro
+++ b/Team32/MyCircuit.pro
@@ -69,16 +69,16 @@ Constraint {
If (Flag_TestCase==3)
Case Circuit1 {
- { Region Ein1 ; Branch {1,2} ; }
+ { Region Ein1 ; Branch {2,1} ; } // 2 1 since 25/8/18 !! (E1Inv)
{ Region R1 ; Branch {2,3} ; }
{ Region IndB~{1} ; Branch {3,4} ; }
- { Region IndA~{1} ; Branch {4,1} ; }
+ { Region IndA~{1} ; Branch {1,4} ; } // 1 4 since 25/8/18 !!
}
Case Circuit2 {
- { Region Ein2 ; Branch {1,2} ; }
+ { Region Ein2 ; Branch {1,2} ; } //(E2dir)
{ Region R2 ; Branch {2,3} ; }
{ Region IndB~{2} ; Branch {3,4} ; }
- { Region IndA~{2} ; Branch {4,1} ; }
+ { Region IndA~{2} ; Branch {1,4} ; } // 4 1 since 25/8/18 !!
}
EndIf
If (Flag_TestCase==4)
@@ -91,7 +91,7 @@ Constraint {
Case Circuit2 {
{ Region R4 ; Branch {1,2} ; }
{ Region IndB~{2} ; Branch {2,3} ; }
- { Region IndA~{2} ; Branch {3,1} ; }
+ { Region IndA~{2} ; Branch {3,1} ; }
}
EndIf
EndIf
@@ -102,13 +102,13 @@ Constraint {
{ Region Ein1 ; Branch {1,2} ; }
{ Region R12 ; Branch {2,3} ; }
{ Region Ind~{1} ; Branch {3,4} ; }
- { Region Ind~{2} ; Branch {4,1} ; }
+ { Region Ind~{2} ; Branch {4,1} ; }
}
EndIf
If (Flag_TestCase==3)
Case Circuit1 {
- { Region Ein1 ; Branch {1,2} ; }
+ { Region Ein1 ; Branch {2,1} ; } // 2 1 since 12/12/18 !! (E1Inv)
{ Region R1 ; Branch {2,3} ; }
{ Region Ind~{1} ; Branch {3,1} ; }
}
diff --git a/Team32/MyItLoopPro.pro b/Team32/MyItLoopPro.pro
index c2b29e8cdb8b56ec86e0c72c9cd79d83b8b6a630..f95639a861cb8d7f50de65b567305e2e8c528dbe 100644
--- a/Team32/MyItLoopPro.pro
+++ b/Team32/MyItLoopPro.pro
@@ -42,8 +42,9 @@ Macro MyItLoopPro
Test[Flag_AdaptRelax==1]{
CopySolution[A,'x_Prev']; //SET : save solution from the previous converged time step at iter 0
- Evaluate[ $res_Prev=$res,
- $relax_Opt=$relax];
+ Evaluate[ $res_Opt=$res,
+ $relax_Opt=$relax,
+ $inc_Opt=0];
}
{
Generate[A];
@@ -60,7 +61,12 @@ Macro MyItLoopPro
// if adapted relaxation is activated
Test[Flag_AdaptRelax==1]
{
- Call MySolveJac_AdaptRelax;
+ If(Flag_AdaptRelax_Version ==1 )
+ Call MySolveJac_AdaptRelax;
+ EndIf
+ If(Flag_AdaptRelax_Version ==2 )
+ Call MySolveJac_AdaptRelax2;
+ EndIf
}
// ... otherwise, if adapted relaxation is not activated
@@ -69,12 +75,13 @@ Macro MyItLoopPro
Evaluate[ $syscount = $syscount + 1,
$relaxcount=1,
$relaxcounttot=$relaxcounttot+1 ];
- Generate[A]; GetResidual[A, $res];
+ Generate[A];
+ GetNormIncrement[A, $inc]; //Check the new increment
// TODO: GetIncrement h_only
// " PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} "
+ GetResidual[A, $res];
+ Test[$iter==1]{Evaluate[$res0=$res];} // compute first residual at step $TimeStep
}
-
- GetNormIncrement[A, $inc]; //Check the new increment
GetNormSolution[A, $sol]; //Check the new solution
Call DisplayRunTimeVar;
@@ -90,51 +97,55 @@ Macro MySolveJac_AdaptRelax2
//NB:solution from the previous converged time step at iter 0 is saved in x_Prev at this point
Evaluate[ $redo=1,
$relax=relax_max,
- $res_Prev=$res,
+ $res_Opt=$res,
$relaxcount=0];
/* get increment dx by solving the system */
Generate[A];
Solve[A];
Evaluate[ $syscount = $syscount + 1];
- CopySolution[A,'x_New']; //SET : save increment dx (found without relaxation $relax=1)
- //TODO: "DX = x_New - x_Prev"
- // or TODO: CopyIncrement[A,'DX'];
+ CopySolution[A,'x_New'];
+ //SET : save increment dx (found without relaxation $relax=relax_max)
+ AddVector[A, 1/relax_max, 'x_New', -1/relax_max, 'x_Prev', 'Delta_x'];
While[ $relax>relax_min-0.5*relax_step &&
$redo>0 ]
// while all relaxation factor have not been tested yet
{
/* new trial solution = x + Frelax * dx */
- // TODO:"x_New = x_Prev + $relax * DX"
+ AddVector[A, 1, 'x_Prev', $relax, 'Delta_x', 'x_New'];
CopySolution['x_New', A]; // RESET: reset new trial solution in the system
-
/* calculate residual with trial solution res= b-A(x_New)*x_New*/
Generate[A]; // regenerate A with new trial solution (x_New)
- GetResidual[A, $res]; //Check the new residual
// TODO: GetIncrement h_only
// " PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} "
+ AddVector[A, 1, 'x_New', -1, 'x_Prev', 'weighted_Delta_x'];
+ CopyIncrement['weighted_Delta_x',A];
+ GetNormIncrement[A, $inc]; //Check the new increment
+ GetResidual[A, $res]; //Check the new residual
Evaluate[ $relaxcount=$relaxcount+1,
$relaxcounttot=$relaxcounttot+1 ];
// compute first residual at step $TimeStep
Test[$iter==1 && $relax==relax_max]{
Evaluate[ $res0=$res,
- $res_Prev=$res,
- $relax_Opt=$relax ];
+ $res_Opt=$res,
+ $relax_Opt=$relax,
+ $inc_Opt=$inc ];
CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one
}
/* // If one wants to be more verbose
- Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Prev ), $relaxcount}, Format
- "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.1g, (better? %g), relaxcount=%02g"];
+ Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Opt ), $relaxcount}, Format
+ "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.5g, (better? %g), relaxcount=%02g (V2)"];
//*/
// if the residual has decreased
- Test[ ($res < $res_Prev ) || $relax>relax_max-0.5*relax_step ]{
+ Test[ ($res < $res_Opt ) || $relax==relax_max ]{
CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one 'x_New is best'
- Evaluate[$res_Prev=$res,
- $relax_Opt=$relax];
+ Evaluate[$res_Opt=$res,
+ $relax_Opt=$relax,
+ $inc_Opt=$inc];
}
{ // otherwise, if the residual has not decreased...
// and if you don't want to check All Factors: Break the try loop prematurely
@@ -145,8 +156,9 @@ Macro MySolveJac_AdaptRelax2
}
CopySolution['x_Opt',A]; //RESET : reload optimal solution
CopySolution[A,'x_Prev']; //SET : save solution from the previous converged iteration
- Evaluate[$relax=$relax_Opt,
- $res=$res_Prev];
+ Evaluate[$res=$res_Opt,
+ $relax=$relax_Opt,
+ $inc=$inc_Opt];
Return
//...............................................................
@@ -158,7 +170,7 @@ Macro MySolveJac_AdaptRelax
Evaluate[ $redo=1,
$relax=relax_max,
- $res_Prev=$res,
+ $res_Opt=$res,
$relaxcount=0];
While[ $relax>relax_min-0.5*relax_step &&
@@ -172,26 +184,29 @@ Macro MySolveJac_AdaptRelax
$relaxcount=$relaxcount+1,
$relaxcounttot=$relaxcounttot+1 ];
Generate[A];
+ GetNormIncrement[A, $inc]; //Check the new increment
GetResidual[A, $res]; //Check the new residual
// compute first residual at step $TimeStep
Test[$iter==1 && $relax==relax_max]{
Evaluate[ $res0=$res,
- $res_Prev=$res,
- $relax_Opt=$relax ];
+ $res_Opt=$res,
+ $relax_Opt=$relax,
+ $inc_Opt=$inc ];
CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one}
}
/* // If one wants to be more verbose
- Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Prev ), $relaxcount}, Format
- "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.1g, (better? %g), relaxcount=%02g"];
+ Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, ($res < $res_Opt ), $relaxcount}, Format
+ "ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relax = %1.5g, (better? %g), relaxcount=%02g (V1)"];
//*/
// if the residual has decreased
- Test[ ($res < $res_Prev ) || $relax>relax_max-0.5*relax_step ]{
+ Test[ ($res < $res_Opt ) || $relax==relax_max ]{
CopySolution[A,'x_Opt']; //SET : save actual solution as the optimal one
- Evaluate[$res_Prev=$res,
- $relax_Opt=$relax];
+ Evaluate[$res_Opt=$res,
+ $relax_Opt=$relax,
+ $inc_Opt=$inc];
}
{ // otherwise, if the residual has not decreased...
// and if you don't want to check All Factors: Break the try loop prematurely
@@ -202,8 +217,9 @@ Macro MySolveJac_AdaptRelax
}
CopySolution['x_Opt',A]; //RESET : reload optimal solution
CopySolution[A,'x_Prev']; //SET : save solution from the previous converged iteration
- Evaluate[$relax=$relax_Opt,
- $res=$res_Prev];
+ Evaluate[$res=$res_Opt,
+ $relax=$relax_Opt,
+ $inc=$inc_Opt];
Return
//...............................................................
@@ -215,12 +231,12 @@ Macro DisplayRunTimeVar
If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
Print[{$TimeStep, $iter, $res, $resL, $resN, $relax, $relaxcount}, Format
- "*** ts=%05g, it=%03g, res_used=%1.5e, (res = %1.5e, resN=%1.5e), relaxopt = %1.3g, relaxcount=%02g ***"];
+ "*** ts=%05g, it=%03g, res_used=%1.5e, (res = %1.5e, resN=%1.5e), relaxopt = %1.5g, relaxcount=%02g ***"];
EndIf
If(Flag_NLRes==NLRES_ITERATIVELOOPPRO)
Print[{$TimeStep, $iter, $res, $res/$res0, $inc, $relax, $relaxcount}, Format
- "*** ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relaxopt = %1.3g, relaxcount=%02g ***"];
+ "*** ts=%05g, it=%03g, abs=%1.5e, rel=%1.5e, dx=%1.5e, relaxopt = %1.5g, relaxcount=%02g ***"];
EndIf
PostOperation[PostOpItLoop_a];
diff --git a/Team32/magstadyn_a.pro b/Team32/magstadyn_a.pro
index 6aaeb800be98f9f195ce559d6f2a8b407582424f..d968275e4a51fc06c7b70689587e8b39f62ef49a 100644
--- a/Team32/magstadyn_a.pro
+++ b/Team32/magstadyn_a.pro
@@ -46,7 +46,10 @@ Function{
Flag_PrintMaps =1,
Flag_LiveLocalPostOp=1,
Flag_LiveGlobalPostOp=1,
- po = "Output_"
+ po = "Output_",
+ Flag_ExtrapolationOrder=0,
+ Thickness=1,
+ Flag_DomainLam=0
];
// Macro myname
@@ -68,6 +71,7 @@ Printf("SymmetryFactor %g", SymmetryFactor);
Include "BH.pro"; // nonlinear BH caracteristic of magnetic material
Include "BH_anhysteretic.pro";
+Include "BH_team32materialdata.pro"; //New add
If(Flag_NL_law==NL_ENERGHYST)
Include "param_EnergHyst.pro";
Else
@@ -85,6 +89,7 @@ Group {
DomainCC = Region[ {Air, AirInf, Inds} ];
DomainC = Region[ { } ];
+ DomainLam = Region[ {} ];
If(!Flag_ConductingCore)
DomainCC += Region[ {Core} ];
@@ -100,6 +105,10 @@ Group {
DomainL = Region[ {Domain,-DomainNL} ];
EndIf
DomainDummy = #12345 ; // Dummy region number for postpro with functions
+
+ If(Flag_DomainLam)
+ DomainLam += Region[ {Core} ];
+ EndIf
}
Function {
@@ -115,11 +124,26 @@ Function {
dhdb [ DomainNL ] = dhdb_EIcore[$1];
*/
- If(Flag_NL_law==NL_ANA ) // Using analytical law
+ If(Flag_NL_law==NL_ANA ) // Using analytical law
+ // From BH.pro ...................................
h[DomainNL] = h_1a[$1]; // $1 => b ={d a}
dhdb[DomainNL] = dhdb_1a[$1];
nu[DomainNL] = nu_1a[$1] ;
dhdb_NL[DomainNL] = dhdb_1a_NL[$1];
+ /*
+ // From BH_anhysteretic.pro.......................
+ h[DomainNL] = h_anhys[$1];
+ dhdb[DomainNL] = dhdb_anhys[$1];
+ nu[DomainNL] = nu_anhys[(SquNorm[$1])];
+ dhdb_NL[DomainNL] = dhdb_anhys_NL[$1];
+ */
+ /*
+ // From BH_team32materialdata.pro..................
+ h[DomainNL] = h_any[$1];
+ dhdb[DomainNL] = dhdb_any[$1];
+ nu[DomainNL] = nu_any[(SquNorm[$1])];
+ dhdb_NL[DomainNL] = dhdb_any_NL[$1];
+ */
EndIf
If(Flag_NL_law==NL_ANA_JA) //Using anhysteretic curve from Jiles-Atherton model
h[DomainNL] = h_anhys[$1];
@@ -373,7 +397,7 @@ Formulation {
// EnergHyst++b
If(Flag_NL_law==NL_ENERGHYST)
For iSub In {1:N}
- { Name J~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
+ { Name X~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
EndFor
// EnergHyst++e
EndIf
@@ -381,17 +405,35 @@ Formulation {
Equation {
If(Flag_NL_law==NL_LIN || Flag_NL_law==NL_ANA || Flag_NL_law==NL_ANA_JA )
+ // nu-Version) with nu[{d a}] * Dof{d a} and dhdb_NL[...]
+ //-------------------------------------------------------
+ /*
Galerkin { [ nu[{d a}] * Dof{d a} , {d a} ] ;
In Domain ; Jacobian Vol ; Integration I1 ; }
If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d ap} ] ;
In DomainNL ; Jacobian Vol ; Integration I1 ; }
Else
- Galerkin { [ (1/$relax) * (SetVariable[(dhdb_NL[{d a}]),
- QuadraturePointIndex[]]{$dhdb_any_temp}) * Dof{d a} , {d a} ] ;
+ Galerkin { [ ((1/$relax) * (SetVariable[(dhdb_NL[{d a}]),
+ QuadraturePointIndex[]]{$dhdb_any_temp})) * Dof{d a} , {d ap} ] ;
In DomainNL ; Jacobian Vol ; Integration I1 ; }
- Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_any_temp}) * {d a} , {d a} ] ;
+ Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_any_temp}) * {d a} , {d ap} ] ;
In DomainNL ; Jacobian Vol ; Integration I1 ; }
+ /*/
+ // h-Version) with {h}=nu[{d a}]*{d a} and dhdb[...]
+ //--------------------------------------------------
+ Galerkin { [ nu[{d a}] * {d a} , {d a} ] ;
+ In Domain ; Jacobian Vol ; Integration I1 ; }
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL [ dhdb[{d a}] * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1 ; }
+ Else
+ Galerkin { [ ((1/$relax) * (SetVariable[(dhdb[{d a}]),
+ QuadraturePointIndex[]]{$dhdb_any_temp})) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1 ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_any_temp}) * {d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1 ; }
+ //*/
EndIf
EndIf
@@ -407,13 +449,27 @@ Formulation {
// BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
Galerkin { [ Dof{h} , {h} ] ;
In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Galerkin { [ -h_Jiles[ {h}[1] ,
- {d a}[1],
- {d a} ]{List[hyst_FeSi]} , {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
-
+
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[( h_Jiles[ {h}[1] ,
+ {d a}[1],
+ {d a} ]{List[hyst_FeSi]}),
+ QuadraturePointIndex[]]{$h_Jiles_temp}) , {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else
+ Galerkin { [ - (SetVariable[ $relax * ( SetVariable[ h_Jiles[ {h}[1] ,
+ {d a}[1],
+ {d a} ]{List[hyst_FeSi]}
+ ,QuadraturePointIndex[]]{$h_Jiles_unrelaxed}
+ )
+ +(1-$relax)* {h},
+ QuadraturePointIndex[]
+ ]{$h_Jiles_temp}) , {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
// NL part ...
- Galerkin { [ Dof{h}, {d a} ] ; // register of Dof => idem
+ //Galerkin { [ Dof{h}, {d a} ] ; // register of Dof => idem
+ Galerkin { [ (GetVariable[ QuadraturePointIndex[]]{$h_Jiles_temp}), {d a} ];
In DomainNL ; Jacobian Vol ; Integration I1p ; }
If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
Galerkin { JacNL[ dhdb_Jiles[ {h} ,
@@ -424,9 +480,9 @@ Formulation {
Galerkin { [ (1/$relax) * (SetVariable[(dhdb_Jiles[{h} ,
{d a} ,
{h}-{h}[1] ]{List[hyst_FeSi]}),
- QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * Dof{d a} , {d a} ];
+ QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * Dof{d a} , {d ap} ];
In DomainNL ; Jacobian Vol ; Integration I1p ; } // kj modif I1 --> I1p here
- Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * {d a} , {d a} ];
+ Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * {d a} , {d ap} ];
In DomainNL ; Jacobian Vol ; Integration I1p ; } // kj modif I1 --> I1p here
EndIf
EndIf
@@ -434,123 +490,284 @@ Formulation {
//*******************************************************
//--------Hysteresis law with EnergHyst model----------
//*******************************************************
- If(Flag_NL_law==NL_ENERGHYST && N!=5) //EnergHyst model law
+ If(Flag_NL_law==NL_ENERGHYST) //EnergHyst model law
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In DomainL ; Jacobian Vol ; Integration I1 ; }
- // update h
- // saving h_Vinch_K in local quantity h
- // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
- Galerkin { [ Dof{h} , {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- //Galerkin { [ -0*Vector[X[],Y[],Z[]] , {h} ] ; In DomainNL ; Jacobian Vol ; Integration I1p ; } //For Debug
- Galerkin { [ - (SetVariable[(h_Vinch_K[ {h} , // new since 2/3/2018: use {h} instead of {h}[1] here (need to init bc by recomputing b from {h} in h_Vinch_K)
- {d a}, {d a}[1],
- {J_1}, {J_1}[1],
- {J_2}, {J_2}[1],
- {J_3}, {J_3}[1] ]{List[param_EnergHyst]}),
- QuadraturePointIndex[]]{$h_Vinch_temp}), {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
-
- // update Jk
- // saving Update_Cell_K in local quantity Jk
- // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
- For iSub In {1:N}
- Galerkin { [ Dof{J~{iSub}} ,{J~{iSub}} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Galerkin { [ - Update_Cell_K[ iSub ,
- (GetVariable[QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J~{iSub}} ,
- {J~{iSub}}[1] ]{List[param_EnergHyst]}, {J~{iSub}} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- EndFor
-
- // NL part ...
- Galerkin { [ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}), {d a} ];
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
-
- If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
- Galerkin { JacNL[ (dhdb_Vinch_K[ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J_1},{J_1}[1],
- {J_2},{J_2}[1],
- {J_3},{J_3}[1]]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Else
- Galerkin { [ (1/$relax) * (SetVariable[(dhdb_Vinch_K[ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J_1},{J_1}[1],
- {J_2},{J_2}[1],
- {J_3},{J_3}[1]]{List[param_EnergHyst]}) ,
- QuadraturePointIndex[]]{$dhdb_Vinch_temp}) * Dof{d a} , {d a} ];
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ If (!Flag_UpdateSeparated)
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // U P D A T E - H
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // saving h_EB in local quantity h
+ // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
+ Galerkin { [ Dof{h} , {h} ] ;
In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_Vinch_temp}) * {d a} , {d a} ];
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- EndIf
- EndIf
- //*******************************************************
- If(Flag_NL_law==NL_ENERGHYST && N==5) //EnergHyst model law // FOR NOW {J_4} and {J_5} have to be added by hand
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1]
+ ]{List[param_EnergHyst]}),
+ QuadraturePointIndex[]]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - (SetVariable[ $relax* ( SetVariable[ h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1]
+ ]{List[param_EnergHyst]}
+ , QuadraturePointIndex[]
+ ]{$h_EB_temp_unrelaxed}
+ )
+ +(1-$relax)*{h},
+ QuadraturePointIndex[]
+ ]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
+ ElseIf (N==3)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}),
+ QuadraturePointIndex[]]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - (SetVariable[ $relax* ( SetVariable[ h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}
+ , QuadraturePointIndex[]
+ ]{$h_EB_temp_unrelaxed}
+ )
+ +(1-$relax)*{h},
+ QuadraturePointIndex[]
+ ]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
+ ElseIf (N==5)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}),
+ QuadraturePointIndex[]]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - (SetVariable[ $relax* ( SetVariable[ h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1], {X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}
+ , QuadraturePointIndex[]
+ ]{$h_EB_temp_unrelaxed}
+ )
+ +(1-$relax)*{h},
+ QuadraturePointIndex[]
+ ]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // U P D A T E - C E L L S
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // saving Cell_EB in local quantity Xk
+ // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
+ For iSub In {1:N}
+ Galerkin { [ Dof{X~{iSub}} ,{X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - Cell_EB[ iSub ,
+ (GetVariable[QuadraturePointIndex[]]{$h_EB_temp}),
+ {X~{iSub}}[1] ]{List[param_EnergHyst]}, {X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - ( $relax* Cell_EB[ iSub ,
+ (GetVariable[QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X~{iSub}}[1] ]{List[param_EnergHyst]}
+ +(1-$relax)*{X~{iSub}} ), {X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
- Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
- In DomainL ; Jacobian Vol ; Integration I1 ; }
+ EndFor
+
+ // NL part ...
+
+ // H FIELD >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ Galerkin { [ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}), {d a} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ // DHDB TENSOR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}),
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
- // update h
- // saving h_Vinch_K in local quantity h
- // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
- Galerkin { [ Dof{h} , {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Galerkin { [ - (SetVariable[(h_Vinch_K[ {h} , // new since 2/3/2018: use {h} instead of {h}[1] here (need to init bc by recomputing b from {h} in h_Vinch_K)
- {d a}, {d a}[1],
- {J_1}, {J_1}[1],
- {J_2}, {J_2}[1],
- {J_3}, {J_3}[1],
- {J_4}, {J_4}[1],
- {J_5}, {J_5}[1] ]{List[param_EnergHyst]}),
- QuadraturePointIndex[]]{$h_Vinch_temp}), {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
-
- // update Jk
- // saving Update_Cell_K in local quantity Jk
- // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
- For iSub In {1:N}
- Galerkin { [ Dof{J~{iSub}} ,{J~{iSub}} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Galerkin { [ - Update_Cell_K[ iSub ,
- (GetVariable[QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J~{iSub}} ,
- {J~{iSub}}[1] ]{List[param_EnergHyst]}, {J~{iSub}} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- EndFor
+ ElseIf (N==3)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}),
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
+ ElseIf (N==5)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}),
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ Else // IF Flag_UpdateSeparated==1
+ // NL part ...
+ // H FIELD >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ Galerkin { [ (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}), {d a} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ ElseIf(N==3)
+ Galerkin { [ (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}), {d a} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ ElseIf(N==5)
+ Galerkin { [ (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}), {d a} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ // DHDB TENSOR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
- // NL part ...
- Galerkin { [ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}), {d a} ];
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
-
- If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
- Galerkin { JacNL[ (dhdb_Vinch_K[ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J_1},{J_1}[1],
- {J_2},{J_2}[1],
- {J_3},{J_3}[1],
- {J_4},{J_4}[1],
- {J_5},{J_5}[1]]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Else
- Galerkin { [ (1/$relax) * (SetVariable[(dhdb_Vinch_K[ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J_1},{J_1}[1],
- {J_2},{J_2}[1],
- {J_3},{J_3}[1],
- {J_4},{J_4}[1],
- {J_5},{J_5}[1]]{List[param_EnergHyst]}) ,
- QuadraturePointIndex[]]{$dhdb_Vinch_temp}) * Dof{d a} , {d a} ];
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
- Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_Vinch_temp}) * {d a} , {d a} ];
- In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ ElseIf (N==3)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
+ ElseIf (N==5)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ EndIf
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
EndIf
-
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
EndIf
//*******************************************************
//*******************************************************
+ Galerkin { DtDof [ ((Thickness^2)/12)*sigma[] * Dof{d a} , {d a} ];
+ In DomainLam; Jacobian Vol; Integration I1; } //kj: Momentally removed because no equivalent in h-form for now
+
Galerkin { [ nxh[], {a} ];
In Boundary ; Jacobian Vol ; Integration I1 ; }
@@ -611,6 +828,59 @@ Formulation {
EndIf
}
}
+
+If (Flag_UpdateSeparated)
+ { Name UpdateH ; Type FemEquation ;
+ Quantity {
+ //{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
+ //{ Name ap ; Type Local ; NameOfSpace Hcurl_a_2D ; } // "ap = aprime usefull to avoid auto-symmetrization of JacNL tensor with getdp"
+ { Name h ; Type Local ; NameOfSpace H_hysteresis ; }
+ /*
+ For iSub In {1:N}
+ { Name X~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
+ EndFor
+ */
+ }
+ Equation {
+ Galerkin { [ Dof{h} ,{h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ /*
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ /*/
+ Galerkin { [ - (GetVariable[ ElementNum[], QuadraturePointIndex[]]{$h_EB_temp_separated}), {h} ];
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ //*/
+ }
+ }
+
+
+ For iSub In {1 : N}
+ { Name UpdateCell~{iSub} ; Type FemEquation ;
+ Quantity {
+ //{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
+ //{ Name ap ; Type Local ; NameOfSpace Hcurl_a_2D ; } // "ap = aprime usefull to avoid auto-symmetrization of JacNL tensor with getdp"
+ //{ Name h ; Type Local ; NameOfSpace H_hysteresis ; }
+ { Name X~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
+ }
+ Equation {
+ Galerkin { [ Dof{X~{iSub}} ,{X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ Galerkin { [ - Cell_EB[ iSub ,
+ (GetVariable[ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X~{iSub}}[1] ]{List[param_EnergHyst]}, {X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ }
+ }
+ EndFor
+
+EndIf
+
+
}
//-----------------------------------------------------------------------------------------------
@@ -631,6 +901,14 @@ Resolution {
If(Flag_AnalysisType==AN_STATIC || Flag_AnalysisType==AN_TIME )
{ Name A ; NameOfFormulation MagStaDyn_a_2D ; }
EndIf
+
+ If (Flag_UpdateSeparated)
+ { Name UpH ; NameOfFormulation UpdateH ; }
+ For iSub In {1 : N}
+ { Name UpCell~{iSub} ; NameOfFormulation UpdateCell~{iSub} ; }
+ EndFor
+ EndIf
+
}
Operation {
@@ -671,15 +949,24 @@ Resolution {
If(Flag_AnalysisType==AN_TIME)
+ SetExtrapolationOrder[Flag_ExtrapolationOrder];
+
// set some runtime variables
Evaluate[$syscount = 0];
Evaluate[$relaxcount=0];
Evaluate[$relaxcounttot=0];
Evaluate[$dnccount=0];
Evaluate[$itertot=0];
+ Evaluate[$itermax=0];
//Save TimeStep 0
SaveSolution[A];
+ If (Flag_UpdateSeparated)
+ InitSolution[UpH]; SaveSolution[UpH];
+ For iSub In {1 : N}
+ InitSolution[UpCell~{iSub}]; SaveSolution[UpCell~{iSub}];
+ EndFor
+ EndIf
If (Flag_LiveLocalPostOp)
PostOperation[Get_LocalFields] ;
EndIf
@@ -688,12 +975,23 @@ Resolution {
EndIf
TimeLoopTheta[time0, timemax, delta_time, 1.]{ // Euler implicit (1) -- Crank-Nicolson (0.5)
+
+ If (Flag_UpdateSeparated)
+ // IMPORTANT to init the current time step in order to get correctly
+ // the value from the previous time step {}[1] in the main System A.
+ InitSolution[UpH]; //SaveSolution[UpH];
+ For iSub In {1 : N}
+ InitSolution[UpCell~{iSub}]; //SaveSolution[UpCell~{iSub}];
+ EndFor
+ EndIf
+
If(!Flag_NL) // Linear case ...
Generate[A]; Solve[A];
EndIf
If(Flag_NL) // NonLinear Case ...
Evaluate[$relax=relax_max];
+
If(Flag_NLRes==NLRES_ITERATIVELOOP)
// I T E R A T I V E . L O O P ...................................
IterativeLoop[ Nb_max_iter, stop_criterion, RF_tuned[]] {
@@ -718,26 +1016,42 @@ Resolution {
If(Flag_NLRes==NLRES_ITERATIVELOOPN)
// I T E R A T I V E . L O O P . N ..............................
- // INIT for h_only case::::::::::::
- /*
+ //:::::::::::::::::::::::::::::::::::::::::::
+ If(Flag_ItLoopNDoFirstIter)
GenerateJac[A] ;
If (!Flag_AdaptRelax)
SolveJac[A] ;
Else
SolveJac_AdaptRelax[A, List[RelaxFac_Lin],TestAllFactors ] ;
EndIf
- //*/
- //:::::::::::::::::::::::::::::::::
-
+ EndIf
+ //:::::::::::::::::::::::::::::::::::::::::::
+ //*****Choose between one of the following possibilities:*****
+ If(Flag_ItLoopNRes==NLITLOOPN_SOLUTION)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ System { { A , Reltol_Mag, Abstol_Mag, Solution MeanL2Norm }} ]{ // 0 : Solution (dx=x-xp; x=x)
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_RESIDUAL)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ System { { A , Reltol_Mag, Abstol_Mag, Residual MeanL2Norm }} ]{ //1 : Residual (dx=res=b-Ax; x=b)
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_RECALCRESIDUAL)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ System { { A , Reltol_Mag, Abstol_Mag, RecalcResidual MeanL2Norm }} ]{ //2 : RecalcResidual (dx=res=b-Ax; x=b)
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_POSTOPX)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ PostOperation { { az_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ // 31 : PostOp unknown field az
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_POSTOPB)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ PostOperation { { b_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ // 32 : PostOp b field
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_POSTOPH)
IterativeLoopN[ Nb_max_iter, RF_tuned[],
- //*****Choose between one of the 3 following possibilities:*****
- //System { { A , Reltol_Mag, Abstol_Mag, Solution MeanL2Norm }} ]{ //1a) //dx=x-xp; x=x
- //System { { A , Reltol_Mag, Abstol_Mag, RecalcResidual MeanL2Norm }} ]{ //1b) //dx=res=b-Ax; x=b
- //System { { A , Reltol_Mag, Abstol_Mag, Residual MeanL2Norm }} ]{ //1c) //dx=res=b-Ax; x=b #(default for square) #CHECK here
- //PostOperation { { az_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //2)
- //PostOperation { { b_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //3)
- PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //4) //(default for t32) #CHECK here //Need the above "INIT" for square with EnergHyst or JA because h_only = 0 at iter 1
- //**************************************************************
+ PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ // 33 : PostOp h field
+ EndIf
+ //**************************************************************
Evaluate[$res = $ResidualN, $resL = $Residual,$resN = $ResidualN,$iter = $Iteration-1];
Test[ $iter>0]{
@@ -768,6 +1082,10 @@ Resolution {
//...............................................................
EndIf
+
+ Test[$iter>$itermax]{
+ Evaluate[$itermax = $iter];
+ }
Evaluate[$itertot = $itertot+$iter];
Test[ (Flag_NLRes==NLRES_ITERATIVELOOPPRO && $iter == Nb_max_iter) ||
(Flag_NLRes==NLRES_ITERATIVELOOP && $res > Abstol_Mag ) ||
@@ -778,9 +1096,15 @@ Resolution {
Evaluate[$dnccount=$dnccount+1];
}
- EndIf
+ EndIf // END NONLINEAR CASE
SaveSolution[A];
+ If (Flag_UpdateSeparated)
+ Generate[UpH]; Solve[UpH]; SaveSolution[UpH];
+ For iSub In {1 : N}
+ Generate[UpCell~{iSub}]; Solve[UpCell~{iSub}]; SaveSolution[UpCell~{iSub}];
+ EndFor
+ EndIf
If (Flag_LiveLocalPostOp)
PostOperation[Get_LocalFields] ;
EndIf
@@ -789,7 +1113,7 @@ Resolution {
PostOperation[Get_GlobalQuantities];
EndIf
//}
- }
+ } // End TIMELOOPTHETA
Print["--------------------------------------------------------------"];
Print[{$syscount}, Format "Total number of linear systems solved: %g"];
Print[{$relaxcounttot}, Format "Total number of relaxation factor tested: %g"];
@@ -797,6 +1121,10 @@ Resolution {
Print[{$dnccount}, Format "Total number of non converged TimeStep: %g"];
Print["--------------------------------------------------------------"];
EndIf // End Flag_AnalysisType==AN_TIME (Time domain)
+
+ Print["syscount relaxcounttot meanrelaxcount iterFEtot meaniterFE maxiterFE dnccount CPUtime"];
+ Print[{$syscount, $relaxcounttot, $relaxcounttot/$syscount, $itertot,$itertot/$TimeStep, $itermax, $dnccount,GetCpuTime[]}, Format "%g %g %g %g %g %g %g %g"];
+
}
}
}
@@ -1028,9 +1356,12 @@ PostOperation Get_GlobalQuantities UsingPost MagStaDyn_a_2D {
PostOperation Get_AllTS UsingPost MagStaDyn_a_2D {
+ ///*
+ Print[ ir, OnElementsOf Inds, File StrCat[Dir,"ir_all",ExtGmsh] ] ;
Print[ az, OnElementsOf Domain, File StrCat[Dir,"a_all",ExtGmsh] ];
Print[ b, OnElementsOf Domain, File StrCat[Dir,"b_all",ExtGmsh] ];
Print[ h, OnElementsOf Domain, File StrCat[Dir,"h_all",ExtGmsh] ];
+ //*/
For k In {1:num_postop_points}
Print[ hb, OnPoint {xpos~{k},ypos~{k},0}, Format TimeTable,
@@ -1039,6 +1370,7 @@ PostOperation Get_AllTS UsingPost MagStaDyn_a_2D {
Print[ I, OnRegion IndA_1, Format Table, File StrCat[Dir,"I1_all",ExtGnuplot] ];
Print[ I, OnRegion IndA_2, Format Table, File StrCat[Dir,"I2_all",ExtGnuplot] ];
+
}
PostOperation {
diff --git a/Team32/magstadyn_avs_3d.pro b/Team32/magstadyn_avs_3d.pro
index fb5f3b384d2b5740a976040e034ba9503d0e9a93..1f0f2e1f0b2eb03eb17ccfc0bcfcec1acdc19d54 100644
--- a/Team32/magstadyn_avs_3d.pro
+++ b/Team32/magstadyn_avs_3d.pro
@@ -403,26 +403,43 @@ Formulation {
// EnergHyst++b
If(Flag_NL_law==NL_ENERGHYST)
For iSub In {1:N}
- { Name J~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
+ { Name X~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
EndFor
// EnergHyst++e
EndIf
}
- Equation {
If(Flag_NL_law==NL_LIN || Flag_NL_law==NL_ANA || Flag_NL_law==NL_ANA_JA )
- Galerkin { [ nu[{d a}] * Dof{d a} , {d a} ] ;
+ // nu-Version) with nu[{d a}] * Dof{d a} and dhdb_NL[...]
+ //-------------------------------------------------------
+ /*
+ Galerkin { [ nu[{d a}] * Dof{d a} , {d a} ] ;
In Domain ; Jacobian Vol ; Integration II ; }
If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
Galerkin { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d ap} ] ;
In DomainNL ; Jacobian Vol ; Integration II ; }
Else
- Galerkin { [ (1/$relax) * (SetVariable[(dhdb_NL[{d a}]),
- QuadraturePointIndex[]]{$dhdb_any_temp}) * Dof{d a} , {d a} ] ;
+ Galerkin { [ ((1/$relax) * (SetVariable[(dhdb_NL[{d a}]),
+ QuadraturePointIndex[]]{$dhdb_any_temp})) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_any_temp}) * {d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II ; }
+ /*/
+ // h-Version) with {h}=nu[{d a}]*{d a} and dhdb[...]
+ //--------------------------------------------------
+ Galerkin { [ nu[{d a}] * {d a} , {d a} ] ;
+ In Domain ; Jacobian Vol ; Integration II ; }
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL [ dhdb[{d a}] * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II ; }
+ Else
+ Galerkin { [ ((1/$relax) * (SetVariable[(dhdb[{d a}]),
+ QuadraturePointIndex[]]{$dhdb_any_temp})) * Dof{d a} , {d ap} ] ;
In DomainNL ; Jacobian Vol ; Integration II ; }
- Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_any_temp}) * {d a} , {d a} ] ;
+ Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_any_temp}) * {d a} , {d ap} ] ;
In DomainNL ; Jacobian Vol ; Integration II ; }
- EndIf
+ //*/
+ EndIf
EndIf
//***************************************************
@@ -437,84 +454,318 @@ Formulation {
// BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
Galerkin { [ Dof{h} , {h} ] ;
In DomainNL ; Jacobian Vol ; Integration II1p ; }
- Galerkin { [ -h_Jiles[{h}[1] ,
- {d a}[1],
- {d a} ]{List[hyst_FeSi]} , {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration II1p ; }
-
+
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[( h_Jiles[ {h}[1] ,
+ {d a}[1],
+ {d a} ]{List[hyst_FeSi]}),
+ QuadraturePointIndex[]]{$h_Jiles_temp}) , {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else
+ Galerkin { [ - (SetVariable[ $relax * ( SetVariable[ h_Jiles[ {h}[1] ,
+ {d a}[1],
+ {d a} ]{List[hyst_FeSi]}
+ ,QuadraturePointIndex[]]{$h_Jiles_unrelaxed}
+ )
+ +(1-$relax)* {h},
+ QuadraturePointIndex[]
+ ]{$h_Jiles_temp}) , {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
// NL part ...
- Galerkin { [ Dof{h}, {d a} ] ; // register of Dof => idem
- In DomainNL ; Jacobian Vol ; Integration II ; }
+ //Galerkin { [ Dof{h}, {d a} ] ; // register of Dof => idem
+ Galerkin { [ (GetVariable[ QuadraturePointIndex[]]{$h_Jiles_temp}), {d a} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
- Galerkin { JacNL[ dhdb_Jiles[{h} ,
- {d a} ,
- {h}-{h}[1] ]{List[hyst_FeSi]} * Dof{d a} , {d ap} ] ;
- In DomainNL ; Jacobian Vol ; Integration II ; }
+ Galerkin { JacNL[ dhdb_Jiles[ {h} ,
+ {d a} ,
+ {h}-{h}[1] ]{List[hyst_FeSi]} * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; } // kj modif II --> II1p here
Else
Galerkin { [ (1/$relax) * (SetVariable[(dhdb_Jiles[{h} ,
{d a} ,
{h}-{h}[1] ]{List[hyst_FeSi]}),
- QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * Dof{d a} , {d a} ];
- In DomainNL ; Jacobian Vol ; Integration II ; }
- Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * {d a} , {d a} ];
- In DomainNL ; Jacobian Vol ; Integration II ; }
+ QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; } // kj modif II --> II1p here
+ Galerkin { [ -(1/$relax) * (GetVariable[ QuadraturePointIndex[]]{$dhdb_Jiles_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; } // kj modif II --> II1p here
EndIf
EndIf
//*******************************************************
//--------Hysteresis law with EnergHyst model----------
//*******************************************************
- If(Flag_NL_law==NL_ENERGHYST ) //EnergHyst model law
+ If(Flag_NL_law==NL_ENERGHYST) //EnergHyst model law
Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
In DomainL ; Jacobian Vol ; Integration II ; }
- // update h
- // saving h_Vinch_K in local quantity h
- // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
- Galerkin { [ Dof{h} , {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration II1p ; }
- Galerkin { [ - (SetVariable[(h_Vinch_K[ {h}[1] ,
- {d a}, {d a}[1],
- {J_1}, {J_1}[1],
- {J_2}, {J_2}[1],
- {J_3}, {J_3}[1] ]{List[param_EnergHyst]}),
- QuadraturePointIndex[]]{$h_Vinch_temp}), {h} ] ;
- In DomainNL ; Jacobian Vol ; Integration II1p ; }
-
- // update Jk
- // saving Update_Cell_K in local quantity Jk
- // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
- For iSub In {1:N}
- Galerkin { [ Dof{J~{iSub}} ,{J~{iSub}} ] ;
- In DomainNL ; Jacobian Vol ; Integration II1p ; } // Integration with II1p or II ???
- Galerkin { [ - Update_Cell_K[ iSub ,
- (GetVariable[QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J~{iSub}} ,
- {J~{iSub}}[1] ]{List[param_EnergHyst]}, {J~{iSub}} ] ;
- In DomainNL ; Jacobian Vol ; Integration II1p ; } // Integration with II1p or II ???
- EndFor
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ If (!Flag_UpdateSeparated)
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // U P D A T E - H
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // saving h_EB in local quantity h
+ // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
+ Galerkin { [ Dof{h} , {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1]
+ ]{List[param_EnergHyst]}),
+ QuadraturePointIndex[]]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - (SetVariable[ $relax* ( SetVariable[ h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1]
+ ]{List[param_EnergHyst]}
+ , QuadraturePointIndex[]
+ ]{$h_EB_temp_unrelaxed}
+ )
+ +(1-$relax)*{h},
+ QuadraturePointIndex[]
+ ]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+ ElseIf (N==3)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}),
+ QuadraturePointIndex[]]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - (SetVariable[ $relax* ( SetVariable[ h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}
+ , QuadraturePointIndex[]
+ ]{$h_EB_temp_unrelaxed}
+ )
+ +(1-$relax)*{h},
+ QuadraturePointIndex[]
+ ]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+ ElseIf (N==5)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}),
+ QuadraturePointIndex[]]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - (SetVariable[ $relax* ( SetVariable[ h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1], {X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}
+ , QuadraturePointIndex[]
+ ]{$h_EB_temp_unrelaxed}
+ )
+ +(1-$relax)*{h},
+ QuadraturePointIndex[]
+ ]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // U P D A T E - C E L L S
+ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ // saving Cell_EB in local quantity Xk
+ // BF is constant per element (Vector + BF_Volume) => 1 integration point is enough
+ For iSub In {1:N}
+ Galerkin { [ Dof{X~{iSub}} ,{X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { [ - Cell_EB[ iSub ,
+ (GetVariable[QuadraturePointIndex[]]{$h_EB_temp}),
+ {X~{iSub}}[1] ]{List[param_EnergHyst]}, {X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ - ( $relax* Cell_EB[ iSub ,
+ (GetVariable[QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X~{iSub}}[1] ]{List[param_EnergHyst]}
+ +(1-$relax)*{X~{iSub}} ), {X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
- // NL part ...
- Galerkin { [ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}), {d a} ];
- In DomainNL ; Jacobian Vol ; Integration II1p ; } // Integration with II1p or II ???
-
- If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
- Galerkin { JacNL[ (dhdb_Vinch_K[ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J_1},{J_1}[1],
- {J_2},{J_2}[1],
- {J_3},{J_3}[1]]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
- In DomainNL ; Jacobian Vol ; Integration II1p ; } // Integration with II1p or II ???
- Else
- Galerkin { [ (1/$relax) * (SetVariable[(dhdb_Vinch_K[ (GetVariable[ QuadraturePointIndex[]]{$h_Vinch_temp}),
- {J_1},{J_1}[1],
- {J_2},{J_2}[1],
- {J_3},{J_3}[1]]{List[param_EnergHyst]}) ,
- QuadraturePointIndex[]]{$dhdb_Vinch_temp}) * Dof{d a} , {d a} ];
- In DomainNL ; Jacobian Vol ; Integration II1p ; } // Integration with II1p or II ???
- Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_Vinch_temp}) * {d a} , {d a} ];
- In DomainNL ; Jacobian Vol ; Integration II1p ; } // Integration with II1p or II ???
+ EndFor
+
+ // NL part ...
+
+ // H FIELD >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ Galerkin { [ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}), {d a} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ // DHDB TENSOR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}),
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+
+ ElseIf (N==3)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}),
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+ ElseIf (N==5)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp}),
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ QuadraturePointIndex[]]{$h_EB_temp_unrelaxed}), // _unrelaxed or not ?
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ Else // IF Flag_UpdateSeparated==1
+ // NL part ...
+ // H FIELD >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ Galerkin { [ (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}), {d a} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ ElseIf(N==3)
+ Galerkin { [ (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}), {d a} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ ElseIf(N==5)
+ Galerkin { [ (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}), {d a} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+ // DHDB TENSOR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ // >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+ //.......................................................................................
+ // TODO: avoid the following duplications if one could group {X_1}[1],{X_2}[1],{X_3}[1],...
+ If (N==1)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+
+ ElseIf (N==3)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+ ElseIf (N==5)
+ If(Flag_NLRes!=NLRES_ITERATIVELOOPPRO)
+ Galerkin { JacNL[ (dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) * Dof{d a} , {d ap} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Else // IF Flag_NLRes==NLRES_ITERATIVELOOPPRO
+ Galerkin { [ (1/$relax) * (SetVariable[(dhdb_EB[ (GetVariable[ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X_1}[1],{X_2}[1],{X_3}[1],{X_4}[1],{X_5}[1]
+ ]{List[param_EnergHyst]}) ,
+ QuadraturePointIndex[]]{$dhdb_EB_temp}) * Dof{d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Galerkin { [ -(1/$relax) * (GetVariable[QuadraturePointIndex[]]{$dhdb_EB_temp}) * {d a} , {d ap} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ EndIf
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The implementation with the asked number of cells (N) is not written yet in the formulation in 'magstadyn_a.pro'"];
+ EndIf
+ //.......................................................................................
+ // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
EndIf
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////
EndIf
//*******************************************************
//*******************************************************
@@ -573,6 +824,52 @@ Formulation {
}
}
+If (Flag_UpdateSeparated)
+ { Name UpdateH ; Type FemEquation ;
+ Quantity {
+ //{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
+ //{ Name ap ; Type Local ; NameOfSpace Hcurl_a_2D ; } // "ap = aprime usefull to avoid auto-symmetrization of JacNL tensor with getdp"
+ { Name h ; Type Local ; NameOfSpace H_hysteresis ; }
+ //{ Name X~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
+ }
+ Equation {
+ Galerkin { [ Dof{h} ,{h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ /*
+ Galerkin { [ - (SetVariable[(h_EB[ {h} , // use {h} or {h}[1] here ? (need to init bc by recomputing b from {h} in h_EB)
+ {d a} ,
+ {X_1}[1],{X_2}[1],{X_3}[1]
+ ]{List[param_EnergHyst]}),
+ ElementNum[],QuadraturePointIndex[]]{$h_EB_temp}), {h} ] ;
+ In DomainNL ; Jacobian Vol ; Integration I1p ; }
+ */
+ Galerkin { [ (GetVariable[ ElementNum[], QuadraturePointIndex[]]{$h_EB_temp_separated}), {h} ];
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ }
+ }
+
+
+ For iSub In {1 : N}
+ { Name UpdateCell~{iSub} ; Type FemEquation ;
+ Quantity {
+ //{ Name a ; Type Local ; NameOfSpace Hcurl_a_2D ; }
+ //{ Name ap ; Type Local ; NameOfSpace Hcurl_a_2D ; } // "ap = aprime usefull to avoid auto-symmetrization of JacNL tensor with getdp"
+ //{ Name h ; Type Local ; NameOfSpace H_hysteresis ; }
+ { Name X~{iSub} ; Type Local ; NameOfSpace Vect~{iSub};}
+ }
+ Equation {
+ Galerkin { [ Dof{X~{iSub}} ,{X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ Galerkin { [ - Cell_EB[ iSub ,
+ (GetVariable[ElementNum[],QuadraturePointIndex[]]{$h_EB_temp_separated}),
+ {X~{iSub}}[1] ]{List[param_EnergHyst]}, {X~{iSub}} ] ;
+ In DomainNL ; Jacobian Vol ; Integration II1p ; }
+ }
+ }
+ EndFor
+
+EndIf
+
}
//-----------------------------------------------------------------------------------------------
@@ -594,12 +891,21 @@ Resolution {
If(Flag_AnalysisType==AN_STATIC || Flag_AnalysisType==AN_TIME )
{ Name A ; NameOfFormulation MagStaDyn_a_3D ; }
EndIf
+
+ If (Flag_UpdateSeparated)
+ { Name UpH ; NameOfFormulation UpdateH ; }
+ For iSub In {1 : N}
+ { Name UpCell~{iSub} ; NameOfFormulation UpdateCell~{iSub} ; }
+ EndFor
+ EndIf
}
Operation {
- CreateDir["res3d/"] ;
+ CreateDir[Dir];
DeleteFile[StrCat[Dir,"I1",ExtGnuplot]];
DeleteFile[StrCat[Dir,"I2",ExtGnuplot]];
- For k In{1:2}
+ DeleteFile[StrCat[Dir,"U1",ExtGnuplot]];
+ DeleteFile[StrCat[Dir,"U2",ExtGnuplot]];
+ For k In{1:num_postop_points}
DeleteFile[StrCat[Dir,Sprintf("hbp_%g",k),ExtGnuplot]];
DeleteFile[StrCat[Dir,Sprintf("bx_%g",k),ExtGnuplot]];
DeleteFile[StrCat[Dir,Sprintf("by_%g",k),ExtGnuplot]];
@@ -627,19 +933,41 @@ Resolution {
If(Flag_AnalysisType==AN_TIME)
+ SetExtrapolationOrder[Flag_ExtrapolationOrder];
+
// set some runtime variables
Evaluate[$syscount = 0];
Evaluate[$relaxcount=0];
Evaluate[$relaxcounttot=0];
Evaluate[$dnccount=0];
Evaluate[$itertot=0];
+ Evaluate[$itermax=0];
//Save TimeStep 0
SaveSolution[A];
+ If (Flag_UpdateSeparated)
+ InitSolution[UpH]; SaveSolution[UpH];
+ For iSub In {1 : N}
+ InitSolution[UpCell~{iSub}]; SaveSolution[UpCell~{iSub}];
+ EndFor
+ EndIf
+ If (Flag_LiveLocalPostOp)
PostOperation[Get_LocalFields] ;
- PostOperation[Get_GlobalQuantities];
+ EndIf
+ If (Flag_LiveGlobalPostOp)
+ PostOperation[Get_GlobalQuantities];
+ EndIf
TimeLoopTheta[time0, timemax, delta_time, 1.]{ // Euler implicit (1) -- Crank-Nicolson (0.5)
+ If (Flag_UpdateSeparated)
+ // IMPORTANT to init the current time step in order to get correctly
+ // the value from the previous time step {}[1] in the main System A.
+ InitSolution[UpH]; //SaveSolution[UpH];
+ For iSub In {1 : N}
+ InitSolution[UpCell~{iSub}]; //SaveSolution[UpCell~{iSub}];
+ EndFor
+ EndIf
+
If(!Flag_NL) // Linear case ...
Generate[A]; Solve[A];
EndIf
@@ -670,25 +998,42 @@ Resolution {
If(Flag_NLRes==NLRES_ITERATIVELOOPN)
// I T E R A T I V E . L O O P . N ..............................
- //:::::::::::::::::::::::::::::::::
- // INIT for h_only case::::::::::::
- /*
+ //:::::::::::::::::::::::::::::::::::::::::::
+ If(Flag_ItLoopNDoFirstIter)
GenerateJac[A] ;
If (!Flag_AdaptRelax)
SolveJac[A] ;
Else
SolveJac_AdaptRelax[A, List[RelaxFac_Lin],TestAllFactors ] ;
EndIf
- //*/
- //:::::::::::::::::::::::::::::::::
-
+ EndIf
+ //:::::::::::::::::::::::::::::::::::::::::::
+ //*****Choose between one of the following possibilities:*****
+ If(Flag_ItLoopNRes==NLITLOOPN_SOLUTION)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ System { { A , Reltol_Mag, Abstol_Mag, Solution MeanL2Norm }} ]{ // 0 : Solution (dx=x-xp; x=x)
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_RESIDUAL)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ System { { A , Reltol_Mag, Abstol_Mag, Residual MeanL2Norm }} ]{ //1 : Residual (dx=res=b-Ax; x=b)
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_RECALCRESIDUAL)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ System { { A , Reltol_Mag, Abstol_Mag, RecalcResidual MeanL2Norm }} ]{ //2 : RecalcResidual (dx=res=b-Ax; x=b)
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_POSTOPX)
IterativeLoopN[ Nb_max_iter, RF_tuned[],
- //*****Choose between one of the 3 following possibilities:*****
- //System { { A , Reltol_Mag, Abstol_Mag, Residual MeanL2Norm }} ]{ //1)
- //PostOperation { { a_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //2)
- //PostOperation { { b_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //3)
- PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ //4) //Need the above "INIT" for square with EnergHyst or JA because h_only = 0 at iter 1
- //**************************************************************
+ PostOperation { { az_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ // 31 : PostOp unknown field az
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_POSTOPB)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ PostOperation { { b_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ // 32 : PostOp b field
+ EndIf
+ If(Flag_ItLoopNRes==NLITLOOPN_POSTOPH)
+ IterativeLoopN[ Nb_max_iter, RF_tuned[],
+ PostOperation { { h_only , Reltol_Mag, Abstol_Mag, MeanL2Norm }} ]{ // 33 : PostOp h field
+ EndIf
+ //**************************************************************
Evaluate[$res = $ResidualN, $resL = $Residual,$resN = $ResidualN,$iter = $Iteration-1];
Test[ $iter>0]{
@@ -720,6 +1065,9 @@ Resolution {
//...............................................................
EndIf
+ Test[$iter>$itermax]{
+ Evaluate[$itermax = $iter];
+ }
Evaluate[$itertot = $itertot+$iter];
Test[ (Flag_NLRes==NLRES_ITERATIVELOOPPRO && $iter == Nb_max_iter) ||
(Flag_NLRes==NLRES_ITERATIVELOOP && $res > Abstol_Mag ) ||
@@ -730,9 +1078,15 @@ Resolution {
Evaluate[$dnccount=$dnccount+1];
}
- EndIf
+ EndIf // END NONLINEAR CASE
SaveSolution[A];
+ If (Flag_UpdateSeparated)
+ Generate[UpH]; Solve[UpH]; SaveSolution[UpH];
+ For iSub In {1 : N}
+ Generate[UpCell~{iSub}]; Solve[UpCell~{iSub}]; SaveSolution[UpCell~{iSub}];
+ EndFor
+ EndIf
If (Flag_LiveLocalPostOp)
PostOperation[Get_LocalFields] ;
EndIf
@@ -741,7 +1095,7 @@ Resolution {
PostOperation[Get_GlobalQuantities];
EndIf
//}
- }
+ } // End TIMELOOPTHETA
Print["--------------------------------------------------------------"];
Print[{$syscount}, Format "Total number of linear systems solved: %g"];
Print[{$relaxcounttot}, Format "Total number of relaxation factor tested: %g"];
@@ -749,6 +1103,10 @@ Resolution {
Print[{$dnccount}, Format "Total number of non converged TimeStep: %g"];
Print["--------------------------------------------------------------"];
EndIf // End Flag_AnalysisType==AN_TIME (Time domain)
+
+ Print["syscount relaxcounttot meanrelaxcount iterFEtot meaniterFE maxiterFE dnccount CPUtime"];
+ Print[{$syscount, $relaxcounttot, $relaxcounttot/$syscount, $itertot,$itertot/$TimeStep, $itermax, $dnccount,GetCpuTime[]}, Format "%g %g %g %g %g %g %g %g"];
+
}
}
}
@@ -957,6 +1315,66 @@ PostProcessing {
}
+PostOperation Get_AllTS_3D UsingPost MagStaDyn_a_3D {
+ /*
+ Print[ ir, OnElementsOf Inds, File StrCat[Dir,"ir_all",ExtGmsh] ] ;
+ Print[ az, OnElementsOf Domain, File StrCat[Dir,"a_all",ExtGmsh] ];
+ Print[ b, OnElementsOf Domain, File StrCat[Dir,"b_all",ExtGmsh] ];
+ Print[ h, OnElementsOf Domain, File StrCat[Dir,"h_all",ExtGmsh] ];
+ */
+
+ Print[ hs, OnElementsOf DomainB, File StrCat[Dir,"hs_all",ExtGmsh] ] ;
+ Print[ js, OnElementsOf DomainB, File StrCat[Dir,"js_all",ExtGmsh] ] ;
+ Print[ b, OnElementsOf Domain, File StrCat[Dir,"b_all",ExtGmsh] ];
+
+ For k In {1:num_postop_points}
+ Print[ hb, OnPoint {xpos~{k},ypos~{k},0}, Format TimeTable,
+ File StrCat[Dir,Sprintf("hbp_%g_all",k),ExtGnuplot] ];
+ EndFor
+
+ Print[ I, OnRegion Ind_1, Format Table, File StrCat[Dir,"I1_all",ExtGnuplot] ];
+ Print[ I, OnRegion Ind_2, Format Table, File StrCat[Dir,"I2_all",ExtGnuplot] ];
+
+}
+
+ PostOperation Get_GlobalQuantities_ALLTS_3D UsingPost MagStaDyn_a_3D {
+
+ If (!Flag_CircuitCoupling)
+ Print[ I, OnRegion Ind_1, Format Table,
+ File > StrCat[Dir,"I1",ExtGnuplot]];
+ Print[ I, OnRegion Ind_2, Format Table,
+ File > StrCat[Dir,"I2",ExtGnuplot]];
+
+ Print[ U, OnRegion Ind_1, Format Table,
+ File > StrCat[Dir,"U1_all",ExtGnuplot]];
+ Print[ U, OnRegion Ind_2, Format Table,
+ File > StrCat[Dir,"U2_all",ExtGnuplot]];
+ EndIf
+
+ If (Flag_CircuitCoupling)
+ Print[ I, OnRegion Ein1, Format Table,
+ File > StrCat[Dir,"I1_all",ExtGnuplot] ];
+
+ Print[ U, OnRegion Ein1, Format Table,
+ File > StrCat[Dir,"U1_all",ExtGnuplot] ];
+
+ If (Flag_TestCase==3)
+ Print[ I, OnRegion Ein2, Format Table, // Iminus because of the sign convention to be in agreement with measurements team 32
+ File > StrCat[Dir,"I2_all",ExtGnuplot]];
+
+ Print[ U, OnRegion Ein2, Format Table,
+ File > StrCat[Dir,"U2_all",ExtGnuplot] ];
+ EndIf
+
+ If (Flag_TestCase==4)
+ Print[ I, OnRegion Ind_2, Format Table, // Iminus because of the sign convention to be in agreement with measurements team 32
+ File > StrCat[Dir,"I2_all",ExtGnuplot] ];
+
+ EndIf
+ EndIf
+}
+
+
PostOperation Get_LocalFields_ALLTS UsingPost MagStaDyn_a_3D {
If(Flag_AnalysisType!=AN_TIME)
diff --git a/Team32/param_EnergHyst.dat b/Team32/param_EnergHyst.dat
index cf500f08d8922cfd709a3c95588c1cbd9ba44778..078d799d09fbbdddb3dd805e6a12b926c79cc3f2 100644
--- a/Team32/param_EnergHyst.dat
+++ b/Team32/param_EnergHyst.dat
@@ -1,10 +1,10 @@
///*//default for t32 #CHECK here
dim00 = (Flag_3Dmodel==1)?3:2;
N00 = 3;
-FLAG_TANORLANG00 = 1;
-FLAG_VARORDIFF00 = 2;
+FLAG_ANHYLAW00 = 2;
+FLAG_APPROACH00 = 2;
FLAG_INVMETHOD00 = 1;
-FLAG_ROOTFINDING1D00 = 0;
+FLAG_JACEVAL00 = 1;
FLAG_MINMETHOD00 = 1;
FLAG_WARNING00 = 0 ;
TOLERANCE_JS00 = 1.e-3;
@@ -23,12 +23,12 @@ FLAG_HOMO00 = 1;
/*
dim00 = (Flag_3Dmodel==1)?3:2;
N00 = 3;
-FLAG_TANORLANG00 = 1;
-FLAG_VARORDIFF00 = 2;
+FLAG_ANHYLAW00 = 1;
+FLAG_APPROACH00 = 2;
FLAG_INVMETHOD00 = 1;
-FLAG_ROOTFINDING1D00 = 0;
+FLAG_JACEVAL00 = 1;
FLAG_MINMETHOD00 = 1;
-FLAG_WARNING00 = 3 ;
+FLAG_WARNING00 = 0 ;
TOLERANCE_JS00 = 1.e-3;
TOLERANCE_000 = 1.e-8;
TOLERANCE_NR00 = 1.e-8;
diff --git a/Team32/param_EnergHyst.pro b/Team32/param_EnergHyst.pro
index 018a91aa8c14c541fec7e2dd30c59ed7f8d0340b..74a3667dbffcd3bdb7035dbb51a8c70e10a75d4e 100644
--- a/Team32/param_EnergHyst.pro
+++ b/Team32/param_EnergHyst.pro
@@ -5,39 +5,27 @@ Function {
dim = dim00;
N = N00; //nbr de cercle
- FLAG_TANORLANG = FLAG_TANORLANG00;
+ FLAG_ANHYLAW = FLAG_ANHYLAW00;
/*
-FLAG_TANORLANG = 1 --> hyperbolic tangent
-FLAG_TANORLANG = 2 --> double langevin function
+FLAG_ANHYLAW = 1 --> hyperbolic tangent
+FLAG_ANHYLAW = 2 --> double langevin function
*/
- FLAG_VARORDIFF = FLAG_VARORDIFF00;
+ FLAG_APPROACH = FLAG_APPROACH00;
/*
-FLAG_VARORDIFF = 1 --> variational approach (Jk)
-FLAG_VARORDIFF = 2 --> simple differential approach (hrk)
-FLAG_VARORDIFF = 3 --> full differential approach (hrk)
+FLAG_APPROACH = 1 --> variational approach (Jk)
+FLAG_APPROACH = 2 --> vector play model approach (hrk)
+FLAG_APPROACH = 3 --> full differential approach (hrk)
*/
FLAG_INVMETHOD = FLAG_INVMETHOD00;
/*
FLAG_INVMETHOD = 1 --> NR_ana (homemade)
FLAG_INVMETHOD = 2 --> NR_num (homemade)
FLAG_INVMETHOD = 3 --> bfgs (homemade)
-FLAG_INVMETHOD = 4 --> hybridsj (gsl)
-FLAG_INVMETHOD = 5 --> hybridj (gsl)
-FLAG_INVMETHOD = 6 --> newton (gsl)
-FLAG_INVMETHOD = 7 --> gnewton (gsl)
*/
- FLAG_ROOTFINDING1D = FLAG_ROOTFINDING1D00;
+ FLAG_JACEVAL = FLAG_JACEVAL00;
/*
-FLAG_ROOTFINDING1D = 0 --> Not done
-FLAG_ROOTFINDING1D = 1 --> (1D minimization) golden section (gsl)
-FLAG_ROOTFINDING1D = 2 --> (1D minimization) brent (gsl)
-FLAG_ROOTFINDING1D = 3 --> (1D minimization) quad golden (gsl)
-FLAG_ROOTFINDING1D = 4 --> (1D root bracketing) bisection (gsl)
-FLAG_ROOTFINDING1D = 5 --> (1D root bracketing) brent (gsl)
-FLAG_ROOTFINDING1D = 6 --> (1D root bracketing) falsepos (gsl)
-FLAG_ROOTFINDING1D = 7 --> (1D root finding with derivatives) newton (gsl)
-FLAG_ROOTFINDING1D = 8 --> (1D root finding with derivatives) secant (gsl)
-FLAG_ROOTFINDING1D = 9 --> (1D root finding with derivatives) steffenson (gsl)
+FLAG_JACEVAL = 1 --> analytical Jacobian
+FLAG_JACEVAL = 2 --> numerical Jacobian
*/
FLAG_MINMETHOD = FLAG_MINMETHOD00;
/*
@@ -47,18 +35,21 @@ FLAG_MINMETHOD = 3 --> conjugate pr (gsl)
FLAG_MINMETHOD = 4 --> bfgs2 (gsl)
FLAG_MINMETHOD = 5 --> bfgs (gsl)
FLAG_MINMETHOD = 6 --> steepest descent (gsl)
+FLAG_MINMETHOD = 11 --> steepest descent+ (homemade)\n"
+FLAG_MINMETHOD = 22 --> conjugate Fletcher-Reeves (homemade)\n"
+FLAG_MINMETHOD = 33 --> conjugate Polak-Ribiere (homemade)\n"
+FLAG_MINMETHOD = 333 --> conjugate Polak-Ribiere+ (homemade)\n"
+FLAG_MINMETHOD = 1999 --> conjugate Dai Yuan 1999 (p.85) (homemade)\n"
+FLAG_MINMETHOD = 2005 --> conjugate Hager Zhang 2005 (p.161) (homemade)\n"
+FLAG_MINMETHOD = 77 --> newton (homemade)\n", ::FLAG_MINMETHOD);
*/
FLAG_WARNING = FLAG_WARNING00 ; // SHOW WARNING
/*
#define FLAG_WARNING_INFO_INV 1
-#define FLAG_WARNING_INFO_ROOTFINDING 2
-#define FLAG_WARNING_INFO_MIN 3
-#define FLAG_WARNING_DISP_INV 11
-#define FLAG_WARNING_DISP_ROOTFINDING 12
-#define FLAG_WARNING_STOP_INV 101
-#define FLAG_WARNING_STOP_ROOTFINDING 102
-#define FLAG_WARNING_ITER 100
+#define FLAG_WARNING_INFO_APPROACH 2
+#define FLAG_WARNING_STOP_INV 10
+#define FLAG_WARNING_DISPABOVEITER 1
*/
TOLERANCE_JS = TOLERANCE_JS00; // SENSITIVE_PARAM (1.e-3) // 1.e-4
TOLERANCE_0 = TOLERANCE_000; // SENSITIVE_PARAM (1.e-7)
@@ -73,14 +64,14 @@ FLAG_MINMETHOD = 6 --> steepest descent (gsl)
// 1.e0 for VinchT & transfo)
FLAG_HOMO = FLAG_HOMO00;
- If(FLAG_TANORLANG==1)
+ If(FLAG_ANHYLAW==1)
// Material M25050A
Ja = 1.22;
ha = 65;
Jb = 0.;
hb = 0.;
EndIf
- If(FLAG_TANORLANG==2)
+ If(FLAG_ANHYLAW==2)
/*
// Material M23535A (default for square) #CHECK here
Ja = 0.595;
@@ -113,86 +104,73 @@ Jb = 0.572282990517
hb = 399.474553851
*/
EndIf
- If(N==3)
-///*
+
+ If(N==1)
+ kappa_1 = 0.;
+ w_1 = 1.;
+
+ ElseIf (N==3)
+/*
// Material M25050A (default for square) #CHECK here
Js_1 = 0.11;
Js_2 = 0.8;
Js_3 = 0.31;
- chi_1 = 0;
- chi_2 = 16;
- chi_3 = 47;
+ kappa_1 = 0;
+ kappa_2 = 16;
+ kappa_3 = 47;
w_1 = 0.090163934426230;
w_2 = 0.655737704918033;
w_3 = 0.254098360655738;
/*/
/* //team32_RD_z_anhyshift (abstract ISEM 2017)
- chi_1 = 0*11.6413;
- chi_2 = 50.1413;
- chi_3 = 107.518;
+ kappa_1 = 0*11.6413;
+ kappa_2 = 50.1413;
+ kappa_3 = 107.518;
w_1 = 0.173524;
w_2 = 0.591535;
w_3 = 0.234941;
- Js_1 = w_1*(Ja+Jb);
- Js_2 = w_2*(Ja+Jb);
- Js_3 = w_3*(Ja+Jb);
//*/
/* //FH(xini=400) // THIS IS PROBLEMATIC (CONVERGENCE PROBLEM ok now) (new try)
- chi_1 = 0.;
- chi_2 = 53.779614004;
- chi_3 = 233.384447699;
+ kappa_1 = 0.;
+ kappa_2 = 53.779614004;
+ kappa_3 = 233.384447699;
w_1 = 0.1048; //+0.1 to converge
w_2 = 0.817943263499; // -0.1 to converge
w_3 = 0.0772371911006;
- Js_1 = w_1*(Ja+Jb);
- Js_2 = w_2*(Ja+Jb);
- Js_3 = w_3*(Ja+Jb);
//*/
/* //FH(xini=400) // NEW 1/9/2017 (chamonix) (default for t32) #CHECK here
- chi_1 = 0.;
- chi_2 = 53.789872;
- chi_3 = 233.910892;
+ kappa_1 = 0.;
+ kappa_2 = 53.789872;
+ kappa_3 = 233.910892;
w_1 = 0.0280702473932;
w_2 = 0.896565361832;
w_3 = 0.0753643907746;
- Js_1 = w_1*(Ja+Jb);
- Js_2 = w_2*(Ja+Jb);
- Js_3 = w_3*(Ja+Jb);
//*/
/*//FH(xini=400) TD // NEW 1/9/2017 (chamonix)
0 0.0565084014449 0.0
1 0.839372763961 65.418004
2 0.104118834594 246.564466
*/
- EndIf
- If(N==1)
- Js_1 = 1;
- Js_2 = 0.;
- Js_3 = 0.;
- chi_1 = 0.;
- chi_2 = 0.;
- chi_3 = 0.;
- w_1 = 1.;
- w_2 = 0.;
- w_3 = 0.;
- EndIf
- If(N==5) //FH(xini=400) // NEW 1/9/2017
- chi_1 = 0.;
- chi_2 = 14.81413;
- chi_3 = 52.317547;
- chi_4 = 102.345013;
- chi_5 = 247.012801;
+///* // NEW 8/8/2018 for team32
+ kappa_1 = 0.;
+ kappa_2 = 54.705531;
+ kappa_3 = 199.248421;
+ w_1 = 0.028745;
+ w_2 = 0.921797;
+ w_3 = 0.049458;
+//*/
+ElseIf (N==5) //FH(xini=400) // NEW 1/9/2017
+ kappa_1 = 0.;
+ kappa_2 = 14.81413;
+ kappa_3 = 52.317547;
+ kappa_4 = 102.345013;
+ kappa_5 = 247.012801;
w_1 = 0.028258057159;
w_2 = 0.087100091228;
w_3 = 0.702372713403;
w_4 = 0.134380357677;
w_5 = 0.0478887805338;
- Js_1 = w_1*(Ja+Jb);
- Js_2 = w_2*(Ja+Jb);
- Js_3 = w_3*(Ja+Jb);
- Js_4 = w_4*(Ja+Jb);
- Js_5 = w_5*(Ja+Jb);
/*//FH(xini=400) TD // NEW 1/9/2017
0 0.0574009326857 0.0
@@ -200,19 +178,17 @@ hb = 399.474553851
2 0.576936798861 62.494953
3 0.189127398303 123.191766
4 0.0646106266745 259.419756*/
- EndIf
+EndIf
- If (N==3)
+ If (N==1)
param_EnergHyst={ dim, N,
Ja, ha, Jb, hb,
- w_1, chi_1,
- w_2, chi_2,
- w_3, chi_3,
+ w_1, kappa_1,
FLAG_INVMETHOD,
- FLAG_ROOTFINDING1D,
- FLAG_VARORDIFF,
+ FLAG_JACEVAL,
+ FLAG_APPROACH,
FLAG_MINMETHOD,
- FLAG_TANORLANG,
+ FLAG_ANHYLAW,
FLAG_WARNING,
TOLERANCE_JS,
TOLERANCE_0,
@@ -225,21 +201,42 @@ hb = 399.474553851
DELTA_0,
FLAG_HOMO
};
- EndIf
-
- If (N==5)
+ ElseIf (N==3)
+ param_EnergHyst={ dim, N,
+ Ja, ha, Jb, hb,
+ w_1, kappa_1,
+ w_2, kappa_2,
+ w_3, kappa_3,
+ FLAG_INVMETHOD,
+ FLAG_JACEVAL,
+ FLAG_APPROACH,
+ FLAG_MINMETHOD,
+ FLAG_ANHYLAW,
+ FLAG_WARNING,
+ TOLERANCE_JS,
+ TOLERANCE_0,
+ TOLERANCE_NR,
+ MAX_ITER_NR,
+ TOLERANCE_OM,
+ MAX_ITER_OM,
+ FLAG_ANA,
+ TOLERANCE_NJ,
+ DELTA_0,
+ FLAG_HOMO
+ };
+ ElseIf (N==5)
param_EnergHyst={ dim, N,
Ja, ha, Jb, hb,
- w_1, chi_1,
- w_2, chi_2,
- w_3, chi_3,
- w_4, chi_4,
- w_5, chi_5,
+ w_1, kappa_1,
+ w_2, kappa_2,
+ w_3, kappa_3,
+ w_4, kappa_4,
+ w_5, kappa_5,
FLAG_INVMETHOD,
- FLAG_ROOTFINDING1D,
- FLAG_VARORDIFF,
+ FLAG_JACEVAL,
+ FLAG_APPROACH,
FLAG_MINMETHOD,
- FLAG_TANORLANG,
+ FLAG_ANHYLAW,
FLAG_WARNING,
TOLERANCE_JS,
TOLERANCE_0,
@@ -252,123 +249,89 @@ hb = 399.474553851
DELTA_0,
FLAG_HOMO
};
+ Else
+ Printf["N=%g is not valid",N];
+ Error["The parameters (w_k, kappa_k) for the asked number of cells (N) is not written yet in 'param_EnergHyst.pro'"];
EndIf
-
//*******************************************************************
// DISPLAY PARAMETERS SETTING
//*******************************************************************
Label_FLAG_INVMETHOD="...";
-Label_FLAG_ROOTFINDING1D="...";
-Label_FLAG_TANORLANG="...";
-Label_FLAG_VARORDIFF="...";
+Label_FLAG_JACEVAL="...";
+Label_FLAG_ANHYLAW="...";
+Label_FLAG_APPROACH="...";
Label_FLAG_MINMETHOD="...";
Label_FLAG_ANA="...";
If (FLAG_INVMETHOD == 1)
- Label_FLAG_INVMETHOD="FLAG_INVMETHOD=1 --> NR_ana (homemade)";
+ Label_FLAG_INVMETHOD="FLAG_INVMETHOD=1 --> NR_ana (homemade)";
EndIf
If (FLAG_INVMETHOD == 2)
- Label_FLAG_INVMETHOD="FLAG_INVMETHOD=2 --> NR_num (homemade)";
+ Label_FLAG_INVMETHOD="FLAG_INVMETHOD=2 --> NR_num (homemade)";
EndIf
If (FLAG_INVMETHOD == 3)
- Label_FLAG_INVMETHOD="FLAG_INVMETHOD=3 --> bfgs (homemade)";
- EndIf
- If (FLAG_INVMETHOD == 4)
- Label_FLAG_INVMETHOD="FLAG_INVMETHOD=4 --> hybridsj (gsl)";
- EndIf
- If (FLAG_INVMETHOD == 5)
- Label_FLAG_INVMETHOD="FLAG_INVMETHOD=5 --> hybridj (gsl)";
- EndIf
- If (FLAG_INVMETHOD == 6)
- Label_FLAG_INVMETHOD="FLAG_INVMETHOD=6 --> newton (gsl)";
- EndIf
- If (FLAG_INVMETHOD == 7)
- Label_FLAG_INVMETHOD="FLAG_INVMETHOD=7 --> gnewton (gsl)";
- EndIf
- If(FLAG_ROOTFINDING1D == 0 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=0 --> not done";
- EndIf
- If(FLAG_ROOTFINDING1D == 1 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=1 --> (1D minimization) golden section (gsl)";
- EndIf
- If(FLAG_ROOTFINDING1D == 2 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=2 --> (1D minimization) brent (gsl)";
- EndIf
- If(FLAG_ROOTFINDING1D == 3 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=3 --> (1D minimization) quad golden (gsl)";
- EndIf
- If(FLAG_ROOTFINDING1D == 4 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=4 --> (1D root bracketing) bisection (gsl)";
- EndIf
- If(FLAG_ROOTFINDING1D == 5 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=5 --> (1D root bracketing) brent (gsl)";
- EndIf
- If(FLAG_ROOTFINDING1D == 6 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=6 --> (1D root bracketing) falsepos (gsl)";
- EndIf
- If(FLAG_ROOTFINDING1D == 7 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=7 --> (1D root finding with derivatives) newton (gsl)";
+ Label_FLAG_INVMETHOD="FLAG_INVMETHOD=3 --> bfgs (homemade)";
EndIf
- If(FLAG_ROOTFINDING1D == 8 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=8 --> (1D root finding with derivatives) secant (gsl)";
+ If(FLAG_JACEVAL == 1)
+ Label_FLAG_JACEVAL="FLAG_JACEVAL=1 --> analytical Jacobian";
EndIf
- If(FLAG_ROOTFINDING1D == 9 && FLAG_INVMETHOD != 0)
- Label_FLAG_ROOTFINDING1D="FLAG_ROOTFINDING1D=9 --> (1D root finding with derivatives) steffenson (gsl)";
+ If(FLAG_JACEVAL == 2)
+ Label_FLAG_JACEVAL="FLAG_JACEVAL=2 --> numerical Jacobian";
EndIf
- If (FLAG_TANORLANG == 1)
- Label_FLAG_TANORLANG="FLAG_TANORLANG=1 --> tangent hyperbolic";
+ If (FLAG_ANHYLAW == 1)
+ Label_FLAG_ANHYLAW="FLAG_ANHYLAW=1 --> tangent hyperbolic";
EndIf
- If (FLAG_TANORLANG == 2)
- Label_FLAG_TANORLANG="FLAG_TANORLANG=2 --> double langevin function";
+ If (FLAG_ANHYLAW == 2)
+ Label_FLAG_ANHYLAW="FLAG_ANHYLAW=2 --> double langevin function";
EndIf
- If (FLAG_VARORDIFF == 1)
- Label_FLAG_VARORDIFF="FLAG_VARORDIFF=1 --> variational approach (Jk)";
+ If (FLAG_APPROACH == 1)
+ Label_FLAG_APPROACH="FLAG_APPROACH=1 --> variational approach (Jk)";
EndIf
- If (FLAG_VARORDIFF == 2)
- Label_FLAG_VARORDIFF="FLAG_VARORDIFF=2 --> simple differential approach (hrk)";
+ If (FLAG_APPROACH == 2)
+ Label_FLAG_APPROACH="FLAG_APPROACH=2 --> vector play model approach (hrk)";
Label_FLAG_MINMETHOD ="...";
EndIf
- If (FLAG_VARORDIFF == 3)
- Label_FLAG_VARORDIFF="FLAG_VARORDIFF=2 --> full differential approach (hrk)";
+ If (FLAG_APPROACH == 3)
+ Label_FLAG_APPROACH="FLAG_APPROACH=3 --> full differential approach (hrk)";
Label_FLAG_MINMETHOD ="...";
EndIf
- If (FLAG_MINMETHOD == 1 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 1 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=1 --> steepest descent (homemade)";
EndIf
- If (FLAG_MINMETHOD == 11 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 11 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11 --> NEW steepest descent naive (homemade)";
EndIf
- If (FLAG_MINMETHOD == 22 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 22 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11 --> NEW conjugate Fletcher-Reeves (homemade)";
EndIf
- If (FLAG_MINMETHOD == 33 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 33 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11 --> NEW conjugate Polak-Ribiere (homemade)";
EndIf
- If (FLAG_MINMETHOD == 333 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 333 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11 --> NEW conjugate Polak-Ribiere+ (homemade)";
EndIf
- If (FLAG_MINMETHOD == 1999 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 1999 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11 --> NEW conjugate Dai Yuan 1999 (p.85) (homemade)";
EndIf
- If (FLAG_MINMETHOD == 2005 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 2005 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11 --> NEW conjugate Hager Zhang 2005 (p.161) (homemade)";
EndIf
- If (FLAG_MINMETHOD == 77 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 77 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=11 --> NEW newton (homemade)";
EndIf
- If (FLAG_MINMETHOD == 2 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 2 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=2 --> conjugate fr (gsl)";
EndIf
- If (FLAG_MINMETHOD == 3 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 3 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=3 --> conjugate pr (gsl)";
EndIf
- If (FLAG_MINMETHOD == 4 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 4 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=4 --> bfgs2 (gsl)";
EndIf
- If (FLAG_MINMETHOD == 5 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 5 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=5 --> bfgs (gsl)";
EndIf
- If (FLAG_MINMETHOD == 6 && FLAG_VARORDIFF==1)
+ If (FLAG_MINMETHOD == 6 && FLAG_APPROACH==1)
Label_FLAG_MINMETHOD ="FLAG_MINMETHOD=6 --> steepest descent (gsl)";
EndIf
If (FLAG_ANA ==0)
@@ -383,16 +346,16 @@ Label_FLAG_ANA="...";
Printf["Number of cells --> %g",N];
Printf["k omega kappa "];
For iSub In {1:N}
- Printf["%g %g %g ",iSub,w~{iSub},chi~{iSub}];
+ Printf["%g %g %g ",iSub,w~{iSub},kappa~{iSub}];
EndFor
Printf[""];
- Printf[StrCat[Label_FLAG_TANORLANG]];
+ Printf[StrCat[Label_FLAG_ANHYLAW]];
Printf["Ja ha Jb hb"];
Printf["%g %g %g %g",Ja,ha,Jb,hb];
Printf[""];
- Printf[StrCat[Label_FLAG_VARORDIFF]];
+ Printf[StrCat[Label_FLAG_APPROACH]];
Printf[StrCat[Label_FLAG_INVMETHOD]];
- Printf[StrCat[Label_FLAG_ROOTFINDING1D]];
+ Printf[StrCat[Label_FLAG_JACEVAL]];
Printf[StrCat[Label_FLAG_MINMETHOD]];
Printf[StrCat[Label_FLAG_ANA]];
Printf["FLAG_HOMO --> %g",FLAG_HOMO];
@@ -402,7 +365,7 @@ Label_FLAG_ANA="...";
Printf[""];
Printf["TOLERANCE_NR --> %g",TOLERANCE_NR];
Printf["MAX_ITER_NR --> %g",MAX_ITER_NR];
- If(FLAG_VARORDIFF==1)
+ If(FLAG_APPROACH==1)
Printf["TOLERANCE_NJ --> %g",TOLERANCE_NJ];
Printf["DELTA_0 --> %g",DELTA_0];
Printf[""];
diff --git a/Team32/square.geo b/Team32/square.geo
old mode 100644
new mode 100755
index c4309cb1ddd9b178d3fa3243d8b9cf88acf131e9..6d9824410cf2edb7d2c9437ad36c9689b838081b
--- a/Team32/square.geo
+++ b/Team32/square.geo
@@ -36,6 +36,7 @@ Physical Line(15) = {1,2,3,4};
Physical Surface(16)={6};
+Recursive Color LightGrey { Surface{ 6}; }
// Post-processing point
@@ -51,3 +52,4 @@ For k In {1:num_postop_points}
EndFor
EndIf
+
diff --git a/Team32/square.pro b/Team32/square.pro
index ae2425702fd62a0cde5580083af3e83056e191b0..d31484d514f808a11247caae0ea03fbf8c718e9c 100644
--- a/Team32/square.pro
+++ b/Team32/square.pro
@@ -24,7 +24,7 @@ getdp square.pro -pos Get_AllTS
// To launch computation and Then
// do PostOp after the end of the simulation (create res/*.pos res/*dat)
-getdp square.pro -solve Analysis -pos Get_AllTS
+getdp square.pro -v 3 -solve Analysis -pos Get_AllTS
//---------
@@ -67,7 +67,11 @@ Include "square_data.geo";
// Output Directory
//-------------------------------------------------------------------------
-Dir = 'res/';
+If (Exists(simulabel))
+ Dir = Sprintf("res%g/",simulabel);
+Else
+ Dir = 'res/';
+EndIf
//Include "param_EnergHyst.dat";
//Dir = Sprintf("res_N%g_C%g/",N00,Flag_TestCase_00);
@@ -131,8 +135,8 @@ Function {
hx[] = Frelax[] * ( hax * Cos[2*Pi*Freq *$Time]
+ haharmx * Cos[2*Pi*freqharm*$Time]) ;
hy[] = Frelax[] * ( hay0
- + hay * (Cos[phase]*Cos[2*Pi*Freq *$Time] + Sin[phase]*Sin[2*Pi*Freq *$Time])
- + haharmy * (Cos[phase]*Cos[2*Pi*freqharm*$Time] + Sin[phase]*Sin[2*Pi*freqharm*$Time]) );
+ + hay * (Cos[phase]*Cos[2*Pi*Freq *$Time] - Sin[phase]*Sin[2*Pi*Freq *$Time])
+ + haharmy * (Cos[phase]*Cos[2*Pi*freqharm*$Time] - Sin[phase]*Sin[2*Pi*freqharm*$Time]) );
/*
hx[] = Frelax[] * ( hax * Complex_MH[1,0]{Freq}
+ haharmx * Complex_MH[1,0]{freqharm}) ;
diff --git a/Team32/square_data.dat b/Team32/square_data.dat
new file mode 100644
index 0000000000000000000000000000000000000000..fef67df03f4e109c7b35238c2d0639251a03ddf8
--- /dev/null
+++ b/Team32/square_data.dat
@@ -0,0 +1,35 @@
+simulabel = 0;
+THICKNESS_00 = 1;
+NumEle_00 = 1;
+mm = 0.001;
+Flag_3Dmodel = 0;
+Flag_AnalysisType_00 = 1;
+freq_00 = 1;
+nstep_00 = 200;
+NbT_00 = 2;
+tmax_00 = 1;
+ha_00 = 150;
+Flag_TestCase_00 = 0;
+hax_00 = 150;
+haharmx_00 = 0;
+hay0_00 = 0;
+hay_00 = 150;
+haharmy_00 = 0;
+freqharm_00 = 0;
+phase_00 = (20/180)*3.141592653589793;
+Flag_NL_law00 = 4;
+Flag_NLRes00 = 1;
+Flag_UpdateSeparated= 1;
+Flag_ItLoopNRes00 = 31;
+Flag_ItLoopNDoFirstIter00 = 0;
+Nb_max_iter00 = 50;
+stop_criterion00 = 1e-05;
+relaxation_factor00 = 1;
+Flag_AdaptRelax00 = 1;
+relax_max00 = 1;
+relax_min00 = 0.1;
+relax_numtest00 = 10;
+TestAllFactors00 = 1;
+Flag_AdaptRelax_Version = 2;
+Flag_ExtrapolationOrder = 0;
+
diff --git a/Team32/square_data.geo b/Team32/square_data.geo
index 801c10f99c002ffce3a41a125550d66cdef3163c..0df6279bcf99d7f2ce62879f51f41f261266847b 100644
--- a/Team32/square_data.geo
+++ b/Team32/square_data.geo
@@ -1,6 +1,11 @@
//-------------------------------------------------------------------------
// Default Parameters
//-------------------------------------------------------------------------
+LOAD_DATA_FILE=0;
+If (LOAD_DATA_FILE && (FileExists("square_data.dat")))
+ Printf('square_data.dat is used');
+ Include "square_data.dat";
+Else
// GEO parameters .............................
THICKNESS_00 = 1; // thickness of the sample
@@ -14,13 +19,13 @@ Flag_AnalysisType_00 = 1; // 0: "static"
// 1: "time domain"
// 2: "frequency domain"
freq_00 = 1; // frequency (no impact as there is no dynamic effect)
-nstep_00 = 200; // number of time step per period (without Flexible DT)
-NbT_00 = 2; //number of periods
+nstep_00 = 100; // number of time step per period (without Flexible DT)
+NbT_00 = 1; //number of periods
tmax_00 = 1; // final time of simulation
// SOURCE parameters ............................
ha_00 = 150; // maximal amplitude of the imposed magnetic field
-Flag_TestCase_00 = 1; // 0: CASE 0 (Default)
+Flag_TestCase_00 = 6; // 0: CASE 0 (Default)
// 1: CASE 1 (1Dx)
// 2: CASE 2 (1Dx+harm)
// 3: CASE 3 (1Dxy)
@@ -45,25 +50,44 @@ Flag_NL_law00 = 4; // 0: "linear",
// 3: "Jiles-Atherton hysteresis model",
// 4: "EnergHyst model"
-Flag_NLRes00 = 0; // 0: use classical IterativeLoop to solve the NL (non linear) problem
+Flag_ExtrapolationOrder = 0; // degree 0 ==> initialization at the previous Time Solution
+ // degree 1 ==> linear extrapolation from the 2 previous Time Solutions
+ // degree 2 ==> quadratic extrapolation from the 3 previous Time Solutions
+
+Flag_NLRes00 = 2; // 0: use classical IterativeLoop to solve the NL (non linear) problem
// 1: use IterativeLoopN to solve the NL (non linear) problem
// 2: solve the NL (non linear) problem "by hand"
// non linear loop default parameters
+Flag_UpdateSeparated= 1;
Nb_max_iter00 = 50; // maximum number of NL (non linear) iterations
stop_criterion00 = 1e-5; // strop criterion for the NL (non linear) iteration
relaxation_factor00 = 1; // default relaxation factor (between ]0,1] )
Flag_AdaptRelax00 = 1; // set to 1 to test different relaxation factors
+Flag_AdaptRelax_Version = 2; // allow to chose between MySolveJac_AdaptRelax or MySolveJac_AdaptRelax2 in 'MyItLoopPro.pro'
+ // (used when Flag_NLRes00==2 and Flag_AdaptRelax00==1)
relax_max00 = 1 ;
relax_min00 = 0.1;
relax_numtest00 = 10;
-TestAllFactors00 = 2; // 0 : try first relaxation factors (from the list) and stops when the residual goes up
+TestAllFactors00 = 1; // 0 : try first relaxation factors (from the list) and stops when the residual goes up
// 1 : try every relaxation factors (from the list) and keep the optimal one
// 2 : find the maximum relaxation factor that decreases the residual:
// - start with the relaxation factor from the previous time step or from the previous iteration
// - the relaxation factor is multiplied by a ratio as long as the residual decreases
// - the relaxation factor is decreased by a ratio until a decreasing residual is found
// [3 : Build a parabola based on three relaxation factors and deduce a minimizing relaxation factor (NOT WORKING!!)]
+Flag_ItLoopNRes00 = 1;// 0 : Solution (dx=x-xp; x=x)
+ // 1 : Residual (dx=res=b-Ax; x=b) (default for square ???)
+ // 2 : RecalcResidual (dx=res=b-Ax; x=b)
+ // 31 : PostOp unknown field (az or f)
+ // 32 : PostOp b field
+ // 33 : PostOp h field (default for t32)
+Flag_ItLoopNDoFirstIter00 = 0;// set to 1 if one wants to force the first NL iteration before evaluating the error with IterativeLoopN
+ // helpfull for square.pro with EnergHyst or JA because PostOp h field = 0 at iter 1, which stops IterativeLoopN prematurely
+
+
+EndIf
+
Reltol_Mag00 = stop_criterion00; // 0 before with IterativeLoopN
Abstol_Mag00 = stop_criterion00;
Reltoldx_Mag00 = 1e-5*stop_criterion00;
@@ -86,6 +110,13 @@ NLRES_ITERATIVELOOP = 0;
NLRES_ITERATIVELOOPN = 1;
NLRES_ITERATIVELOOPPRO = 2;
+NLITLOOPN_SOLUTION = 0;
+NLITLOOPN_RESIDUAL = 1;
+NLITLOOPN_RECALCRESIDUAL = 2;
+NLITLOOPN_POSTOPX = 31;
+NLITLOOPN_POSTOPB = 32;
+NLITLOOPN_POSTOPH = 33;
+
AN_STATIC = 0;
AN_TIME = 1;
AN_FREQUENCY = 2;
@@ -183,6 +214,23 @@ DefineConstant[
Help Str["- Use 'IterativeLoop' to solve the NL (non linear) problem with classical IterativeLoop Operation",
"- Use 'IterativeLoopN' to solve the NL (non linear) problem with IterativeLoopN Operation",
"- Use 'IterativeLoopPro' to solve the NL (non linear) problem 'by hand'"]},
+ Flag_ItLoopNRes = { (Flag_NLRes!=NLRES_ITERATIVELOOPN)? NLITLOOPN_SOLUTION:Flag_ItLoopNRes00 , Choices {
+ NLITLOOPN_SOLUTION ="Solution",
+ NLITLOOPN_RESIDUAL ="RecalcResidual",
+ NLITLOOPN_RECALCRESIDUAL ="Residual",
+ NLITLOOPN_POSTOPX ="PostOp Potential",
+ NLITLOOPN_POSTOPB ="PostOp b field",
+ NLITLOOPN_POSTOPH ="PostOp h field"},
+ Name StrCat[ppNL, "1Error Calculation based on ..."], Highlight Str[colNL1], Visible Flag_NL ,ReadOnly (Flag_NLRes!=NLRES_ITERATIVELOOPN),
+ Help Str["- the solution (dx=x-xp; x=x)",
+ "- the residual (dx=res=b-Ax; x=b)",
+ "- the recalculated residual (dx=res=b-Ax; x=b)",
+ "- the post-operated potential field (az or phi)",
+ "- the post-operated flux density field (b)",
+ "- the post-operated magnetic field (h)"]},
+ Flag_ItLoopNDoFirstIter = {Flag_ItLoopNDoFirstIter00, Choices{0,1}, Name StrCat[ppNL, "1First iteration done before error calculation"],
+ Highlight Str[colNL1], Visible (Flag_NL && Flag_NLRes==NLRES_ITERATIVELOOPN),
+ Help Str["Set to 1 to force the first NL iteration to be done before evaluating the error with IterativeLoopN"]},
Nb_max_iter = {Nb_max_iter00, Name StrCat[ppNL, "1Nb max iter"],
Highlight Str[colNL4], Visible Flag_NL}
diff --git a/Team32/t32.geo b/Team32/t32.geo
index 7473533284098b6ec64ef90d854aed3e652219e7..bf922e910642404e343571d6750871b4394f0ea3 100644
--- a/Team32/t32.geo
+++ b/Team32/t32.geo
@@ -19,10 +19,10 @@ lc2 = wcore/5;
/*
-lc0=lc0/2;
-lc1=lc1/2;
-lc2=lc2/2;
-*/
+lc0=lc0/3;
+lc1=lc1/3;
+lc2=lc2/3;
+//*/
lcri = lc0*4; // Pi*Rint/30;
lcro = lc0*4; // Pi*Rext/30;
@@ -271,32 +271,32 @@ If(Flag_3Dmodel==0)
//=================================================
// Physical regions for FE analysis (2D)
//=================================================
- Physical Surface("Core", CORE) = surf_ECore[];
+ Physical Surface(CORE) = surf_ECore[];
nb_surf_coil = #surf_Coil[];
- Physical Surface("Coil right (left side)", COIL) = surf_Coil[{0:nb_surf_coil-1:4}];
- Physical Surface("Coil right (right side)", COIL+1) = surf_Coil[{1:nb_surf_coil-1:4}];
- Physical Surface("Coil left (right side)", COIL+2) = surf_Coil[{2:nb_surf_coil-1:4}];
- Physical Surface("Coil left (left side)", COIL+3) = surf_Coil[{3:nb_surf_coil-1:4}];
+ Physical Surface(COIL) = surf_Coil[{0:nb_surf_coil-1:4}];
+ Physical Surface(COIL+1) = surf_Coil[{1:nb_surf_coil-1:4}];
+ Physical Surface(COIL+2) = surf_Coil[{2:nb_surf_coil-1:4}];
+ Physical Surface(COIL+3) = surf_Coil[{3:nb_surf_coil-1:4}];
- Physical Surface("Air", AIR) = surf_Air[];
- Physical Surface("Air infinity shell", AIRINF) = surf_AirInf[];
- Physical Line("Outer boundary", SURF_AIROUT) = ln_bnd[];
+ Physical Surface(AIR) = surf_Air[];
+ Physical Surface(AIRINF) = surf_AirInf[];
+ Physical Line(SURF_AIROUT) = ln_bnd[];
- //=================================================
+ //=================================================
// Some colors... for aesthetics :-)
//=================================================
Reverse Surface {surf_ECore[{2,3}], surf_Coil[{2,3}],surf_Air[{2,3}],surf_AirInf[1]};
- Recursive Color SkyBlue { Surface{surf_Air[]};}
- Recursive Color Blue { Surface{surf_AirInf[]};}
- Recursive Color SteelBlue { Surface{ surf_ECore[]}; }
- Recursive Color Red { Surface{surf_Coil[{0,1}]}; }
- Recursive Color Yellow { Surface{surf_Coil[{2,3}]}; }
+ Recursive Color White { Surface{surf_Air[]};}
+ Recursive Color White { Surface{surf_AirInf[]};}
+ Recursive Color LightGrey { Surface{ surf_ECore[]}; }
+ Recursive Color LightBlue { Surface{surf_Coil[{0,1}]}; }
+ Recursive Color Pink { Surface{surf_Coil[{2,3}]}; }
If(!Flag_Symmetry)
Reverse Surface {surf_ECore[{4,5}], surf_Coil[{4,5}],surf_Air[{4,5}],surf_AirInf[2]};
- Recursive Color Red { Surface{surf_Coil[{4,5}]}; }
- Recursive Color Yellow { Surface{surf_Coil[{6,7}]}; }
+ Recursive Color LightBlue { Surface{surf_Coil[{4,5}]}; }
+ Recursive Color Pink { Surface{surf_Coil[{6,7}]}; }
EndIf
EndIf
@@ -316,13 +316,13 @@ If(Flag_3Dmodel==1)
Physical Surface("Skin Core", SKINCORE) = bnd_Core[];
Physical Volume("Coil right", COIL+0) = vol_CoilR[];
- Physical Volume("Coil left", COIL+1) = vol_CoilL[];
+ Physical Volume("Coil left",COIL+1) = vol_CoilL[];
- Physical Volume("Leg in Coil right",LEG_INCOIL+0) = vol_inCoilR[]; // for source computation
+ Physical Volume("Leg in Coil right", LEG_INCOIL+0) = vol_inCoilR[]; // for source computation
Physical Volume("Leg in Coil left", LEG_INCOIL+1) = vol_inCoilL[];
nb_surf_coil = #surf_Coil[];
- Physical Surface("Electrode Coil right",ELECCOIL+0) = surf_Coil[{1:nb_surf_coil-1:4}]; // e.g. the one of the right
+ Physical Surface("Electrode Coil right", ELECCOIL+0) = surf_Coil[{1:nb_surf_coil-1:4}]; // e.g. the one of the right
Physical Surface("Electrode Coil left", ELECCOIL+1) = surf_Coil[{2:nb_surf_coil-1:4}];
//Right
@@ -336,13 +336,13 @@ If(Flag_3Dmodel==1)
Physical Surface("Skin Coil right", SKINCOIL+0) = bnd_CoilR[];
bnd_CoilR_plus_hole[] = Abs(CombinedBoundary{Volume{vol_CoilR[], vol_inCoilR[]};});
- bnd_inCoilR[] = Abs(CombinedBoundary{Volume{vol_inCoilR[]};});
+ bnd_inCoilR[] = Abs(CombinedBoundary{Volume{vol_inCoilR[]};});
bnd_inCoilR[] -= bnd_CoilR_plus_hole[];
-Physical Surface("Skin Coil right- only around cut", SKINCOIL_CUT+0) = bnd_inCoilR[];
+ Physical Surface("Skin Coil right - only cut",SKINCOIL_CUT+0) = bnd_inCoilR[];
bnd_CoilR_plus_hole[] -= Abs(CombinedBoundary{Volume{vol_CoilR[]};});
bnd_CoilR_plus_hole[] -= surf_cut_xy[];
- Physical Surface("Cut Coil right", COIL_CUT+0)= bnd_CoilR_plus_hole[1]; //cut down (top at 0)
+ Physical Surface("Cut Coil right",COIL_CUT+0)= bnd_CoilR_plus_hole[1]; //cut down (top at 0)
If (Flag_Symmetry3D ==1)
surf_cut_xz[] += bnd_CoilR_plus_hole[1];
@@ -356,25 +356,25 @@ Physical Surface("Skin Coil right- only around cut", SKINCOIL_CUT+0) = bnd_inCoi
surf_cut_xz[] += bnd_CoilL[{4,6,9,12,16}];
bnd_CoilL[] -= surf_cut_xz[];
EndIf
- Physical Surface("Skin Coil left", SKINCOIL+1) = bnd_CoilL[];
+ Physical Surface("Skin Coil left", SKINCOIL+1) = bnd_CoilL[];
bnd_CoilL_plus_hole[] = Abs(CombinedBoundary{Volume{vol_CoilL[], vol_inCoilL[]};});
- bnd_inCoilL[] = Abs(CombinedBoundary{Volume{vol_inCoilL[]};});
+ bnd_inCoilL[] = Abs(CombinedBoundary{Volume{vol_inCoilL[]};});
bnd_inCoilL[] -= bnd_CoilL_plus_hole[];
- Physical Surface("Skin Coil left- only around cut", SKINCOIL_CUT+1) = bnd_inCoilL[];
+ Physical Surface("Skin Coil left - only cut",SKINCOIL_CUT+1) = bnd_inCoilL[];
bnd_CoilL_plus_hole[] -= Abs(CombinedBoundary{Volume{vol_CoilL[]};});
bnd_CoilL_plus_hole[] -= surf_cut_xy[];
- Physical Surface("Cut Coil left", COIL_CUT+1)= bnd_CoilL_plus_hole[1]; //cut down (top at 0)
+ Physical Surface("Cut Coil left",COIL_CUT+1)= bnd_CoilL_plus_hole[1]; //cut down (top at 0)
If (Flag_Symmetry3D ==1)
surf_cut_xz[] += bnd_CoilL_plus_hole[1];
EndIf
- Physical Volume("Air", AIR) = vol_Air[];
- Physical Surface("Symmetry cut XY", SURF_CUTXY) = surf_cut_xy[];
- Physical Surface("Symmetry cut XZ", SURF_CUTXZ) = surf_cut_xz[];
+ Physical Volume("Air",AIR) = vol_Air[];
+ Physical Surface("cut XY", SURF_CUTXY) = surf_cut_xy[];
+ Physical Surface("cut XZ", SURF_CUTXZ) = surf_cut_xz[];
all_vol[] = Volume '*';
bnd_all[] = Abs(CombinedBoundary{Volume{all_vol[]};});
@@ -385,14 +385,14 @@ Physical Surface("Skin Coil right- only around cut", SKINCOIL_CUT+0) = bnd_inCoi
//=================================================
// Some colors... for aesthetics :-)
//=================================================
- Recursive Color SkyBlue { Volume{vol_Air[]};}
- Recursive Color SteelBlue { Volume{ vol_Core[]}; }
+ Recursive Color White { Volume{vol_Air[]};}
+ Recursive Color LightGrey { Volume{ vol_Core[]}; }
Recursive Color Red { Volume{vol_CoilR[]}; }
Recursive Color Yellow { Volume{vol_CoilL[]}; }
EndIf
-
+aaa=newp;
// Post-processing point
For k In {1:num_postop_points}
Point(newp) = {xpos~{k},ypos~{k}, 0 }; // for visu
diff --git a/Team32/t32.pro b/Team32/t32.pro
index d5811f0bf72d7822a9242c2e6bdc688c1ecf58c0..25d561e79bd852e0a5a82e0bab8c5ca9f8cc69c9 100644
--- a/Team32/t32.pro
+++ b/Team32/t32.pro
@@ -25,6 +25,7 @@ getdp t32.pro -pos Get_AllTS
// To launch computation and Then
// do PostOp after the end of the simulation (create res/*.pos res/*dat)
getdp t32.pro -solve Analysis -pos Get_AllTS
+getdp t32.pro -solve Analysis -pos Get_AllTS_3D
//---------
@@ -157,7 +158,11 @@ Function {
xc2_2 = -xc1_2;
SurfCoil[] = SurfaceArea[]{ELECCOIL+0} ;
-
+/*
+ vDir[Ind_2] = Vector[0,0,1];
+ vDir[Ind_1] = Vector[0,0,-1];
+*/
+///*
vDir[Ind_2] = (Fabs[Z[]]>=Lz/2 && Fabs[X[]-x_1]<=wcore/2 ) ? Vector[-1,0,0]:(
(Fabs[Z[]]<Lz/2)? Vector[0, 0, (X[]>x_1)?-1:1 ]:
-Vector[ -Sin[Atan2[Z[]-Lz/2, X[]-((X[]>x_1)?xc1_1:xc1_2) ]#1], 0, Cos[#1]] ) ;
@@ -165,7 +170,7 @@ Function {
vDir[Ind_1] = (Fabs[Z[]]>=Lz/2 && Fabs[X[]-x_2]<=wcore/2 ) ? Vector[-1,0,0]:(
(Fabs[Z[]]<Lz/2)? Vector[0, 0, (X[]>x_2)?-1:1 ]:
-Vector[ -Sin[Atan2[Z[]-Lz/2, X[]-((X[]>x_2)?xc2_2:xc2_1) ]#1], 0, Cos[#1]] ) ;
-
+//*/
EndIf
js0[] = NbWires[]/SurfCoil[]*vDir[] ;
@@ -215,8 +220,8 @@ Function {
//pA = (Flag_AnalysisType==AN_STATIC) ? Pi/2: 0.; // Not Used
//IA[] = F_Sin_wt_p[]{2*Pi*Freq, pA}; // DomainB
- I1[] = ((Flag_NL_law==NL_JA || Flag_NL_law==NL_ENERGHYST)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_1}+V5[]);
- I2[] = ((Flag_NL_law==NL_JA || Flag_NL_law==NL_ENERGHYST)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_2}+V5[]);
+ I1[] = ((Flag_AnalysisType==AN_TIME)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_1}+V5[]);
+ I2[] = ((Flag_AnalysisType==AN_TIME)? Frelax[]:1) * (F_Sin_wt_p[]{2*Pi*Freq, Phase_2}+V5[]);
}
diff --git a/Team32/t32_data.geo b/Team32/t32_data.geo
index bd2eb7d0a0558732ec97c850ade28dd8034e4b7e..0b014c2879523660cd5384e09c4d5ea717dba3f6 100644
--- a/Team32/t32_data.geo
+++ b/Team32/t32_data.geo
@@ -27,20 +27,30 @@ Flag_TestCase00 = 3; //1: - case 1: windings series connected with sinusoidal su
//4: - case 4: secondary winding on load
// RESOLUTION parameters .........................
-Flag_NL_law00 = 4; // 0: "linear",
+Flag_NL_law00 = 1; // 0: "linear",
// 1: "analytical",
// 2: "anhysteretic part of JA-model",
// 3: "Jiles-Atherton hysteresis model",
// 4: "EnergHyst model"
+Flag_ExtrapolationOrder = 2; // degree 0 ==> initialization at the previous Time Solution
+ // degree 1 ==> linear extrapolation from the 2 previous Time Solutions
+ // degree 2 ==> quadratic extrapolation from the 3 previous Time Solutions
+
// non linear loop default parameters
Flag_NLRes00 = 1; // 0: use classical IterativeLoop to solve the NL (non linear) problem
// 1: use IterativeLoopN to solve the NL (non linear) problem
// 2: solve the NL (non linear) problem "by hand"
+
+Flag_DomainLam00= 0;
+
+Flag_UpdateSeparated= 1;
Nb_max_iter00 = 50; // maximum number of NL (non linear) iterations
stop_criterion00 = 1e-5; // strop criterion for the NL (non linear) iteration
relaxation_factor00 = 1; // default relaxation factor (between ]0,1] )
Flag_AdaptRelax00 = 1; // set to 1 to test different relaxation factors
+Flag_AdaptRelax_Version = 2; // allow to chose between MySolveJac_AdaptRelax or MySolveJac_AdaptRelax2 in 'MyItLoopPro.pro'
+ // (used when Flag_NLRes00==2 and Flag_AdaptRelax00==1)
relax_max00 = 1 ;
relax_min00 = 0.1;
relax_numtest00 = 10;
@@ -51,6 +61,15 @@ TestAllFactors00 = 1; // 0 : try first relaxation factors (from the list) an
// - the relaxation factor is multiplied by a ratio as long as the residual decreases
// - the relaxation factor is decreased by a ratio until a decreasing residual is found
// [3 : Build a parabola based on three relaxation factors and deduce a minimizing relaxation factor (NOT WORKING!!)]
+Flag_ItLoopNRes00 = 33;// 0 : Solution (dx=x-xp; x=x)
+ // 1 : Residual (dx=res=b-Ax; x=b) (default for square ???)
+ // 2 : RecalcResidual (dx=res=b-Ax; x=b)
+ // 31 : PostOp unknown field (az or f)
+ // 32 : PostOp b field
+ // 33 : PostOp h field (default for t32)
+Flag_ItLoopNDoFirstIter00 = 0;// set to 1 if one wants to force the first NL iteration before evaluating the error with IterativeLoopN
+ // helpfull for square.pro with EnergHyst or JA because PostOp h field = 0 at iter 1, which stops IterativeLoopN prematurely
+
Reltol_Mag00 = stop_criterion00; // 0 before with IterativeLoopN
Abstol_Mag00 = stop_criterion00;
Reltoldx_Mag00 = 1e-5*stop_criterion00;
@@ -73,6 +92,13 @@ NLRES_ITERATIVELOOP = 0;
NLRES_ITERATIVELOOPN = 1;
NLRES_ITERATIVELOOPPRO = 2;
+NLITLOOPN_SOLUTION = 0;
+NLITLOOPN_RESIDUAL = 1;
+NLITLOOPN_RECALCRESIDUAL = 2;
+NLITLOOPN_POSTOPX = 31;
+NLITLOOPN_POSTOPB = 32;
+NLITLOOPN_POSTOPH = 33;
+
AN_STATIC = 0;
AN_TIME = 1;
AN_FREQUENCY = 2;
@@ -129,6 +155,8 @@ DefineConstant[
//Printf("====> SymmetryFactor=%g", SymmetryFactor);
];
+ Printf("====> SymmetryFactor=%g", SymmetryFactor);
+
// Geometric dimensions........................ (subMenu ppDIM)
close_menu = 1;
DefineConstant[
@@ -144,7 +172,9 @@ DefineConstant[
Highlight Str[colorro]},
hcore = {hcoil+2*wcore, Name StrCat[ppDIM, "2Core height [m]"], ReadOnly 1,
Highlight Str[colorro]},
- AxialLength = {cteaux*2.4*mm, Name StrCat[ppDIM, "0Length along z-axis [m]"], Highlight Str[colorpp]}
+ AxialLength = {cteaux*2.4*mm, Name StrCat[ppDIM, "0Length along z-axis [m]"], Highlight Str[colorpp]},
+ Thickness = {0.48*mm, Name StrCat[ppDIM, "1Thickness of lamination [m]"],
+ Highlight Str[colorpp],Closed close_menu}
];
// Excitation Source .......................... (SubMenu ppTC)
@@ -183,6 +213,8 @@ DefineConstant[
Highlight Str[colAC2],Visible (Flag_AnalysisType==AN_TIME)}
NbSteps = {nstep_00, Name StrCat[ppAC, "Number of steps"],
Highlight Str[colAC2], Visible (Flag_AnalysisType==AN_TIME)}
+ Flag_DomainLam = {Flag_DomainLam00, Choices{0,1},
+ Name StrCat[ppAC, "1Activate DomainLam"], Visible 1}
];
// Material Law ............................... (SubMenu ppML)
@@ -209,6 +241,25 @@ DefineConstant[
Help Str["- Use 'IterativeLoop' to solve the NL (non linear) problem with classical IterativeLoop Operation",
"- Use 'IterativeLoopN' to solve the NL (non linear) problem with IterativeLoopN Operation",
"- Use 'IterativeLoopPro' to solve the NL (non linear) problem 'by hand'"]},
+ Flag_ItLoopNRes = { (Flag_NLRes!=NLRES_ITERATIVELOOPN)? NLITLOOPN_SOLUTION:Flag_ItLoopNRes00 , Choices {
+ NLITLOOPN_SOLUTION ="Solution",
+ NLITLOOPN_RESIDUAL ="RecalcResidual",
+ NLITLOOPN_RECALCRESIDUAL ="Residual",
+ NLITLOOPN_POSTOPX ="PostOp Potential",
+ NLITLOOPN_POSTOPB ="PostOp b field",
+ NLITLOOPN_POSTOPH ="PostOp h field"},
+ Name StrCat[ppNL, "1Error Calculation based on ..."], Highlight Str[colNL1], Visible Flag_NL ,ReadOnly (Flag_NLRes!=NLRES_ITERATIVELOOPN),
+ Help Str["- the solution (dx=x-xp; x=x)",
+ "- the residual (dx=res=b-Ax; x=b)",
+ "- the recalculated residual (dx=res=b-Ax; x=b)",
+ "- the post-operated potential field (az or phi)",
+ "- the post-operated flux density field (b)",
+ "- the post-operated magnetic field (h)"]},
+ Flag_ItLoopNDoFirstIter = {Flag_ItLoopNDoFirstIter00, Choices{0,1}, Name StrCat[ppNL, "1First iteration done before error calculation"],
+ Highlight Str[colNL1], Visible (Flag_NL && Flag_NLRes==NLRES_ITERATIVELOOPN),
+ Help Str["Set to 1 to force the first NL iteration to be done before evaluating the error with IterativeLoopN"]},
+
+
Nb_max_iter = {Nb_max_iter00, Name StrCat[ppNL, "1Nb max iter"],
Highlight Str[colNL4], Visible Flag_NL}
stop_criterion = {stop_criterion00, Name StrCat[ppNL, "2stop criterion"],
@@ -276,7 +327,8 @@ mu0 = 4.e-7 * Pi ;
sigma_al = 3.72e7 ; // conductivity of aluminum [S/m]
sigma_cu = 5.77e7 ; // conductivity of copper [S/m]
-sigma_core = 2.5e7/1000;
+//sigma_core = 2.5e7/1000;
+sigma_core = 1.78e6; // for team 32
sigma_coil = 5.9e7;
mur_fe = 1000; // linear case