DianaTutorial: steadystate.py

File steadystate.py, 2.5 KB (added by miha, 16 years ago)
Line 
1import diana
2try:
3    # get diana main and load Hafke reactor model
4    dmain=diana.GetDianaMain()
5    mmanger=dmain.GetModelManager();
6    model=mmanger.CreateModel(diana.CAPE_CONTINUOUS, "HafkeReactor.so");
7    model.Initialize();
8    eso=model.GetActiveESO();
9    # initialize model variables and parameters
10    esopar=eso.GetParameters();
11    eso.SetAllVariables([4.3796022e+01, 6.4482417e+01, 1.3579274e+03,
12                         5.6989992e-01, 3.6715772e+02, 3.2913593e+02]);
13    esopar["q"].SetValue(1.84e-06);
14    esopar["tkzu"].SetValue(2.9165e+02);
15    esopar["tzu"].SetValue(2.9365e+02);
16    esopar["c_zu"].SetValue([3.0e+03, 0.0e+0, 0.0e+0, 6.54e+0]);
17    esopar["c_f_ges"].SetValue(6.54e+0);
18    esopar["dhr"].SetValue([3.025e+02, 3.859e+02, 9.51e+01]);
19    esopar["e"].SetValue([1.055e+02, 1.262e+02, 8.833e+01,
20                          5.569e+01, 4.504e+01]);
21    esopar["k0"].SetValue([3.7e+13, 4.83e+14, 1.65e+10,
22                           3.72e+05, 3.8333e+04]);
23    esopar["rgas"].SetValue(8.314e-03);
24    esopar["kwfk"].SetValue(1.4e-01);
25    esopar["vk"].SetValue(8.0e-06);
26    esopar["rho_cp"].SetValue(4.106e+03);
27    esopar["vr"].SetValue(2.9e-03);
28    esopar["qknormal"].SetValue(8.0e-06);
29    # create steady state continuation solver
30    sfactory=dmain.GetSolverFactory();
31    conti=sfactory.CreateSolver(diana.CAPE_CONTI, model, "sstate.so");
32    conti.Initialize();
33    conpar=conti.GetParameters();
34    # create reporting interface and specify report variables
35    ri=dmain.CreateReportingInterface("basic");
36    conti.SetReportingInterface(ri);
37    ri.Add(esopar["qknormal"]);
38    ri.Add(eso.GetStateVariables().ItemByName("tk"));
39    ri.Add(conpar["Stability"]);
40    ri.Add(conpar["ConditionCurrent"]);
41    # set continuation parameters
42    conti.AddFreeParameter("qknormal", 0.5e-5, 5e-5);
43    conpar["VerboseLevel"].SetValue(3);
44    conpar["InitialDirection"].SetValue(1.0);
45    conpar["Parameterisation"].SetValue("PseudoArclength");
46    conpar["Predictor"].SetValue("Tangent");
47    conpar["InitialStepSize"].SetValue(0.1);
48    conpar["ConditionCheck"].SetValue(diana.SteadyStateNone);
49    conpar["StabilityCheck"].SetValue(True);
50    conpar["MaxStepSize"].SetValue(50.0);
51    conpar["MaxStepsAmount"].SetValue(2000);
52    # continuation
53    conti.Continuate();
54    # create directory and store continuation results
55    import os
56    if not os.path.isdir('results'): os.mkdir('results');
57    ri.WriteDataMatlab("results/steadystate.m");
58except diana.ECapeUser, exc:
59    print exc
60