# 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

# import dolfin




class TimeIntResiduals(unittest.TestCase):

    def setUp(self):

        meshlvl = 1
        geodata = 'mesh/karman2D-rotcyl-bm_geo_cntrlbc.json'
        meshfile = 'mesh/karman2D-rotcyl_lvl{0}.xml.gz'.format(meshlvl)
        physregs = 'mesh/karman2D-rotcyl_lvl{0}_facet_region.xml.gz'.\
            format(meshlvl)
        self.meshparams = dict(strtomeshfile=meshfile,
                               movingwallcntrl=False,
                               strtophysicalregions=physregs,
                               strtobcsobs=geodata)
        self.nu = 1e-3
        self.charvel = 0.2
        self.scheme = 'TH'
        self.problem = 'gen_bccont'

    def test_residuals(self):

# 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


# import dolfin




class TimeIntResiduals(unittest.TestCase):

    def setUp(self):

        meshlvl = 1
        geodata = 'mesh/karman2D-rotcyl-bm_geo_cntrlbc.json'
        meshfile = 'mesh/karman2D-rotcyl_lvl{0}.xml.gz'.format(meshlvl)
        physregs = 'mesh/karman2D-rotcyl_lvl{0}_facet_region.xml.gz'.\
            format(meshlvl)
        self.meshparams = dict(strtomeshfile=meshfile,
                               movingwallcntrl=False,
                               strtophysicalregions=physregs,
                               strtobcsobs=geodata)
        self.nu = 1e-3
        self.charvel = 0.2
        self.scheme = 'TH'
        self.problem = 'gen_bccont'

    def test_residuals(self):

        femp, stokesmatsc, rhsd = \