DianaICSolver: AdsColumn.py

import diana
try:
  N = 10;
  dmain=diana.GetDianaMain()

  # get Diana model manger
  mmanger=dmain.GetModelManager();

  # create model
  model=mmanger.CreateModel(diana.CAPE_CONTINUOUS, "Adscolumn.so");
  model.Initialize();
  eso=model.GetActiveESO();
  esopar=eso.GetParameters();
  esopar["n"].SetValue(N);
  ri=dmain.CreateReportingInterface("basic");
  ri.SetComponentName("data");
  sfactory = dmain.GetSolverFactory();

  # residuals with initial guess
  print "Initial guess"
  print "States:", eso.GetAllVariables();
  print "Derivatives:", eso.GetAllDerivatives();
  print "Residuals:", eso.GetAllResiduals();

  # create consistent initial conditions solver
  icsolver=sfactory.CreateSolver(diana.CAPE_IC, model, "icsolver.so");
  icsolver.Initialize();
  icpar=icsolver.GetParameters();
  icpar["VerboseLevel"].SetValue(1);
  icsolver.Solve();

  # residuals with initial guess
  print "Consistent initial conditions"
  print "States:", eso.GetAllVariables();
  print "Derivatives:", eso.GetAllDerivatives();
  print "Residuals:", eso.GetAllResiduals();

# create solver
  solver=sfactory.CreateSolver(diana.CAPE_DAE, model, "ida.so");
  solver.Initialize();
  solpar=solver.GetParameters();
  solver.SetReportingInterface(ri);
  ri.Add(solpar["T"]);

  for i in range(1,N+1):
        name = "c_a["+str(i)+"]";
        ri.Add(eso.GetStateVariables().ItemByName(name));
  for i in range(1,N+1):
        name = "c_b["+str(i)+"]";
        ri.Add(eso.GetStateVariables().ItemByName(name));
  for i in range(1,N+1):
        name = "q_a["+str(i)+"]";
        ri.Add(eso.GetStateVariables().ItemByName(name));
  for i in range(1,N+1):
        name = "q_b["+str(i)+"]";
        ri.Add(eso.GetStateVariables().ItemByName(name));

  esopar["cina"].SetValue(8);
  esopar["cinb"].SetValue(8);
  esopar["v"].SetValue(1);
  esopar["l"].SetValue(20);
  esopar["n"].SetValue(N);
  solpar["VerboseLevel"].SetValue(1);
  solpar["Tend"].SetValue(2);
  solver.Solve();
  esopar["cina"].SetValue(0);
  esopar["cinb"].SetValue(0);
  solpar["Tend"].SetValue(100);
  solver.Solve();

# create directory and store results
  import os
  if not os.path.isdir('results'):
          os.mkdir('results');
  ri.WriteData("results/AdsColumn");
  ri.WriteDataMatlab("results/AdsColumn.m");

except diana.ECapeUser, exc:
  print repr(exc)