instationary.py

from dune.codegen.generation import (include_file,
                                     preamble,
                                     )
from dune.codegen.pdelab.driver import (get_form_ident,
                                        get_trial_element,
                                        is_linear,
                                        name_initree,
                                        preprocess_leaf_data,
                                        )
from dune.codegen.pdelab.driver.gridfunctionspace import (name_trial_gfs,
                                                          type_range,
                                                          )
from dune.codegen.pdelab.driver.gridoperator import (name_gridoperator,
                                                     type_gridoperator,
                                                     name_localoperator,
                                                     )
from dune.codegen.pdelab.driver.constraints import (has_dirichlet_constraints,
                                                    name_bctype_grid_function,
                                                    name_constraintscontainer,
                                                    )
from dune.codegen.pdelab.driver.interpolate import (interpolate_dirichlet_data,
                                                    name_boundary_grid_function,
                                                    )
from dune.codegen.pdelab.driver.solve import (print_matrix,
                                              print_residual,
                                              name_linearsolver,
                                              name_stationarynonlinearproblemsolver,
                                              name_vector,
                                              type_linearsolver,
                                              type_stationarynonlinearproblemssolver,
                                              type_vector,
                                              )
from dune.codegen.pdelab.driver.vtk import (name_vtk_sequence_writer,
                                            visualize_initial_condition,
                                            name_predicate)
from dune.codegen.options import (get_form_option,
                                  get_option,
                                  )


def solve_instationary():
    # Create time loop
    if get_form_option('matrix_free'):
        raise NotImplementedError("Instationary matrix free not implemented!")
    else:
        time_loop()

    print_residual()
    print_matrix()


@preamble(section="instat")
def time_loop():
    ini = name_initree()
    lop = name_localoperator(get_form_ident())
    time = name_time()
    element = get_trial_element()
    vector_type = type_vector(get_form_ident())
    vector = name_vector(get_form_ident())
    interpolate_dirichlet_data(vector)
    gfs = name_trial_gfs()

    is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
    assemble_new_constraints = ""
    if has_dirichlet_constraints(is_dirichlet):
        bctype = name_bctype_grid_function(element, is_dirichlet)
        cc = name_constraintscontainer()
        assemble_new_constraints = ("  // Assemble constraints for new time step\n"
                                    "  {}.setTime({}+dt);\n"
                                    "  Dune::PDELab::constraints({}, {}, {});\n"
                                    "\n".format(lop, time, bctype, gfs, cc)
                                    )

    # Choose between explicit and implicit time stepping
    explicit = get_option('explicit_time_stepping')
    if explicit:
        osm = name_explicitonestepmethod()
        apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
    else:
        osm = name_onestepmethod()
    if has_dirichlet_constraints(is_dirichlet):
        dirichlet = preprocess_leaf_data(element, "interpolate_expression")
        boundary = name_boundary_grid_function(element, dirichlet)
        apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
    else:
        apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)

    # Setup visualization
    visualize_initial_condition()
    vtk_sequence_writer = name_vtk_sequence_writer()

    predicate = name_predicate()

    return ["",
            "double T = {}.get<double>(\"instat.T\", 1.0);".format(ini),
            "double dt = {}.get<double>(\"instat.dt\", 0.1);".format(ini),
            "int step_number(0);"
            "int output_every_nth = {}.get<int>(\"instat.output_every_nth\", 1);".format(ini),
            "while (time<T-1e-8){",
            "{}".format(assemble_new_constraints),
            "  // Do time step",
            "  {} {}new({});".format(vector_type, vector, vector),
            "  {}".format(apply_call),
            "",
            "  // Accept new time step",
            "  {} = {}new;".format(vector, vector),
            "  time += dt;",
            "",
            "  step_number += 1;",
            "  if (step_number%output_every_nth == 0){",
            "    // Output to VTK File",
            "    {}.vtkWriter()->clear();".format(vtk_sequence_writer),
            "    Dune::PDELab::addSolutionToVTKWriter(vtkSequenceWriter, {}, {},".format(gfs, vector),
            "                                         Dune::PDELab::vtk::defaultNameScheme(), {});".format(predicate),
            "    {}.write({}, Dune::VTK::appendedraw);".format(vtk_sequence_writer, time),
            "  }",
            "}",
            ""]


@preamble(section="init")
def define_time(name):
    return "double {} = 0.0;".format(name)


def name_time():
    define_time("time")
    return "time"


