SfePy NTC

Source code for sfepy.discrete.problem

import os
import os.path as op
import time
from copy import copy

import numpy as nm

from sfepy.base.base import dict_from_keys_init, select_by_names
from sfepy.base.base import output, get_default, Struct, IndexedStruct
import sfepy.base.ioutils as io
from sfepy.base.conf import ProblemConf, get_standard_keywords
from sfepy.base.conf import transform_variables, transform_materials
from functions import Functions
from sfepy.discrete.fem.mesh import Mesh
from sfepy.discrete.common.fields import fields_from_conf
from variables import Variables, Variable
from materials import Materials, Material
from equations import Equations
from integrals import Integrals
from sfepy.discrete.state import State
from sfepy.discrete.conditions import Conditions
from sfepy.discrete.evaluate import create_evaluable, eval_equations
from sfepy.discrete.fem import fea
from sfepy.solvers.ts import TimeStepper
from sfepy.discrete.evaluate import BasicEvaluator, LCBCEvaluator
from sfepy.solvers import Solver
from sfepy.solvers.ls import ScipyDirect
from sfepy.solvers.nls import Newton

##
# 29.01.2006, c
[docs]class Problem(Struct): """ Problem definition, the top-level class holding all data necessary to solve a problem. It can be constructed from a :class:`ProblemConf` instance using `Problem.from_conf()` or directly from a problem description file using `Problem.from_conf_file()` For interactive use, the constructor requires only the `equations`, `nls` and `ls` keyword arguments. """ @staticmethod
[docs] def from_conf_file(conf_filename, required=None, other=None, init_fields=True, init_equations=True, init_solvers=True): _required, _other = get_standard_keywords() if required is None: required = _required if other is None: other = _other conf = ProblemConf.from_file(conf_filename, required, other) obj = Problem.from_conf(conf, init_fields=init_fields, init_equations=init_equations, init_solvers=init_solvers) return obj
@staticmethod
[docs] def from_conf(conf, init_fields=True, init_equations=True, init_solvers=True): if conf.options.get('absolute_mesh_path', False): conf_dir = None else: conf_dir = op.dirname(conf.funmod.__file__) functions = Functions.from_conf(conf.functions) if conf.get('filename_mesh') is not None: from sfepy.discrete.fem.domain import FEDomain mesh = Mesh.from_file(conf.filename_mesh, prefix_dir=conf_dir) domain = FEDomain(mesh.name, mesh) if conf.options.get('ulf', False): domain.mesh.coors_act = domain.mesh.coors.copy() elif conf.get('filename_domain') is not None: from sfepy.discrete.iga.domain import IGDomain domain = IGDomain.from_file(conf.filename_domain) else: raise ValueError('missing filename_mesh or filename_domain!') obj = Problem('problem_from_conf', conf=conf, functions=functions, domain=domain, auto_conf=False, auto_solvers=False) obj.set_regions(conf.regions, obj.functions) obj.clear_equations() if init_fields: obj.set_fields(conf.fields) if init_equations: obj.set_equations(conf.equations, user={'ts' : obj.ts}) if init_solvers: obj.set_solvers(conf.solvers, conf.options) return obj
def __init__(self, name, conf=None, functions=None, domain=None, fields=None, equations=None, auto_conf=True, nls=None, ls=None, ts=None, auto_solvers=True): self.name = name self.conf = conf self.functions = functions self.reset() self.ts = get_default(ts, self.get_default_ts()) if auto_conf: if equations is None: raise ValueError('missing equations in auto_conf mode!') if fields is None: variables = equations.variables fields = {} for field in [var.get_field() for var in variables]: fields[field.name] = field if domain is None: domain = fields.values()[0].domain if conf is None: self.conf = Struct(ebcs={}, epbcs={}, lcbcs={}) self.equations = equations self.fields = fields self.domain = domain if auto_solvers: if ls is None: ls = ScipyDirect({}) if nls is None: nls = Newton({}, lin_solver=ls) ev = self.get_evaluator() nls.fun = ev.eval_residual nls.fun_grad = ev.eval_tangent_matrix self.set_solvers_instances(ls=ls, nls=nls) self.setup_output()
[docs] def reset(self): if hasattr(self.conf, 'options'): self.setup_hooks(self.conf.options) else: self.setup_hooks() self.mtx_a = None self.solvers = None self.clear_equations()
[docs] def setup_hooks(self, options=None): """ Setup various hooks (user-defined functions), as given in `options`. Supported hooks: - `matrix_hook` - check/modify tangent matrix in each nonlinear solver iteration - `nls_iter_hook` - called prior to every iteration of nonlinear solver, if the solver supports that - takes the Problem instance (`self`) as the first argument """ hook_names = ['nls_iter_hook', 'matrix_hook'] for hook_name in hook_names: setattr(self, hook_name, None) if options is not None: hook = options.get(hook_name, None) if hook is not None: hook = self.conf.get_function(hook) setattr(self, hook_name, hook) iter_hook = self.nls_iter_hook if iter_hook is not None: self.nls_iter_hook = lambda *args, **kwargs: \ iter_hook(self, *args, **kwargs)
[docs] def copy(self, name=None): """ Make a copy of Problem. """ if name is None: name = self.name + '_copy' obj = Problem(name, conf=self.conf, functions=self.functions, domain=self.domain, fields=self.fields, equations=self.equations, auto_conf=False, auto_solvers=False) obj.ebcs = self.ebcs obj.epbcs = self.epbcs obj.lcbcs = self.lcbcs obj.set_solvers(self.conf.solvers, self.conf.options) obj.setup_output(output_filename_trunk=self.ofn_trunk, output_dir=self.output_dir, output_format=self.output_format, file_per_var=self.file_per_var, linearization=self.linearization) return obj
[docs] def create_subproblem(self, var_names, known_var_names): """ Create a sub-problem with equations containing only terms with the given virtual variables. Parameters ---------- var_names : list The list of names of virtual variables. known_var_names : list The list of names of (already) known state variables. Returns ------- subpb : Problem instance The sub-problem. """ subpb = Problem(self.name + '_' + '_'.join(var_names), conf=self.conf, functions=self.functions, domain=self.domain, fields=self.fields, auto_conf=False, auto_solvers=False) subpb.set_solvers(self.conf.solvers, self.conf.options) subeqs = self.equations.create_subequations(var_names, known_var_names) subpb.set_equations_instance(subeqs, keep_solvers=True) return subpb
[docs] def setup_default_output(self, conf=None, options=None): """ Provide default values to `Problem.setup_output()` from `conf.options` and `options`. """ conf = get_default(conf, self.conf) if options and getattr(options, 'output_filename_trunk', None): default_output_dir, of = op.split(options.output_filename_trunk) default_trunk = io.get_trunk(of) else: default_trunk = None default_output_dir = conf.options.get('output_dir', None) if options and getattr(options, 'output_format', None): default_output_format = options.output_format else: default_output_format = conf.options.get('output_format', None) default_file_per_var = conf.options.get('file_per_var', None) default_float_format = conf.options.get('float_format', None) default_linearization = Struct(kind='strip') self.setup_output(output_filename_trunk=default_trunk, output_dir=default_output_dir, file_per_var=default_file_per_var, output_format=default_output_format, float_format=default_float_format, linearization=default_linearization)
[docs] def setup_output(self, output_filename_trunk=None, output_dir=None, output_format=None, float_format=None, file_per_var=None, linearization=None): """ Sets output options to given values, or uses the defaults for each argument that is None. """ self.output_modes = {'vtk' : 'sequence', 'h5' : 'single'} self.ofn_trunk = get_default(output_filename_trunk, io.get_trunk(self.domain.name)) self.set_output_dir(output_dir) self.output_format = get_default(output_format, 'vtk') self.float_format = get_default(float_format, None) self.file_per_var = get_default(file_per_var, False) self.linearization = get_default(linearization, Struct(kind='strip')) if ((self.output_format == 'h5') and (self.linearization.kind == 'adaptive')): self.linearization.kind = None
[docs] def set_output_dir(self, output_dir=None): """ Set the directory for output files. The directory is created if it does not exist. """ self.output_dir = get_default(output_dir, os.curdir) if self.output_dir and not op.exists(self.output_dir): os.makedirs(self.output_dir)
[docs] def set_regions(self, conf_regions=None, conf_materials=None, functions=None): conf_regions = get_default(conf_regions, self.conf.regions) functions = get_default(functions, self.functions) self.domain.create_regions(conf_regions, functions)
[docs] def set_materials(self, conf_materials=None): """ Set definition of materials. """ self.conf_materials = get_default(conf_materials, self.conf.materials)
[docs] def select_materials(self, material_names, only_conf=False): if type(material_names) == dict: conf_materials = transform_materials(material_names) else: conf_materials = select_by_names(self.conf.materials, material_names) if not only_conf: self.set_materials(conf_materials) return conf_materials
[docs] def set_fields(self, conf_fields=None): conf_fields = get_default(conf_fields, self.conf.fields) self.fields = fields_from_conf(conf_fields, self.domain.regions)
[docs] def set_variables(self, conf_variables=None): """ Set definition of variables. """ self.conf_variables = get_default(conf_variables, self.conf.variables) self.reset()
[docs] def select_variables(self, variable_names, only_conf=False): if type(variable_names) == dict: conf_variables = transform_variables(variable_names) else: conf_variables = select_by_names(self.conf.variables, variable_names) if not only_conf: self.set_variables(conf_variables) return conf_variables
[docs] def clear_equations(self): self.integrals = None self.equations = None self.ebcs = None self.epbcs = None self.lcbcs = None
[docs] def set_equations(self, conf_equations=None, user=None, keep_solvers=False, make_virtual=False): """ Set equations of the problem using the `equations` problem description entry. Fields and Regions have to be already set. """ conf_equations = get_default(conf_equations, self.conf.get('equations', None)) self.set_variables() variables = Variables.from_conf(self.conf_variables, self.fields) self.set_materials() materials = Materials.from_conf(self.conf_materials, self.functions) self.integrals = self.get_integrals() equations = Equations.from_conf(conf_equations, variables, self.domain.regions, materials, self.integrals, user=user) self.equations = equations if not keep_solvers: self.solvers = None
[docs] def set_equations_instance(self, equations, keep_solvers=False): """ Set equations of the problem to `equations`. """ self.mtx_a = None self.clear_equations() self.equations = equations if not keep_solvers: self.solvers = None
[docs] def set_solvers(self, conf_solvers=None, options=None): """ Choose which solvers should be used. If solvers are not set in `options`, use first suitable in `conf_solvers`. """ conf_solvers = get_default(conf_solvers, self.conf.solvers) self.solver_confs = {} for key, val in conf_solvers.iteritems(): self.solver_confs[val.name] = val def _find_suitable(prefix): for key, val in self.solver_confs.iteritems(): if val.kind.find(prefix) == 0: return val return None def _get_solver_conf(kind): try: key = options[kind] conf = self.solver_confs[key] except: conf = _find_suitable(kind + '.') return conf self.ts_conf = _get_solver_conf('ts') if self.ts_conf is None: self.ts_conf = Struct(name='no ts', kind='ts.stationary') self.nls_conf = _get_solver_conf('nls') self.ls_conf = _get_solver_conf('ls') info = 'using solvers:' if self.ts_conf: info += '\n ts: %s' % self.ts_conf.name if self.nls_conf: info += '\n nls: %s' % self.nls_conf.name if self.ls_conf: info += '\n ls: %s' % self.ls_conf.name if info != 'using solvers:': output(info)
[docs] def set_solvers_instances(self, ls=None, nls=None): """ Set the instances of linear and nonlinear solvers that will be used in `Problem.solve()` call. """ if (ls is not None) and (nls is not None): if not (nls.lin_solver is ls): raise ValueError('linear solver not used in nonlinear!') self.solvers = Struct(name='solvers', ls=ls, nls=nls) if nls is not None: self.nls_status = get_default(nls.status, IndexedStruct())
[docs] def get_solver_conf(self, name): return self.solver_confs[name]
[docs] def get_default_ts(self, t0=None, t1=None, dt=None, n_step=None, step=None): t0 = get_default(t0, 0.0) t1 = get_default(t1, 1.0) dt = get_default(dt, 1.0) n_step = get_default(n_step, 1) ts = TimeStepper(t0, t1, dt, n_step, step=step) return ts
[docs] def get_integrals(self, names=None): """ Get integrals, initialized from problem configuration if available. Parameters ---------- names : list, optional If given, only the named integrals are returned. Returns ------- integrals : Integrals instance The requested integrals. """ conf_integrals = self.conf.get('integrals', {}) integrals = Integrals.from_conf(conf_integrals) if names is not None: integrals.update([integrals[ii] for ii in names if ii in integrals.names]) return integrals
[docs] def update_time_stepper(self, ts): if ts is not None: self.ts = ts
[docs] def update_materials(self, ts=None, mode='normal', verbose=True): """ Update materials used in equations. Parameters ---------- ts : TimeStepper instance The time stepper. mode : 'normal', 'update' or 'force' The update mode, see :func:`Material.time_update() <sfepy.discrete.materials.Material.time_update()>`. verbose : bool If False, reduce verbosity. """ if self.equations is not None: self.update_time_stepper(ts) self.equations.time_update_materials(self.ts, mode=mode, problem=self, verbose=verbose)
[docs] def update_equations(self, ts=None, ebcs=None, epbcs=None, lcbcs=None, functions=None, create_matrix=False): """ Update equations for current time step. The tangent matrix graph is automatically recomputed if the set of active essential or periodic boundary conditions changed w.r.t. the previous time step. Parameters ---------- ts : TimeStepper instance, optional The time stepper. If not given, `self.ts` is used. ebcs : Conditions instance, optional The essential (Dirichlet) boundary conditions. If not given, `self.ebcs` are used. epbcs : Conditions instance, optional The periodic boundary conditions. If not given, `self.epbcs` are used. lcbcs : Conditions instance, optional The linear combination boundary conditions. If not given, `self.lcbcs` are used. functions : Functions instance, optional The user functions for boundary conditions, materials, etc. If not given, `self.functions` are used. """ self.update_time_stepper(ts) functions = get_default(functions, self.functions) graph_changed = self.equations.time_update(self.ts, ebcs, epbcs, lcbcs, functions, self) self.graph_changed = graph_changed if graph_changed or (self.mtx_a is None) or create_matrix: self.mtx_a = self.equations.create_matrix_graph() ## import sfepy.base.plotutils as plu ## plu.spy(self.mtx_a) ## plu.plt.show()
[docs] def set_bcs(self, ebcs=None, epbcs=None, lcbcs=None): """ Update boundary conditions. """ if isinstance(ebcs, Conditions): self.ebcs = ebcs else: conf_ebc = get_default(ebcs, self.conf.ebcs) self.ebcs = Conditions.from_conf(conf_ebc, self.domain.regions) if isinstance(epbcs, Conditions): self.epbcs = epbcs else: conf_epbc = get_default(epbcs, self.conf.epbcs) self.epbcs = Conditions.from_conf(conf_epbc, self.domain.regions) if isinstance(lcbcs, Conditions): self.lcbcs = lcbcs else: conf_lcbc = get_default(lcbcs, self.conf.lcbcs) self.lcbcs = Conditions.from_conf(conf_lcbc, self.domain.regions)
[docs] def time_update(self, ts=None, ebcs=None, epbcs=None, lcbcs=None, functions=None, create_matrix=False): self.set_bcs(ebcs, epbcs, lcbcs) self.update_equations(ts, self.ebcs, self.epbcs, self.lcbcs, functions, create_matrix)
[docs] def setup_ic(self, conf_ics=None, functions=None): conf_ics = get_default(conf_ics, self.conf.ics) ics = Conditions.from_conf(conf_ics, self.domain.regions) functions = get_default(functions, self.functions) self.equations.setup_initial_conditions(ics, functions)
[docs] def select_bcs(self, ebc_names=None, epbc_names=None, lcbc_names=None, create_matrix=False): if ebc_names is not None: conf_ebc = select_by_names(self.conf.ebcs, ebc_names) else: conf_ebc = None if epbc_names is not None: conf_epbc = select_by_names(self.conf.epbcs, epbc_names) else: conf_epbc = None if lcbc_names is not None: conf_lcbc = select_by_names(self.conf.lcbcs, lcbc_names) else: conf_lcbc = None self.set_bcs(conf_ebc, conf_epbc, conf_lcbc) self.update_equations(self.ts, self.ebcs, self.epbcs, self.lcbcs, self.functions, create_matrix)
[docs] def get_timestepper(self): return self.ts
[docs] def create_state(self): return State(self.equations.variables)
[docs] def get_mesh_coors(self): return self.domain.get_mesh_coors()
[docs] def set_mesh_coors(self, coors, update_fields=False, actual=False, clear_all=True): """ Set mesh coordinates. Parameters ---------- coors : array The new coordinates. update_fields : bool If True, update also coordinates of fields. actual : bool If True, update the actual configuration coordinates, otherwise the undeformed configuration ones. """ fea.set_mesh_coors(self.domain, self.fields, coors, update_fields=update_fields, actual=actual, clear_all=clear_all)
[docs] def refine_uniformly(self, level): """ Refine the mesh uniformly `level`-times. Notes ----- This operation resets almost everything (fields, equations, ...) - it is roughly equivalent to creating a new Problem instance with the refined mesh. """ if level == 0: return domain = self.domain for ii in range(level): domain = domain.refine() self.domain = domain self.set_regions(self.conf.regions, self.functions) self.clear_equations() self.set_fields(self.conf.fields) self.set_equations(self.conf.equations, user={'ts' : self.ts})
[docs] def get_dim(self, get_sym=False): """Returns mesh dimension, symmetric tensor dimension (if `get_sym` is True). """ dim = self.domain.mesh.dim if get_sym: return dim, (dim + 1) * dim / 2 else: return dim
[docs] def init_time(self, ts): self.update_time_stepper(ts) self.equations.init_time(ts) self.update_materials(mode='force')
[docs] def advance(self, ts=None): self.update_time_stepper(ts) self.equations.advance(self.ts)
[docs] def save_state(self, filename, state=None, out=None, fill_value=None, post_process_hook=None, linearization=None, file_per_var=False, **kwargs): """ Parameters ---------- file_per_var : bool or None If True, data of each variable are stored in a separate file. If None, it is set to the application option value. linearization : Struct or None The linearization configuration for higher order approximations. If its kind is 'adaptive', `file_per_var` is assumed True. """ linearization = get_default(linearization, self.linearization) if linearization.kind != 'adaptive': file_per_var = get_default(file_per_var, self.file_per_var) else: file_per_var = True extend = not file_per_var if (out is None) and (state is not None): out = state.create_output_dict(fill_value=fill_value, extend=extend, linearization=linearization) if post_process_hook is not None: out = post_process_hook(out, self, state, extend=extend) if linearization.kind == 'adaptive': for key, val in out.iteritems(): mesh = val.get('mesh', self.domain.mesh) aux = io.edit_filename(filename, suffix='_' + val.var_name) mesh.write(aux, io='auto', out={key : val}, float_format=self.float_format, **kwargs) if hasattr(val, 'levels'): output('max. refinement per group:', val.levels) elif file_per_var: meshes = {} if self.equations is None: varnames = {} for key, val in out.iteritems(): varnames[val.var_name] = 1 varnames = varnames.keys() outvars = self.create_variables(varnames) itervars = outvars.__iter__ else: itervars = self.equations.variables.iter_state for var in itervars(): rname = var.field.region.name if meshes.has_key(rname): mesh = meshes[rname] else: mesh = Mesh.from_region(var.field.region, self.domain.mesh, localize=True, is_surface=var.is_surface) meshes[rname] = mesh vout = {} for key, val in out.iteritems(): try: if val.var_name == var.name: vout[key] = val except AttributeError: msg = 'missing var_name attribute in output!' raise ValueError(msg) aux = io.edit_filename(filename, suffix='_' + var.name) mesh.write(aux, io='auto', out=vout, float_format=self.float_format, **kwargs) else: mesh = out.pop('__mesh__', self.domain.mesh) mesh.write(filename, io='auto', out=out, float_format=self.float_format, **kwargs)
[docs] def save_ebc(self, filename, ebcs=None, epbcs=None, force=True, default=0.0): """ Save essential boundary conditions as state variables. Parameters ---------- filename : str The output file name. ebcs : Conditions instance, optional The essential (Dirichlet) boundary conditions. If not given, `self.conf.ebcs` are used. epbcs : Conditions instance, optional The periodic boundary conditions. If not given, `self.conf.epbcs` are used. force : bool If True, sequential nonzero values are forced to individual `ebcs` so that the conditions are visible even when zero. default : float The default constant value of state vector. """ output('saving ebc...') variables = self.get_variables(auto_create=True) if ebcs is None: ebcs = Conditions.from_conf(self.conf.ebcs, self.domain.regions) if epbcs is None: epbcs = Conditions.from_conf(self.conf.epbcs, self.domain.regions) try: variables.equation_mapping(ebcs, epbcs, self.ts, self.functions, problem=self) except: output('cannot make equation mapping!') raise state = State(variables) state.fill(default) if force: vals = dict_from_keys_init(variables.state) for ii, key in enumerate(vals.iterkeys()): vals[key] = ii + 1 state.apply_ebc(force_values=vals) else: state.apply_ebc() out = state.create_output_dict(extend=True) self.save_state(filename, out=out, fill_value=default) output('...done')
[docs] def save_regions(self, filename_trunk, region_names=None): """ Save regions as meshes. Parameters ---------- filename_trunk : str The output filename without suffix. region_names : list, optional If given, only the listed regions are saved. """ filename = '%s.mesh' % filename_trunk self.domain.save_regions(filename, region_names=region_names)
[docs] def save_regions_as_groups(self, filename_trunk, region_names=None): """ Save regions in a single mesh but mark them by using different element/node group numbers. See :func:`Domain.save_regions_as_groups() <sfepy.discrete.fem.domain.Domain.save_regions_as_groups()>` for more details. Parameters ---------- filename_trunk : str The output filename without suffix. region_names : list, optional If given, only the listed regions are saved. """ filename = '%s.%s' % (filename_trunk, self.output_format) self.domain.save_regions_as_groups(filename, region_names=region_names)
[docs] def save_field_meshes(self, filename_trunk): output('saving field meshes...') for field in self.fields: output(field.name) field.write_mesh(filename_trunk + '_%s') output('...done')
[docs] def get_evaluator(self, reuse=False): """ Either create a new Evaluator instance (reuse == False), or return an existing instance, created in a preceding call to Problem.init_solvers(). """ if reuse: try: ev = self.evaluator except AttributeError: raise AttributeError('call Problem.init_solvers() or'\ ' set reuse to False!') else: if self.equations.variables.has_lcbc: ev = LCBCEvaluator(self, matrix_hook=self.matrix_hook) else: ev = BasicEvaluator(self, matrix_hook=self.matrix_hook) self.evaluator = ev return ev
[docs] def init_solvers(self, nls_status=None, ls_conf=None, nls_conf=None, mtx=None, presolve=False): """Create and initialize solvers.""" ls_conf = get_default(ls_conf, self.ls_conf, 'you must set linear solver!') nls_conf = get_default(nls_conf, self.nls_conf, 'you must set nonlinear solver!') if presolve: tt = time.clock() ls = Solver.any_from_conf(ls_conf, mtx=mtx, presolve=presolve, problem=self) if presolve: tt = time.clock() - tt output('presolve: %.2f [s]' % tt) ev = self.get_evaluator() if self.conf.options.get('ulf', False): self.nls_iter_hook = ev.new_ulf_iteration nls = Solver.any_from_conf(nls_conf, fun=ev.eval_residual, fun_grad=ev.eval_tangent_matrix, lin_solver=ls, iter_hook=self.nls_iter_hook, status=nls_status, problem=self) self.set_solvers_instances(ls=ls, nls=nls)
[docs] def get_solvers(self): return getattr(self, 'solvers', None)
[docs] def is_linear(self): nls_conf = get_default(None, self.nls_conf, 'you must set nonlinear solver!') aux = Solver.any_from_conf(nls_conf) if aux.conf.get('is_linear', False): return True else: return False
[docs] def set_linear(self, is_linear): nls_conf = get_default(None, self.nls_conf, 'you must set nonlinear solver!') nls_conf.is_linear = is_linear
[docs] def solve(self, state0=None, nls_status=None, ls_conf=None, nls_conf=None, force_values=None, var_data=None): """Solve self.equations in current time step. Parameters ---------- var_data : dict A dictionary of {variable_name : data vector} used to initialize parameter variables. """ solvers = self.get_solvers() if solvers is None: self.init_solvers(nls_status, ls_conf, nls_conf) solvers = self.get_solvers() if state0 is None: state0 = State(self.equations.variables) else: if isinstance(state0, nm.ndarray): state0 = State(self.equations.variables, vec=state0) self.equations.set_data(var_data, ignore_unknown=True) self.update_materials() state0.apply_ebc(force_values=force_values) vec0 = state0.get_reduced() self.nls_status = get_default(nls_status, self.nls_status) vec = solvers.nls(vec0, status=self.nls_status) state = state0.copy(preserve_caches=True) state.set_reduced(vec, preserve_caches=True) return state
[docs] def create_evaluable(self, expression, try_equations=True, auto_init=False, preserve_caches=False, copy_materials=True, integrals=None, ebcs=None, epbcs=None, lcbcs=None, ts=None, functions=None, mode='eval', var_dict=None, strip_variables=True, extra_args=None, verbose=True, **kwargs): """ Create evaluable object (equations and corresponding variables) from the `expression` string. Convenience function calling :func:`create_evaluable() <sfepy.discrete.evaluate.create_evaluable()>` with defaults provided by the Problem instance `self`. The evaluable can be repeatedly evaluated by calling :func:`eval_equations() <sfepy.discrete.evaluate.eval_equations()>`, e.g. for different values of variables. Parameters ---------- expression : str The expression to evaluate. try_equations : bool Try to get variables from `self.equations`. If this fails, variables can either be provided in `var_dict`, as keyword arguments, or are created automatically according to the expression. auto_init : bool Set values of all variables to all zeros. preserve_caches : bool If True, do not invalidate evaluate caches of variables. copy_materials : bool Work with a copy of `self.equations.materials` instead of reusing them. Safe but can be slow. integrals : Integrals instance, optional The integrals to be used. Automatically created as needed if not given. ebcs : Conditions instance, optional The essential (Dirichlet) boundary conditions for 'weak' mode. If not given, `self.ebcs` are used. epbcs : Conditions instance, optional The periodic boundary conditions for 'weak' mode. If not given, `self.epbcs` are used. lcbcs : Conditions instance, optional The linear combination boundary conditions for 'weak' mode. If not given, `self.lcbcs` are used. ts : TimeStepper instance, optional The time stepper. If not given, `self.ts` is used. functions : Functions instance, optional The user functions for boundary conditions, materials etc. If not given, `self.functions` are used. mode : one of 'eval', 'el_avg', 'qp', 'weak' The evaluation mode - 'weak' means the finite element assembling, 'qp' requests the values in quadrature points, 'el_avg' element averages and 'eval' means integration over each term region. var_dict : dict, optional The variables (dictionary of (variable name) : (Variable instance)) to be used in the expression. Use this if the name of a variable conflicts with one of the parameters of this method. strip_variables : bool If False, the variables in `var_dict` or `kwargs` not present in the expression are added to the actual variables as a context. extra_args : dict, optional Extra arguments to be passed to terms in the expression. verbose : bool If False, reduce verbosity. **kwargs : keyword arguments Additional variables can be passed as keyword arguments, see `var_dict`. Returns ------- equations : Equations instance The equations that can be evaluated. variables : Variables instance The corresponding variables. Set their values and use :func:`eval_equations() <sfepy.discrete.evaluate.eval_equations()>`. Examples -------- `problem` is Problem instance. >>> out = problem.create_evaluable('dq_state_in_volume_qp.i1.Omega(u)') >>> equations, variables = out `vec` is a vector of coefficients compatible with the field of 'u' - let's use all ones. >>> vec = nm.ones((variables['u'].n_dof,), dtype=nm.float64) >>> variables['u'].set_data(vec) >>> vec_qp = eval_equations(equations, variables, mode='qp') Try another vector: >>> vec = 3 * nm.ones((variables['u'].n_dof,), dtype=nm.float64) >>> variables['u'].set_data(vec) >>> vec_qp = eval_equations(equations, variables, mode='qp') """ from sfepy.discrete.equations import get_expression_arg_names variables = get_default(var_dict, {}) var_context = get_default(var_dict, {}) if try_equations and self.equations is not None: # Make a copy, so that possible variable caches are preserved. for key, var in self.equations.variables.as_dict().iteritems(): if key in variables: continue var = var.copy(name=key) if not preserve_caches: var.clear_evaluate_cache() variables[key] = var elif var_dict is None: possible_var_names = get_expression_arg_names(expression) variables = self.create_variables(possible_var_names) materials = self.get_materials() if materials is not None: if copy_materials: materials = materials.semideep_copy() else: materials = Materials(objs=materials._objs) else: possible_mat_names = get_expression_arg_names(expression) materials = self.create_materials(possible_mat_names) _kwargs = copy(kwargs) for key, val in kwargs.iteritems(): if isinstance(val, Variable): if val.name != key: msg = 'inconsistent variable name! (%s == %s)' \ % (val.name, key) raise ValueError(msg) var_context[key] = variables[key] = val.copy(name=key) _kwargs.pop(key) elif isinstance(val, Material): if val.name != key: msg = 'inconsistent material name! (%s == %s)' \ % (val.name, key) raise ValueError(msg) materials[val.name] = val _kwargs.pop(key) kwargs = _kwargs ebcs = get_default(ebcs, self.ebcs) epbcs = get_default(epbcs, self.epbcs) lcbcs = get_default(lcbcs, self.lcbcs) ts = get_default(ts, self.get_timestepper()) functions = get_default(functions, self.functions) integrals = get_default(integrals, self.get_integrals()) out = create_evaluable(expression, self.fields, materials, variables.itervalues(), integrals, ebcs=ebcs, epbcs=epbcs, lcbcs=lcbcs, ts=ts, functions=functions, auto_init=auto_init, mode=mode, extra_args=extra_args, verbose=verbose, kwargs=kwargs) if not strip_variables: variables = out[1] variables.extend([var for var in var_context.itervalues() if var not in variables]) equations = out[0] mode = 'update' if not copy_materials else 'normal' equations.time_update_materials(self.ts, mode=mode, problem=self, verbose=verbose) return out
[docs] def evaluate(self, expression, try_equations=True, auto_init=False, preserve_caches=False, copy_materials=True, integrals=None, ebcs=None, epbcs=None, lcbcs=None, ts=None, functions=None, mode='eval', dw_mode='vector', term_mode=None, var_dict=None, strip_variables=True, ret_variables=False, verbose=True, extra_args=None, **kwargs): """ Evaluate an expression, convenience wrapper of :func:`Problem.create_evaluable` and :func:`eval_equations() <sfepy.discrete.evaluate.eval_equations>`. Parameters ---------- dw_mode : 'vector' or 'matrix' The assembling mode for 'weak' evaluation mode. term_mode : str The term call mode - some terms support different call modes and depending on the call mode different values are returned. ret_variables : bool If True, return the variables that were created to evaluate the expression. other : arguments See docstrings of :func:`Problem.create_evaluable()`. Returns ------- out : array The result of the evaluation. variables : Variables instance The variables that were created to evaluate the expression. Only provided if `ret_variables` is True. """ aux = self.create_evaluable(expression, try_equations=try_equations, auto_init=auto_init, preserve_caches=preserve_caches, copy_materials=copy_materials, integrals=integrals, ebcs=ebcs, epbcs=epbcs, lcbcs=lcbcs, ts=ts, functions=functions, mode=mode, var_dict=var_dict, strip_variables=strip_variables, extra_args=extra_args, verbose=verbose, **kwargs) equations, variables = aux out = eval_equations(equations, variables, preserve_caches=preserve_caches, mode=mode, dw_mode=dw_mode, term_mode=term_mode) if ret_variables: out = (out, variables) return out
[docs] def eval_equations(self, names=None, preserve_caches=False, mode='eval', dw_mode='vector', term_mode=None, verbose=True): """ Evaluate (some of) the problem's equations, convenience wrapper of :func:`eval_equations() <sfepy.discrete.evaluate.eval_equations>`. Parameters ---------- names : str or sequence of str, optional Evaluate only equations of the given name(s). preserve_caches : bool If True, do not invalidate evaluate caches of variables. mode : one of 'eval', 'el_avg', 'qp', 'weak' The evaluation mode - 'weak' means the finite element assembling, 'qp' requests the values in quadrature points, 'el_avg' element averages and 'eval' means integration over each term region. dw_mode : 'vector' or 'matrix' The assembling mode for 'weak' evaluation mode. term_mode : str The term call mode - some terms support different call modes and depending on the call mode different values are returned. verbose : bool If False, reduce verbosity. Returns ------- out : dict or result The evaluation result. In 'weak' mode it is the vector or sparse matrix, depending on `dw_mode`. Otherwise, it is a dict of results with equation names as keys or a single result for a single equation. """ return eval_equations(self.equations, self.equations.variables, names=names, preserve_caches=preserve_caches, mode=mode, dw_mode=dw_mode, term_mode=term_mode, verbose=verbose)
[docs] def get_time_solver(self, ts_conf=None, **kwargs): """ Create and return a TimeSteppingSolver instance. Notes ----- Also sets `self.ts` attribute. """ ts_conf = get_default(ts_conf, self.ts_conf, 'you must set time-stepping solver!') ts_solver = Solver.any_from_conf(ts_conf, problem=self, **kwargs) self.ts = ts_solver.ts return ts_solver
[docs] def get_materials(self): if self.equations is not None: materials = self.equations.materials else: materials = None return materials
[docs] def create_materials(self, mat_names=None): """ Create materials with names in `mat_names`. Their definitions have to be present in `self.conf.materials`. Notes ----- This method does not change `self.equations`, so it should not have any side effects. """ if mat_names is not None: conf_materials = self.select_materials(mat_names, only_conf=True) else: conf_materials = self.conf.materials materials = Materials.from_conf(conf_materials, self.functions) return materials
[docs] def init_variables(self, state): """Initialize variables with history.""" self.equations.variables.init_state(state)
[docs] def get_variables(self, auto_create=False): if self.equations is not None: variables = self.equations.variables elif auto_create: variables = self.create_variables() else: variables = None return variables
[docs] def create_variables(self, var_names=None): """ Create variables with names in `var_names`. Their definitions have to be present in `self.conf.variables`. Notes ----- This method does not change `self.equations`, so it should not have any side effects. """ if var_names is not None: conf_variables = self.select_variables(var_names, only_conf=True) else: conf_variables = self.conf.variables variables = Variables.from_conf(conf_variables, self.fields) return variables
[docs] def get_output_name(self, suffix=None, extra=None, mode=None): """ Return default output file name, based on the output directory, output format, step suffix and mode. If present, the extra string is put just before the output format suffix. """ out = op.join(self.output_dir, self.ofn_trunk) if suffix is not None: if mode is None: mode = self.output_modes[self.output_format] if mode == 'sequence': out = '.'.join((out, suffix)) if extra is not None: out = '.'.join((out, extra, self.output_format)) else: out = '.'.join((out, self.output_format)) return out
[docs] def remove_bcs(self): """ Convenience function to remove boundary conditions. """ self.time_update(ebcs={}, epbcs={}, lcbcs={})