Source code for eval_ns_forms

#!/usr/bin/env python
"""
Operators present in the FE discretization of (adjoint) Navier-Stokes terms.
"""
from __future__ import print_function
from __future__ import absolute_import
import sympy as s
from six.moves import range

[docs]def create_scalar(name, n_ep): vec = s.zeros(n_ep, 1) for ip in range(n_ep): vec[ip,0] = '%s%d' % (name, ip) return vec
[docs]def create_vector(name, n_ep, dim): """ordering is DOF-by-DOF""" vec = s.zeros(dim * n_ep, 1) for ii in range(dim): for ip in range(n_ep): vec[n_ep*ii+ip,0] = '%s%d%d' % (name, ii, ip) return vec
[docs]def create_scalar_base(name, n_ep): phi = s.zeros(1, n_ep) for ip in range(n_ep): phi[0,ip] = '%s%d' % (name, ip) return phi
[docs]def create_vector_base(name, phic, dim): n_ep = phic.shape[1] phi = s.zeros(dim, dim * n_ep) indx = [] for ii in range(dim): phi[ii,n_ep*ii:n_ep*(ii+1)] = phic indx.append(ii) return phi, indx
[docs]def create_scalar_base_grad(name, phic, dim): n_ep = phic.shape[1] gc = s.zeros(dim, n_ep) for ii in range(dim): for ip in range(n_ep): gc[ii,ip] = '%s%d%d' % (name, ii, ip) return gc
[docs]def create_vector_base_grad(name, gc, transpose=False): dim, n_ep = gc.shape g = s.zeros(dim * dim, dim * n_ep) indx = [] if transpose: for ir in range(dim): for ic in range(dim): g[dim*ir+ic,n_ep*ic:n_ep*(ic+1)] = gc[ir,:] indx.append((ic, ir)) else: for ir in range(dim): for ic in range(dim): g[dim*ir+ic,n_ep*ir:n_ep*(ir+1)] = gc[ic,:] indx.append((ir, ic)) return g, indx
[docs]def create_u_operator(u, transpose=False): dim = u.shape[0] op_u = s.zeros(dim * dim, dim) if transpose: for ir in range(dim): for ic in range(dim): op_u[dim*ir+ic,ic] = u[ir] else: for ii in range(dim): op_u[dim*ii:dim*(ii+1),ii] = u return op_u
[docs]def grad_vector_to_matrix(name, gv): dim2 = gv.shape[0] dim = int(s.sqrt(dim2)) gm = s.zeros(dim, dim) for ir in range(dim): for ic in range(dim): gm[ir,ic] = gv[dim*ir+ic,0] return gm
[docs]def substitute_continuous(expr, names, u, phi): pu = phi * u for ii in range(phi.rows): expr = expr.subs(pu[ii,0], names[ii]) return expr
[docs]def create_vector_var_data(name, phi, vindx, g, gt, vgindx, u): gu = g * u gum = grad_vector_to_matrix('gum', gu) print('g %s:\n' % name, gum) gut = gt * u gutm = grad_vector_to_matrix('gutm', gut) print('gt %s:\n' % name, gutm) pu = phi * u names = ['c%s%d' % (name, indx) for indx in vindx] cu = substitute_continuous(pu, names, u, phi) print('continuous %s:\n' % name, cu) gnames = ['cg%s%d_%d' % (name, indx[0], indx[1]) for indx in vgindx] cgu = substitute_continuous(gu, gnames, u, g) cgum = grad_vector_to_matrix('gum', cgu) print('continuous g %s:\n' % name, cgum) cgut = substitute_continuous(gut, gnames, u, g) cgutm = grad_vector_to_matrix('gutm', cgut) print('continuous gt %s:\n' % name, cgutm) op_u = create_u_operator(cu) print(op_u) op_ut = create_u_operator(cu, transpose=True) print(op_ut) out = { 'g' : gu, 'g_m' : gum, 'q' : pu, 'c' : cu, 'cg' : cgu, 'cg_m' : cgum, 'cgt' : cgut, 'cgt_m' : cgutm, 'op' : op_u, 'opt' : op_ut, 'names' : names, 'gnames' : gnames, } return out
[docs]def create_scalar_var_data(name, phi, g, u): gu = g * u pu = phi * u names = ['c%s' % name] cu = substitute_continuous(pu, names, u, phi) print('continuous %s:\n' % name, cu) gnames = ['cg%s_%d' % (name, ii) for ii in range(g.shape[0])] cgu = substitute_continuous(gu, gnames, u, g) print('continuous g %s:\n' % name, cgu) op_gu = create_u_operator(cgu) print(op_gu) out = { 'g' : gu, 'q' : pu, 'c' : cu, 'cg' : cgu, 'gop' : op_gu, 'names' : names, 'gnames' : gnames, } return out
[docs]def main(): n_ep = 3 dim = 2 u = create_vector('u', n_ep, dim) v = create_vector('v', n_ep, dim) b = create_vector('b', n_ep, dim) p = create_scalar('p', n_ep) q = create_scalar('q', n_ep) r = create_scalar('r', n_ep) print(u) print(v) phic = create_scalar_base('phic', n_ep) phi, vindx = create_vector_base('phi', phic, dim) gc = create_scalar_base_grad('gc', phic, dim) g, vgindx = create_vector_base_grad('g', gc) gt, aux = create_vector_base_grad('gt', gc, transpose=True) print(phi) print(phic) print(gc) print(g) print(gt) ud = create_vector_var_data('u', phi, vindx, g, gt, vgindx, u) vd = create_vector_var_data('v', phi, vindx, g, gt, vgindx, v) bd = create_vector_var_data('b', phi, vindx, g, gt, vgindx, b) pd = create_scalar_var_data('p', phic, gc, p) qd = create_scalar_var_data('q', phic, gc, q) rd = create_scalar_var_data('r', phic, gc, r) print(list(ud.keys())) assert bool(bd['op'].T * g == bd['opt'].T * gt) assert bool(bd['opt'].T * g == bd['op'].T * gt) assert bool(bd['cgt_m'] == bd['cg_m'].T) print('((b * grad) u), v)') form1 = vd['c'].T * bd['op'].T * ud['cg'] form2 = vd['c'].T * bd['opt'].T * ud['cgt'] print(form1) print(form2) print(bool(form1 == form2)) print('((v * grad) u), b)') form1 = vd['c'].T * bd['op'].T * ud['cgt'] form2 = vd['c'].T * bd['opt'].T * ud['cg'] print(form1) print(form2) print(bool(form1 == form2)) print('((u * grad) v), b)') form1 = vd['cgt'].T * bd['op'] * ud['c'] form2 = vd['cg'].T * bd['opt'] * ud['c'] print(form1) print(form2) print(bool(form1 == form2)) print('((b * grad) v), u)') form1 = vd['cg'].T * bd['op'] * ud['c'] form2 = vd['cgt'].T * bd['opt'] * ud['c'] print(form1) print(form2) print(bool(form1 == form2)) print('((v * grad) b), u)') form1 = vd['c'].T * bd['cgt_m'] * ud['c'] form2 = vd['c'].T * bd['cg_m'].T * ud['c'] print(form1) print(form2) print(bool(form1 == form2)) print('((b * grad) u), (b * grad) v)') form1 = vd['cg'].T * bd['op'] * bd['op'].T * ud['cg'] print(form1) print('((u * grad) b), (b * grad) v)') form1 = vd['cg'].T * bd['op'] * bd['cg_m'] * ud['c'] print(form1) print('(grad p, (b * grad) v)') form1 = vd['cg'].T * bd['op'] * pd['cg'] print(form1) print('(grad q, (b * grad) u)') form1 = qd['cg'].T * bd['op'].T * ud['cg'] print(form1) print('(grad q, (u * grad) b)') form1 = qd['cg'].T * bd['cg_m'] * ud['c'] print(form1) print('(grad r, (u * grad) v)') form1 = vd['cgt'].T * rd['gop'] * ud['c'] print(form1) return ud, vd, bd, pd, qd, rd
if __name__ == '__main__': ud, vd, bd, pd, qd, rd = main()