Python sympy.Abs() Examples

def eliminate(rxns, wrt):
        """ Linear combination coefficients for elimination of a substance

        Parameters
        ----------
        rxns : iterable of Equilibrium instances
        wrt : str (substance key)

        Examples
        --------
        >>> e1 = Equilibrium({'Cd+2': 4, 'H2O': 4}, {'Cd4(OH)4+4': 1, 'H+': 4}, 10**-32.5)
        >>> e2 = Equilibrium({'Cd(OH)2(s)': 1}, {'Cd+2': 1, 'OH-': 2}, 10**-14.4)
        >>> Equilibrium.eliminate([e1, e2], 'Cd+2')
        [1, 4]
        >>> print(1*e1 + 4*e2)
        4 Cd(OH)2(s) + 4 H2O = Cd4(OH)4+4 + 4 H+ + 8 OH-; 7.94e-91

        """
        import sympy
        viol = [r.net_stoich([wrt])[0] for r in rxns]
        factors = defaultdict(int)
        for v in viol:
            for f in sympy.primefactors(v):
                factors[f] = max(factors[f], sympy.Abs(v//f))
        rcd = reduce(mul, (k**v for k, v in factors.items()))
        viol[0] *= -1
        return [rcd//v for v in viol] 

def test_events():
    # use bouncing ball to test events work

    # simulate in block diagram
    int_opts = block_diagram.DEFAULT_INTEGRATOR_OPTIONS.copy()
    int_opts['rtol'] = 1E-12
    int_opts['atol'] = 1E-15
    int_opts['nsteps'] = 1000
    int_opts['max_step'] = 2**-3
    x = x1, x2 = Array(dynamicsymbols('x_1:3'))
    mu, g = sp.symbols('mu g')
    constants = {mu: 0.8, g: 9.81}
    ic = np.r_[10, 15]
    sys = SwitchedSystem(
        x1, Array([0]),
        state_equations=r_[x2, -g],
        state_update_equation=r_[sp.Abs(x1), -mu*x2],
        state=x,
        constants_values=constants,
        initial_condition=ic
    )
    bd = BlockDiagram(sys)
    res = bd.simulate(5, integrator_options=int_opts)

    # compute actual impact time
    tvar = dynamicsymbols._t
    impact_eq = (x2*tvar - g*tvar**2/2 + x1).subs(
        {x1: ic[0], x2: ic[1], g: 9.81}
    )
    t_impact = sp.solve(impact_eq, tvar)[-1]

    # make sure simulation actually changes velocity sign around impact
    abs_diff_impact = np.abs(res.t - t_impact)
    impact_idx = np.where(abs_diff_impact == np.min(abs_diff_impact))[0]
    assert np.sign(res.x[impact_idx-1, 1]) != np.sign(res.x[impact_idx+1, 1]) 

def monkey_patch_function(expr):
    """
    Name says it all
    Monkey patch any functions with bad func names (abs -> Abs)
    If anything else needs to be fixed, goes here
    """
    return expr.replace(sympify("abs(x)").func, Abs) 

def initialize_damp(damp, nbl, spacing, abc_type="damp", fs=False):
    """
    Initialize damping field with an absorbing boundary layer.

    Parameters
    ----------
    damp : Function
        The damping field for absorbing boundary condition.
    nbl : int
        Number of points in the damping layer.
    spacing :
        Grid spacing coefficient.
    mask : bool, optional
        whether the dampening is a mask or layer.
        mask => 1 inside the domain and decreases in the layer
        not mask => 0 inside the domain and increase in the layer
    """
    dampcoeff = 1.5 * np.log(1.0 / 0.001) / (nbl)

    eqs = [Eq(damp, 1.0)] if abc_type == "mask" else []
    for d in damp.dimensions:
        if not fs or d is not damp.dimensions[-1]:
            # left
            dim_l = SubDimension.left(name='abc_%s_l' % d.name, parent=d,
                                      thickness=nbl)
            pos = Abs((nbl - (dim_l - d.symbolic_min) + 1) / float(nbl))
            val = dampcoeff * (pos - sin(2*np.pi*pos)/(2*np.pi))
            val = -val if abc_type == "mask" else val
            eqs += [Inc(damp.subs({d: dim_l}), val/d.spacing)]
        # right
        dim_r = SubDimension.right(name='abc_%s_r' % d.name, parent=d,
                                   thickness=nbl)
        pos = Abs((nbl - (d.symbolic_max - dim_r) + 1) / float(nbl))
        val = dampcoeff * (pos - sin(2*np.pi*pos)/(2*np.pi))
        val = -val if abc_type == "mask" else val
        eqs += [Inc(damp.subs({d: dim_r}), val/d.spacing)]

    Operator(eqs, name='initdamp')() 

def norm(f, order=2):
    """
    Compute the norm of a Function.

    Parameters
    ----------
    f : Function
        Input Function.
    order : int, optional
        The order of the norm. Defaults to 2.
    """
    kwargs = {}
    if f.is_TimeFunction and f._time_buffering:
        kwargs[f.time_dim.max_name] = f._time_size - 1

    # Protect SparseFunctions from accessing duplicated (out-of-domain) data,
    # otherwise we would eventually be summing more than expected
    p, eqns = f.guard() if f.is_SparseFunction else (f, [])

    s = dv.types.Scalar(name='sum', dtype=f.dtype)

    with MPIReduction(f) as mr:
        op = dv.Operator([dv.Eq(s, 0.0)] +
                         eqns +
                         [dv.Inc(s, Abs(Pow(p, order))), dv.Eq(mr.n[0], s)],
                         name='norm%d' % order)
        op.apply(**kwargs)

    v = Pow(mr.v, 1/order)

    return f.dtype(v) 

def convert_comp(comp):
    if comp.group():
        return convert_expr(comp.group().expr())
    elif comp.abs_group():
        return sympy.Abs(convert_expr(comp.abs_group().expr()), evaluate=False)
    elif comp.atom():
        return convert_atom(comp.atom())
    elif comp.frac():
        return convert_frac(comp.frac())
    elif comp.func():
        return convert_func(comp.func()) 

def setup_w_over_q(wOverQ, w, qmin, qmax, npad, sigma=0):
    """
    Initialise spatially variable w/Q field used to implement attenuation and
    absorb outgoing waves at the edges of the model. Uses Devito Operator.

    Parameters
    ----------
    wOverQ : Function, required
        The omega over Q field used to implement attenuation in the model,
        and the absorbing boundary condition for outgoing waves.
    w : float32, required
        center angular frequency, e.g. peak frequency of Ricker source wavelet
        used for modeling.
    qmin : float32, required
        Q value at the edge of the model. Typically set to 0.1 to strongly
        attenuate outgoing waves.
    qmax : float32, required
        Q value in the interior of the model. Typically set to 100 as a
        reasonable and physically meaningful Q value.
    npad : int, required
        Number of points in the absorbing boundary region. Note that we expect
        this to be the same on all sides of the model.
    sigma : float32, optional, defaults to None
        sigma value for call to scipy gaussian smoother, default 5.
    """
    # sanity checks
    assert w > 0, "supplied w value [%f] must be positive" % (w)
    assert qmin > 0, "supplied qmin value [%f] must be positive" % (qmin)
    assert qmax > 0, "supplied qmax value [%f] must be positive" % (qmax)
    assert npad > 0, "supplied npad value [%f] must be positive" % (npad)
    for n in wOverQ.grid.shape:
        if n - 2*npad < 1:
            raise ValueError("2 * npad must not exceed dimension size!")

    lqmin = np.log(qmin)
    lqmax = np.log(qmax)

    # 1. Get distance to closest boundary in all dimensions
    # 2. Logarithmic variation between qmin, qmax across the absorbing boundary
    eqs = [Eq(wOverQ, 1)]
    for d in wOverQ.dimensions:
        # left
        dim_l = SubDimension.left(name='abc_%s_l' % d.name, parent=d,
                                  thickness=npad)
        pos = Abs(dim_l - d.symbolic_min) / float(npad)
        eqs.append(Eq(wOverQ.subs({d: dim_l}), Min(wOverQ.subs({d: dim_l}), pos)))
        # right
        dim_r = SubDimension.right(name='abc_%s_r' % d.name, parent=d,
                                   thickness=npad)
        pos = Abs(d.symbolic_max - dim_r) / float(npad)
        eqs.append(Eq(wOverQ.subs({d: dim_r}), Min(wOverQ.subs({d: dim_r}), pos)))

    eqs.append(Eq(wOverQ, w / exp(lqmin + wOverQ * (lqmax - lqmin))))
    # 2020.05.04 currently does not support spatial smoothing of the Q field
    # due to MPI weirdness in reassignment of the numpy array
    Operator(eqs, name='WOverQ_Operator')()