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_function,
name_constraintscontainer,
)
from dune.codegen.pdelab.driver.interpolate import (interpolate_dirichlet_data,
name_boundary_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_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_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"
@preamble(section="instat")
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"
@preamble(section="gridoperator")
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"
@preamble(section="instat")
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"
@preamble(section="instat")
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"