Source code for simudo.physics.steppers
import logging
import numpy as np
from cached_property import cached_property
from ..fem import (
AdaptiveStepper, ConstantStepperMixin, NewtonSolver, NewtonSolverMaxDu)
from ..util import SetattrInitMixin
from ..io.output_writer import OutputWriter # FIXME
[docs]class NewtonBailoutException(Exception):
pass
[docs]class NonequilibriumCoupledStepper(AdaptiveStepper):
'''Adaptively solve a coupled poisson-drift diffusion-optical problem.
The optical problem can be solved self-consistently with the drift-diffusion problem.
Parameters
----------
solution: :py:class:`~.problem_data.ProblemData`
Solution on which to operate. This solution object will be
progressively re-solved at each of the target parameter values.
selfconsistent_optics: bool
If True, the optical problem will be re-solved
at every Newton iteration, allowing for problems where optical absorption
and carrier concentrations are interdependent.
output_writer: :py:class:`~.output_writer.OutputWriter`
If ``output_writer`` is a string, it is used as a filename to be
passed to the default OutputWriter object. If ``output_writer`` is
an object inherited from OutputWriter, this object will be used
instead. If None, then no output is written. Output will be
written after each target_value is solved.
'''
def __init__(self, **kwargs):
output_writer = kwargs.get('output_writer', None)
if output_writer is not None:
if isinstance(output_writer, str):
output_writer = OutputWriter(output_writer, plot_iv=True)
output_writer.stepper = self
kwargs['output_writer'] = output_writer
super().__init__(**kwargs)
@cached_property
def to_save_objects(self):
pdd = self.solution.pdd
r = pdd.get_to_save()
r.update(self.solution.optical.get_to_save())
return {k: v.value for k, v in r.items()
if v.solver_save}
[docs] def user_make_solver(self, solution):
sl = solution
pdd = solution.pdd
# previous_successful = len(self.get_last_successful(1)) > 0
# reset optical auxiliary functions
if self.selfconsistent_optics:
for o in solution.optical.fields:
o.update_output()
o.update_input()
# exercise caution when there's no solution to backtrack to
be_extra_careful = not self.get_last_successful(1)
def omega_cb(self):
if not be_extra_careful:
return 1.0
du = self.solution.du_norm
return 1.0
# PDD solver
solver = NewtonSolver.from_nice_obj(pdd)
solver.parameters.update(
maximum_iterations=200 if be_extra_careful else 25,
omega_cb=omega_cb,
extra_iterations=3 if be_extra_careful else 0,
minimum_iterations=3,
relative_tolerance=1e-4,
absolute_tolerance=1)
def _assign_scalar(subfunc, value):
fsr.assign_scalar(
subfunc.magnitude, value.m_as(subfunc.units))
fsr = pdd.mesh_util.function_subspace_registry
U = sl.unit_registry
space = pdd.mixed_function_helper.solution_mixed_space
# du_clip_mf = space.make_function()
# du_clip = du_clip_mf.function
# du_clip_split = du_clip_mf.split()
#
# du_clip.vector()[:] = np.inf
# FIXME: this should probably be elsewhere
# for k, v in du_clip_split.items():
# if k.endswith("/delta_w"):
# _assign_scalar(v, 0.005*U.eV)
# if k == 'poisson/phi':
# _assign_scalar(v, 0.005*U.volt)
# solver.parameters.update(
# maximum_du=du_clip.vector()[:],
# )
tol_func = pdd.mixed_function_helper.make_tolerance_function()
solver.parameters.update(
absolute_tolerance_vector=tol_func.function.vector()[:]
)
def run_optical_subsolvers(slv):
for o in solution.optical.fields:
o.update_input()
s = NewtonSolver.from_nice_obj(o)
# FIXME: these parameters are disgraceful
s.parameters.update(
maximum_iterations=2,
extra_iterations=2,
omega_cb=lambda self: 1.0)
tol_func = o.mixed_function_helper.make_tolerance_function()
s.parameters.update(
absolute_tolerance_vector=tol_func.function.vector()[:]
)
s.logger = logging.getLogger('newton.optical')
s.solve()
o.update_output()
def bailout_if_error_too_large(slv):
if (slv.solution.du_norm > 1e40) and (slv.iteration > 20):
raise NewtonBailoutException()
def rebase_w(slv):
for b in pdd.bands:
if hasattr(b, 'mixedqfl_do_rebase_w'):
b.mixedqfl_do_rebase_w()
# these optical hooks could be handled with code below.
if self.selfconsistent_optics:
solver.user_before_first_iteration_hooks.append(run_optical_subsolvers)
solver.user_post_iteration_hooks.append(run_optical_subsolvers)
def _append_not_None(lst, x):
if x is not None:
lst.append(x)
# register calculations to be done for each process.
for proc in self.solution.pdd.electro_optical_processes:
_append_not_None(solver.user_pre_iteration_hooks, proc.pre_iteration_hook)
_append_not_None(solver.user_before_first_iteration_hooks, proc.pre_first_iteration_hook)
_append_not_None(solver.user_post_iteration_hooks, proc.post_iteration_hook)
solver.user_before_first_iteration_hooks.append(bailout_if_error_too_large)
solver.user_post_iteration_hooks.append(bailout_if_error_too_large)
solver.user_pre_iteration_hooks.append(rebase_w)
return solver
[docs] def user_solver(self, solution, parameter):
logging.info("%%%%%%%%%% adaptive solve parameter={} step_size={}".format(parameter, self.step_size))
try:
return super().user_solver(solution, parameter)
except NewtonBailoutException:
return False
[docs]class NonequilibriumCoupledConstantStepper(
ConstantStepperMixin,
NonequilibriumCoupledStepper):
pass
[docs]class VoltageStepper(NonequilibriumCoupledConstantStepper):
'''Starting from an existing solution, which may be at thermal
equilibrium or with previously applied bias of illumination, ramp the
bias at one or more contacts.'''
update_parameter_success_factor = 1.5
update_parameter_failure_factor = 0.5
step_size = 0.1
[docs] def user_apply_parameter_to_solution(self, solution, parameter_value):
super().user_apply_parameter_to_solution(solution, parameter_value)
[docs]class OpticalIntensityAdaptiveStepper(NonequilibriumCoupledConstantStepper):
'''Adaptively increase the optical intensity by increasing
:py:attr:`.optical.Optical.Phi_scale`.
'''
update_parameter_success_factor = 3
update_parameter_failure_factor = 0.2
step_size = 1e-30
@cached_property
def constants(self):
return [self.solution.optical.Phi_scale]
@cached_property
def parameter_unit(self):
return self.unit_registry.dimensionless
@cached_property
def parameter_target_values(self):
return [1.0]