# Path: dolfin_navier_scipy/residual_checks.py
# def get_imex_res(V=None, outflowds=None, gradvsymmtrc=True, nu=None,
#                  implscheme='crni', explscheme='abtw'):
#     """ define the residual for an IMEX/AB2 time discretization
# 
#     """
#     if not implscheme == 'crni':
#         raise NotImplementedError()
# 
#     if explscheme == 'abtw':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (1.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx -
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'heun':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx +
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'eule':
#         def convform(cvo=None, cvt=None, phi=None):
#             return inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx
# 
#     def imex_res(vel, pres, dt, lastvel=None, othervel=None, phi=None):
#         if phi is None:
#             phi = dolfin.TestFunction(V)
# 
#         diffvel = .5*(vel+lastvel)  # Crank-Nicolson
#         diffrm = nu*inner(grad(diffvel)+grad(diffvel).T, grad(phi))*dx
#         if gradvsymmtrc:
#             nvec = dolfin.FacetNormal(V.mesh())
#             diffrm = diffrm - (nu*inner(grad(diffvel).T*nvec, phi))*outflowds
#         cnvfrm = convform(cvo=lastvel, cvt=othervel, phi=phi)
# 
#         pfrm = -inner(pres, div(phi))*dx
#         dtprt = 1./dt*dolfin.assemble(inner(vel, phi)*dx) \
#             - 1./dt*dolfin.assemble(inner(lastvel, phi)*dx)
#         res = dolfin.assemble(diffrm+cnvfrm+pfrm) + dtprt
# 
#         # import numpy as np
#         # nfc_c = dolfin.assemble(cnvfrm).get_local()
#         # print(np.linalg.norm(nfc_c), nfc_c[0], nfc_c.size)
#         # print('debggng')
#         return res
# 
#     return imex_res


import unittest
import numpy as np
import dolfin_navier_scipy.stokes_navier_utils as snu
import dolfin_navier_scipy.problem_setups as dnsps
import dolfin_navier_scipy.dolfin_to_sparrays as dts
from dolfin_navier_scipy.residual_checks import get_imex_res
            dnsps.get_sysmats(problem=self.problem, nu=self.nu,
                              bccontrol=False, charvel=self.charvel,
                              scheme=self.scheme, mergerhs=True,
                              meshparams=self.meshparams)

        # setting some parameters

        t0 = 0.0
        tE = .1
        Nts = 2
        tips = dict(t0=t0, tE=tE, Nts=Nts)

        soldict = stokesmatsc  # containing A, J, JT
        soldict.update(femp)  # adding V, Q, invinds, diribcs
        soldict.update(tips)  # adding time integration params
        soldict.update(fv=rhsd['fv'], fp=rhsd['fp'],
                       treat_nonl_explicit=True,
                       return_vp_dict=True,
                       no_data_caching=True,
                       start_ssstokes=True)

        vpdct = snu.solve_nse(**soldict)
        M, A, JT = stokesmatsc['M'], stokesmatsc['A'], stokesmatsc['JT']
        fv = rhsd['fv']
        V, invinds = femp['V'], femp['invinds']
        dt = (tE - t0) / Nts
        tm = (tE - t0) / 2

        reschkdict = dict(V=V, nu=self.nu,
                          gradvsymmtrc=True, outflowds=femp['outflowds'])

# Path: dolfin_navier_scipy/residual_checks.py
# def get_imex_res(V=None, outflowds=None, gradvsymmtrc=True, nu=None,
#                  implscheme='crni', explscheme='abtw'):
#     """ define the residual for an IMEX/AB2 time discretization
# 
#     """
#     if not implscheme == 'crni':
#         raise NotImplementedError()
# 
#     if explscheme == 'abtw':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (1.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx -
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'heun':
#         def convform(cvo=None, cvt=None, phi=None):
#             return (.5*inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx +
#                     .5*inner(dolfin.dot(cvt, nabla_grad(cvt)), phi)*dx)
# 
#     elif explscheme == 'eule':
#         def convform(cvo=None, cvt=None, phi=None):
#             return inner(dolfin.dot(cvo, nabla_grad(cvo)), phi)*dx
# 
#     def imex_res(vel, pres, dt, lastvel=None, othervel=None, phi=None):
#         if phi is None:
#             phi = dolfin.TestFunction(V)
# 
#         diffvel = .5*(vel+lastvel)  # Crank-Nicolson
#         diffrm = nu*inner(grad(diffvel)+grad(diffvel).T, grad(phi))*dx
#         if gradvsymmtrc:
#             nvec = dolfin.FacetNormal(V.mesh())
#             diffrm = diffrm - (nu*inner(grad(diffvel).T*nvec, phi))*outflowds
#         cnvfrm = convform(cvo=lastvel, cvt=othervel, phi=phi)
# 
#         pfrm = -inner(pres, div(phi))*dx
#         dtprt = 1./dt*dolfin.assemble(inner(vel, phi)*dx) \
#             - 1./dt*dolfin.assemble(inner(lastvel, phi)*dx)
#         res = dolfin.assemble(diffrm+cnvfrm+pfrm) + dtprt
# 
#         # import numpy as np
#         # nfc_c = dolfin.assemble(cnvfrm).get_local()
#         # print(np.linalg.norm(nfc_c), nfc_c[0], nfc_c.size)
#         # print('debggng')
#         return res
# 
#     return imex_res


# Path: tests/test_units_residuals.py
import unittest
import numpy as np
import dolfin_navier_scipy.stokes_navier_utils as snu
import dolfin_navier_scipy.problem_setups as dnsps
import dolfin_navier_scipy.dolfin_to_sparrays as dts
from dolfin_navier_scipy.residual_checks import get_imex_res

            dnsps.get_sysmats(problem=self.problem, nu=self.nu,
                              bccontrol=False, charvel=self.charvel,
                              scheme=self.scheme, mergerhs=True,
                              meshparams=self.meshparams)

        # setting some parameters

        t0 = 0.0
        tE = .1
        Nts = 2
        tips = dict(t0=t0, tE=tE, Nts=Nts)

        soldict = stokesmatsc  # containing A, J, JT
        soldict.update(femp)  # adding V, Q, invinds, diribcs
        soldict.update(tips)  # adding time integration params
        soldict.update(fv=rhsd['fv'], fp=rhsd['fp'],
                       treat_nonl_explicit=True,
                       return_vp_dict=True,
                       no_data_caching=True,
                       start_ssstokes=True)

        vpdct = snu.solve_nse(**soldict)
        M, A, JT = stokesmatsc['M'], stokesmatsc['A'], stokesmatsc['JT']
        fv = rhsd['fv']
        V, invinds = femp['V'], femp['invinds']
        dt = (tE - t0) / Nts
        tm = (tE - t0) / 2

        reschkdict = dict(V=V, nu=self.nu,
                          gradvsymmtrc=True, outflowds=femp['outflowds'])

        euleres = get_imex_res(explscheme='eule', **reschkdict)