@class_member(classtag="driver_block")
def typedef_timesteppingmethod(name):
    r_type = type_range()
    explicit = get_option('explicit_time_stepping')
    order = get_option('time_stepping_order')
    if explicit:
        if order == 1:
            return "using {} = Dune::PDELab::ExplicitEulerParameter<{}>;".format(name, r_type)
        elif order == 2:
            return "using {} = Dune::PDELab::HeunParameter<{}>;".format(name, r_type)
        elif order == 3:
            return "using {} = Dune::PDELab::Shu3Parameter<{}>;".format(name, r_type)
        elif order == 4:
            return "using {} = Dune::PDELab::RK4Parameter<{}>;".format(name, r_type)
        else:
            raise NotImplementedError("Time stepping method not supported")
    else:
        if order == 1:
            return "using {} = Dune::PDELab::OneStepThetaParameter<{}>;".format(name, r_type)
        elif order == 2:
            return "using {} = Dune::PDELab::Alexander2Parameter<{}>;".format(name, r_type)
        elif order == 3:
            return "using {} = Dune::PDELab::Alexander3Parameter<{}>;".format(name, r_type)


def type_timesteppingmethod():
    typedef_timesteppingmethod("TSM")
    return "TSM"


@preamble(section="instat")
def define_timesteppingmethod(name):
    tsm_type = type_timesteppingmethod()
    explicit = get_option('explicit_time_stepping')
    if explicit:
        return "{} {};".format(tsm_type, name)
    else:
        order = get_option('time_stepping_order')
        if order == 1:
            ini = name_initree()
            return "{} {}({}.get<double>(\"instat.theta\",1.0));".format(tsm_type, name, ini)
        else:
            return "{} {};".format(tsm_type, name)


def name_timesteppingmethod():
    define_timesteppingmethod("tsm")
    return "tsm"


@class_member(classtag="driver_block")
def typedef_instationarygridoperator(name):
    include_file("dune/pdelab/gridoperator/onestep.hh", filetag="driver")
    go_type = type_gridoperator(get_form_ident())
    mass_go_type = type_gridoperator("mass")
    explicit = get_option('explicit_time_stepping')
    if explicit:
        return "using {} = Dune::PDELab::OneStepGridOperator<{},{},false>;".format(name, go_type, mass_go_type)
    else:
        return "using {} = Dune::PDELab::OneStepGridOperator<{},{}>;".format(name, go_type, mass_go_type)


def type_instationarygridoperator():
    typedef_instationarygridoperator("IGO")
    return "IGO"


@preamble(section="gridoperator")
def define_instationarygridoperator(name):
    igo_type = type_instationarygridoperator()
    go = name_gridoperator(get_form_ident())
    mass_go = name_gridoperator("mass")
    return "{} {}({}, {});".format(igo_type, name, go, mass_go)


def name_instationarygridoperator():
    define_instationarygridoperator("igo")
    return "igo"


@class_member(classtag="driver_block")
def typedef_onestepmethod(name):
    r_type = type_range()
    igo_type = type_instationarygridoperator()
    snp_type = type_stationarynonlinearproblemssolver(igo_type)
    vector_type = type_vector(get_form_ident())
    return "using {} = Dune::PDELab::OneStepMethod<{}, {}, {}, {}, {}>;".format(name, r_type, igo_type, snp_type, vector_type, vector_type)


def type_onestepmethod():
    typedef_onestepmethod("OSM")
    return "OSM"


@preamble(section="instat")
def define_onestepmethod(name):
    ilptype = type_onestepmethod()
    tsm = name_timesteppingmethod()
    igo_type = type_instationarygridoperator()
    igo = name_instationarygridoperator()
    snp = name_stationarynonlinearproblemsolver(igo_type, igo)
    return "{} {}({},{},{});".format(ilptype, name, tsm, igo, snp)


def name_onestepmethod():
    define_onestepmethod("osm")
    return "osm"


@class_member(classtag="driver_block")
def typedef_explicitonestepmethod(name):
    r_type = type_range()
    igo_type = type_instationarygridoperator()
    ls_type = type_linearsolver()
    vector_type = type_vector(get_form_ident())
    return "using {} = Dune::PDELab::ExplicitOneStepMethod<{}, {}, {}, {}>;".format(name, r_type, igo_type, ls_type, vector_type)


def type_explicitonestepmethod():
    typedef_explicitonestepmethod("EOSM")
    return "EOSM"


@preamble(section="instat")
def define_explicitonestepmethod(name):
    eosm_type = type_explicitonestepmethod()
    tsm = name_timesteppingmethod()
    igo = name_instationarygridoperator()
    ls = name_linearsolver()
    return "{} {}({}, {}, {});".format(eosm_type, name, tsm, igo, ls)


def name_explicitonestepmethod():
    define_explicitonestepmethod("eosm")
    return "eosm"