From bbfcc2b522ad4830df755133d287cb228f61d516 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Oct 2015 15:30:38 -0500 Subject: [PATCH] Bring over contrib/ --- contrib/maxima/abarbanel-pml.mac | 128 +++ contrib/maxima/bgk-flow.mac | 83 ++ contrib/maxima/bilinearmap.mac | 40 + contrib/maxima/cns.mac | 45 + contrib/maxima/eclean.mac | 193 ++++ contrib/maxima/ecleanbase.mac | 38 + contrib/maxima/ehclean.mac | 29 + contrib/maxima/ehcleanbase.mac | 19 + contrib/maxima/em_units.mac | 37 + contrib/maxima/euler.mac | 172 ++++ contrib/maxima/gedney-pml.mac | 115 +++ contrib/maxima/maxwell.mac | 74 ++ contrib/maxima/maxwellbase.mac | 79 ++ contrib/maxima/myhelpers.mac | 142 +++ contrib/maxima/pml.mac | 73 ++ contrib/maxima/wave.mac | 99 ++ contrib/notes/dumka3.cpp | 393 ++++++++ contrib/notes/em-cheat.tm | 50 + contrib/notes/fluxfan.fig | 40 + contrib/notes/hedge-notes.tm | 1530 ++++++++++++++++++++++++++++++ contrib/notes/ip-primal-nu.tm | 89 ++ contrib/notes/ip-primal.tm | 87 ++ contrib/notes/mystyle.ts | 32 + 23 files changed, 3587 insertions(+) create mode 100644 contrib/maxima/abarbanel-pml.mac create mode 100644 contrib/maxima/bgk-flow.mac create mode 100644 contrib/maxima/bilinearmap.mac create mode 100644 contrib/maxima/cns.mac create mode 100644 contrib/maxima/eclean.mac create mode 100644 contrib/maxima/ecleanbase.mac create mode 100644 contrib/maxima/ehclean.mac create mode 100644 contrib/maxima/ehcleanbase.mac create mode 100644 contrib/maxima/em_units.mac create mode 100644 contrib/maxima/euler.mac create mode 100644 contrib/maxima/gedney-pml.mac create mode 100644 contrib/maxima/maxwell.mac create mode 100644 contrib/maxima/maxwellbase.mac create mode 100644 contrib/maxima/myhelpers.mac create mode 100644 contrib/maxima/pml.mac create mode 100644 contrib/maxima/wave.mac create mode 100644 contrib/notes/dumka3.cpp create mode 100644 contrib/notes/em-cheat.tm create mode 100644 contrib/notes/fluxfan.fig create mode 100644 contrib/notes/hedge-notes.tm create mode 100644 contrib/notes/ip-primal-nu.tm create mode 100644 contrib/notes/ip-primal.tm create mode 100644 contrib/notes/mystyle.ts diff --git a/contrib/maxima/abarbanel-pml.mac b/contrib/maxima/abarbanel-pml.mac new file mode 100644 index 00000000..879e8363 --- /dev/null +++ b/contrib/maxima/abarbanel-pml.mac @@ -0,0 +1,128 @@ +kill(all); +load("eigen"); +load("itensor"); +load("diag"); + +/* +See + +S. Abarbanel und D. Gottlieb, “On the construction and analysis of absorbing +layers in CEM,” Applied Numerical Mathematics, vol. 27, 1998, S. 331-340. +(eq 3.7-3.11) + +E. Turkel und A. Yefet, “Absorbing PML +boundary layers for wave-like equations,” +Applied Numerical Mathematics, vol. 27, +1998, S. 533-557. +(eq. 4.10) + +*/ + +/* -------------------------------------------------------------------------- */ +/* Variable declarations */ +/* -------------------------------------------------------------------------- */ + +coords:[x,y,z]; +allvars:append([t],coords); + +max_E:makelist(concat(E,i),i,coords); +max_H:makelist(concat(H,i),i,coords); +max_w:append(max_E,max_H); + +max_P:makelist(concat(P,i),i,coords); +max_Q:makelist(concat(Q,i),i,coords); +aux_w:append(max_P, max_Q); + +depends(max_w,allvars); +depends(aux_w,allvars); + +sig:makelist(concat(s,i),i,coords); + +depends(sig, coords); + +/* -------------------------------------------------------------------------- */ +/* Utilities */ +/* -------------------------------------------------------------------------- */ + +crossfunc(f):=makelist( + sum(sum( + levi_civita([i,j,k])*f(j,k), + j,1,3),k,1,3),i,1,3)$ + +curl(x):=crossfunc(lambda([i,j], diff(x[j], coords[i]))); + +shift(l, amount):=makelist(l[mod(i-amount-1, length(l))+1], i, 1, length(l)); + +norm(x):=sum(x[i]^2, i, 1, length(x)); + + +/* -------------------------------------------------------------------------- */ +/* Operator building */ +/* -------------------------------------------------------------------------- */ + +/* +This here is *NOT* in conservation form + u_t + A u_x = 0, +but in ODE form + u_t = (- A) u_x. + ---------- rhs +*/ + +max_rhs:append(1/epsilon*curl(max_H),-1/mu*curl(max_E)); + + +pml_rhs:makelist(0,i,1,6); +aux_eqn:makelist(0,i,1,6); + +for mx:1 thru 3 do + block([my,mz,sgn], + my:mod(mx-1+1,3)+1, + mz:mod(mx-1+2,3)+1, + assert(levi_civita([mx,my,mz])=1), + pml_rhs[mx]:pml_rhs[mx] + +1/epsilon*diff(max_H[mz],coords[my]) + -sig[my]/epsilon*(2*max_E[mx]+max_P[mx]) + , + pml_rhs[my]:pml_rhs[my] + -1/epsilon*diff(max_H[mz],coords[mx]) + -sig[mx]/epsilon*(2*max_E[my]+max_P[my]) + , + pml_rhs[3+mz]:pml_rhs[3+mz] + +1/mu*diff(max_E[mx],coords[my]) + -1/mu*diff(max_E[my],coords[mx]) + +1/mu*diff(sig[mx]/epsilon,coords[mx])*max_Q[mx] + +1/mu*diff(sig[my]/epsilon,coords[my])*max_Q[my] + , + aux_eqn[mx]:aux_eqn[mx] + -diff(max_P[mx],t) + +sig[my]/epsilon*max_E[mx] + , + aux_eqn[my]:aux_eqn[my] + +sig[mx]/epsilon*max_E[my] + , + aux_eqn[3+mx]:aux_eqn[3+mx] + -diff(max_Q[mx],t)-sig[mx]/epsilon*max_Q[mx] + -max_E[my] -max_E[mz] + ); + +pml_eqn:makelist(diff(max_w[i],t)=pml_rhs[i], i, 1, 6); +aux_eqn:makelist(aux_eqn[i]=0, i, 1, 6); +print(expand(covect(pml_eqn))); +print(expand(covect(aux_eqn))); +slist:[ + /*Qx=-Q[y],Qy=Q[x],*/Qz=0, + /*Px=P[y],Py=-P[x],*/Pz=0, + sz=0, + Hx=0,Hy=0,Ez=0, + epsilon=1,mu=1]; + +print(expand(covect(subst(slist,pml_eqn)))); +print(expand(covect(subst(slist,aux_eqn)))); + +load("em_units.mac"); +u_P:u_sigma*u_E/u_epsilon0*u_t; +u_Q:u_E*u_t; + +assert(u_E/u_t - u_sigma/u_epsilon0*u_P=0); +assert(u_H/u_t - (1/u_mu0)*(u_sigma/u_epsilon0)/u_x*u_Q=0); + diff --git a/contrib/maxima/bgk-flow.mac b/contrib/maxima/bgk-flow.mac new file mode 100644 index 00000000..72cbbf8f --- /dev/null +++ b/contrib/maxima/bgk-flow.mac @@ -0,0 +1,83 @@ +kill(all); + +d:2; + +load("myhelpers.mac"); +assume(sqrt_rt>0); +assume(tau>0); +assume(r>0); + +coords:[x,y,z]; +n:makelist(concat(n, coords[k]), k, 1, d); + +sqrt_rt:r; +/* +w0:a1; +w1:a2; +w2:a3; +w3:a4; +w4:a5; +w5:a6; +*/ + +/* ------------------------------------------------------------------------ */ +/* do not edit -- automatically generated by BGKFlowOperator */ + +d:2; + +vars: [w0, w1, w2, w3, w4, w5]; + +depends(vars, coords); + +bgk_neg_rhss: [ +sqrt_rt*(diff(w1, x) + diff(w2, y)), +sqrt_rt*(diff(w0, x) + diff(w3, y) + sqrt(2)*diff(w4, x)), +sqrt_rt*(diff(w0, y) + diff(w3, x) + sqrt(2)*diff(w5, y)), +sqrt_rt*(diff(w1, y) + diff(w2, x)) + 1/tau*(w3 + (-1)*w1*w2/w0), +sqrt(2)*sqrt_rt*diff(w1, x) + 1/tau*(w4 + (-1)*w1*w1/(sqrt(2)*w0)), +sqrt(2)*sqrt_rt*diff(w2, y) + 1/tau*(w5 + (-1)*w2*w2/(sqrt(2)*w0)) +]; + +/* ------------------------------------------------------------------------ */ + + +bgk_A:genmatrix( + lambda([i,j], + sum(n[k]*diff(bgk_neg_rhss[i], diff(vars[j], coords[k])), k, 1, d)), + length(vars), length(vars)); + +/* diagonalize system ------------------------------------------------------- */ +[bgk_V, bgk_D, bgk_invV]:hypdiagonalize(bgk_A); + +/* auxiliary variables ------------------------------------------------------ */ +bgk_wm:makelist(concat(vars[i], m), i, 1, length(vars)); +bgk_wp:makelist(concat(vars[i], p), i, 1, length(vars)); +bgk_sminw:hypsimp(bgk_invV.bgk_wm); +bgk_spinw:hypsimp(bgk_invV.bgk_wp); + +/* rankine-hugoniot flux ---------------------------------------------------- */ +bgk_sflux:hyp_upwind_flux( + [-sqrt(3)*sqrt_rt,-sqrt_rt,0,sqrt_rt,sqrt(3)*sqrt_rt], + bgk_D); + +bgk_wflux:fullhypsimp(bgk_V.ev(bgk_sflux, [sm=bgk_sminw,sp=bgk_spinw])); + +bgk_stronglocalpart:subst(r=rm, bgk_A).bgk_wm; + +bgk_strongwflux:fullhypsimp(bgk_stronglocalpart-bgk_wflux); + +bgk_strongwflux:fullhypsimp(subst([rm=r,rp=r], bgk_strongwflux)); + +/* display flux expression -------------------------------------------------- */ +display2d:false; + +dispsubst:append( + makelist(concat(w, i, p)=w[i][ext], i, 0, length(vars)-1), + makelist(concat(w, i, m)=w[i][int], i, 0, length(vars)-1), + makelist(n[i]=normal[i-1], i, 1, d) + ); + +print(makelist( + subst(dispsubst, bgk_strongwflux[i,1]), + i, 1, length(bgk_strongwflux))); + diff --git a/contrib/maxima/bilinearmap.mac b/contrib/maxima/bilinearmap.mac new file mode 100644 index 00000000..41a2ba26 --- /dev/null +++ b/contrib/maxima/bilinearmap.mac @@ -0,0 +1,40 @@ +kill(all); +lflatten(l):=apply(append, l); +dims:3; +origpoints: [[-1],[1]]; +for i:2 thru dims do + origpoints:lflatten( + makelist( + makelist( + append(opi, [j]) + ,j, [-1, 1] + ), + opi, origpoints)); +points : makelist(ascii(97+i), i, 0, length(origpoints)-1); + +vars:makelist(concat(v,ascii(97+i)), i, 0, dims-1); +mapargs:append( + [1], vars, + lflatten( + makelist( + makelist(vars[i]*vars[j], j, i+1, length(vars)), + i, 1, length(vars) + ) + ) +); + +/* Idea: (x,y)^T = MAT*(1,u,v,uv) */ +mapmat:genmatrix(mm, dims, length(mapargs)); + +f:mapmat.mapargs; + +myjacobian(f, vars):=genmatrix( + lambda([i,j], diff(f[i], vars[j])), + length(f), length(vars)); + +fjac:myjacobian(f, vars); + +print("bilinear map jacobian:"); +print(fjac); +print("bilinear map jacobian det:"); +print(expand(determinant(fjac))); diff --git a/contrib/maxima/cns.mac b/contrib/maxima/cns.mac new file mode 100644 index 00000000..7bf5640f --- /dev/null +++ b/contrib/maxima/cns.mac @@ -0,0 +1,45 @@ +/* compare formulation in JSH/TW with + * http://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Compressible_flow_of_Newtonian_fluids + */ + +kill(all); +d:2; + +div(vec):=sum(diff(vec[i], coords[i]), i, 1, length(vec)); +grad(f):=makelist(diff(f, coords[i]), i, 1, d); +vlaplace(vec):=makelist(div(grad(vec[i])), i, 1, length(vec)); + +uvec:makelist([u,v,w][i], i, 1, d); +coords:makelist([x,y,z][i], i, 1, d); + +depends(uvec, coords); + +dudx:genmatrix( + lambda([i,j], diff(uvec[i], coords[j])), + d, d); + +delta(i, j):=if i = j then 1 else 0; + +tau:mu*genmatrix( + lambda([i,j], + dudx[i,j] + dudx[j,i] - 2/3 * delta(i,j) * mat_trace(dudx)), + d, d); + +rhou_rhs:makelist(div(tau[i]), i, 1, d); + +muv:0; +rhou_rhs2:(1/3*mu+muv)*grad(div(uvec)) + mu*vlaplace(uvec); + +/* +print(rhou_rhs); +print(rhou_rhs2); +*/ +print(ratsimp(rhou_rhs-rhou_rhs2)); + +e_rhs:div(uvec.tau); + +e_rhs2:(1/3*mu+muv)*grad(div(uvec)) + mu*makelist(div(uvec*grad(vec[i])), i, 1, d); + +print(e_rhs); +print(e_rhs2); +print(ratsimp(e_rhs-e_rhs2)); diff --git a/contrib/maxima/eclean.mac b/contrib/maxima/eclean.mac new file mode 100644 index 00000000..3c82f081 --- /dev/null +++ b/contrib/maxima/eclean.mac @@ -0,0 +1,193 @@ +kill(all); + +load("ecleanbase.mac"); + +print(extract_diagonal(eclean_D)); + +/* ------------------------------------------------------------------------- */ +/* flux */ +/* ------------------------------------------------------------------------- */ + +eclean_Dinc:ratsubst(c,1/(sqrt(%epsilon)*sqrt(%mu)), eclean_D); +eclean_sflux:hyp_upwind_flux([-%chi*c,-c,0,c,%chi*c], eclean_Dinc); + +/* +print("e-clean system flux in terms of characteristic variables:"); +print(eclean_sflux); +*/ + +/* FIXME: eclean_V should not depend on epsilon and mu, but it does + For now, make cp and cm equal. */ + +eclean_sflux:subst( + [cp=1/(sqrt(%epsilon)*sqrt(%mu)), + cm=1/(sqrt(%epsilon)*sqrt(%mu)), + %chip=%chi, + %chim=%chi], + eclean_sflux); + +eclean_wflux:fullhypsimp(eclean_V.ev(eclean_sflux, + [sm=eclean_sminw,sp=eclean_spinw])); + +/* +print("e-clean system weak flux in terms of physical variables:"); +print(eclean_wflux); +*/ + +eclean_strongwflux:eclean_A.eclean_wm - eclean_wflux; + +print("e-clean system strong flux in terms of physical variables:"); +print(eclean_strongwflux); + +eclean_maxstrongdiff:fullhypsimp( + eclean_strongwflux + -append(max_strongwflux,matrix([0])) + ); + +/* +print("e-clean additional strong-form flux wrt Maxwell system:"); +print(eclean_maxstrongdiff); +*/ + +eclean_simple_maxstrongdiff:1/2*append( + %chi*c*n*(%phi[m]-%phi[p] - n.(max_Em-max_Ep)), + [0,0,0], + [c*%chi*(-(%phi[m]-%phi[p]) + n.(max_Em-max_Ep))]); + +assert(norm_2_squared(hypsimp( + eclean_maxstrongdiff + - ev(eclean_simple_maxstrongdiff, [c=1/sqrt(%epsilon*%mu)])))=0); + +/* ------------------------------------------------------------------------- */ +/* Radiation BCs */ +/* ------------------------------------------------------------------------- */ + +eclean_radbdryspinw:makelist( + if eclean_D[i,i] >= 0 then eclean_sminw[i,1] else 0, +i, 1, length(eclean_D))$ + +/* +print("Radiation boundary condition for E-divclean system:"); +print(fullhypsimp(eclean_V.eclean_radbdryspinw)); +*/ + +/* ------------------------------------------------------------------------- */ +/* Better PEC */ +/* ------------------------------------------------------------------------- */ +/* PEC means: n x E = 0. (no tangential field) + +For normal Maxwell PEC, you prescribe E+ = - E-, even though the +normal component is wrong (it should be equal to the E- normal component). +That doesn't matter because the Maxwell flux only looks at the tangential +component. But here, it suddenly ends up mattering, so we need to prescribe +the normal component correctly. +*/ + +eclean_simplepecwp:covect(append(max_pecbdrywp, [%phi[m]])); +eclean_betterpecwp:covect(append( + normalcomp(max_Em)-tangentialcomp(max_Em), max_Hm, [-%phi[m]])); +eclean_betterpecsp:fullhypsimp( + eclean_invV.ev(eclean_betterpecwp, + [Em=eclean_Emins, Hm=eclean_Hmins, %phi[m]=eclean_phimins])); + +/* +print("Better PEC condition in characteristic variables:"); +print(eclean_betterpecsp); +*/ + +print("Better PEC condition in physical variables:"); +print(eclean_betterpecwp); + + +eclean_flux_at_betterpec:fullhypsimp(ev(append(max_strongwflux,matrix([0]))/*eclean_strongwflux*/, + [Ep=subrange(eclean_betterpecwp,1,3), + Hp=subrange(eclean_betterpecwp,4,6), + %phi[p]=eclean_betterpecwp[7,1]])); + +eclean_flux_at_simplepec:fullhypsimp(ev(append(max_strongwflux,matrix([0]))/*eclean_strongwflux*/, + [Ep=subrange(eclean_simplepecwp,1,3), + Hp=subrange(eclean_simplepecwp,4,6), + %phi[p]=eclean_simplepecwp[7,1]])); + +assert(norm_2_squared(hypsimp(eclean_flux_at_betterpec-eclean_flux_at_simplepec))=0); + +/* ------------------------------------------------------------------------- */ + +eclean_munzpecwithphi(phival):=fullhypsimp(eclean_invV. + vstack(vstack(-eclean_Emins, eclean_Hmins), [phival])); + +/* +print("Munz et al PEC boundary with phi+=0:"); +print(eclean_munzpecwithphi(0)); +print("Munz et al PEC boundary with phi+=phi-:"); +print(eclean_munzpecwithphi(eclean_phimins)); +print("Munz et al PEC boundary with phi+=-phi-:"); +print(eclean_munzpecwithphi(-eclean_phimins)); +*/ + +/* ------------------------------------------------------------------------- */ +/* chi-radiation + PEC BC */ +/* ------------------------------------------------------------------------- */ +eclean_chiradbdryspinw:[ + -eclean_sminw[3,1], + -eclean_sminw[4,1], + -eclean_sminw[1,1], + -eclean_sminw[2,1], + 0*eclean_sminw[5,1], + eclean_sminw[6,1], + eclean_sminw[7,1] + ]$ + +eclean_chiradwp:fullhypsimp(eclean_V.eclean_chiradbdryspinw); + +print("PEC+chirad BC for div E cleaning system"); +print(eclean_chiradwp); + +eclean_simple_expr_chiradwp:append( + -max_Em+3/2*n*(n.max_Em) + 1/2*%phi[m]*n, + max_Hm, + 1/2*[%phi[m]+max_Em.n] + ); + +eclean_diff_chiradwp:hypsimp(eclean_chiradwp + - subst([c=1/sqrt(%epsilon*%mu)], eclean_simple_expr_chiradwp)); + +assert(norm_2_squared(eclean_diff_chiradwp)=0); + +/* +print("Limit of PEC+chirad BC as %phi, %chi -> 0"); +print(fullhypsimp(limit(limit( + subst([%phi[m]=phi],eclean_chiradwp),phi,0),%chi,0))); + */ + +eclean_strongw_chirad_flux:fullhypsimp(ev(eclean_strongwflux, + [Ep=subrange(eclean_chiradwp,1,3), + Hp=subrange(eclean_chiradwp,4,6), + %phi[p]=eclean_chiradwp[7,1]])); +print("Flux at PEC+chirad:"); +print(eclean_strongw_chirad_flux); + +eclean_strongw_pec_flux:fullhypsimp(ev(eclean_strongwflux, + [Ep=subrange(eclean_betterpecwp,1,3), + Hp=subrange(eclean_betterpecwp,4,6), + %phi[p]=eclean_betterpecwp[7,1]])); +print("Flux at pure PEC:"); +print(eclean_strongw_pec_flux); + +eclean_strongw_flux_diff:fullhypsimp( + eclean_strongw_pec_flux-eclean_strongw_chirad_flux); + +print("PEC Flux - Chirad Flux"); +print(eclean_strongw_flux_diff); + + +/* +f1:subst([ + Normal(0)=nx,Normal(1)=ny,Normal(2)=nz, + Int[0]=Em[1],Int[1]=Em[2],Int[2]=Em[3], + Ext[0]=Ep[1],Ext[1]=Ep[2],Ext[2]=Ep[3], + Int[3]=Hm[1],Int[4]=Hm[2],Int[5]=Hm[3], + Ext[3]=Hp[1],Ext[4]=Hp[2],Ext[5]=Hp[3], + Int[6]=%phi[m],Ext[6]=%phi[p]] + ,f0) +*/ diff --git a/contrib/maxima/ecleanbase.mac b/contrib/maxima/ecleanbase.mac new file mode 100644 index 00000000..6996c0e9 --- /dev/null +++ b/contrib/maxima/ecleanbase.mac @@ -0,0 +1,38 @@ +load("maxwellbase.mac"); + +assume(%chi>0); + +/* +eclean_Asimp:blockmat( + max_Asimp, + vstack(epsinv*muinv*covect(n),constmatrix(3,1,0)), + hstack(%chi^2*n,constmatrix(1,3,0)), + zeromatrix(1,1) +); +*/ + +eclean_Asimp:blockmat( + max_Asimp, + vstack(%chi*sqrt(epsinv*muinv)*covect(n),constmatrix(3,1,0)), + hstack(%chi*sqrt(epsinv*muinv)*n,constmatrix(1,3,0)), + zeromatrix(1,1) +); + +[eclean_V, eclean_D, eclean_invV]:max_invsubst(hypdiagonalize(eclean_Asimp)); +eclean_A:max_invsubst(eclean_Asimp); + +eclean_wm:vstack(max_wm,[%phi[m]]); +eclean_wp:vstack(max_wp,[%phi[p]]); + +eclean_sm:makelist(sm[i],i,1,length(eclean_D)); +eclean_sp:makelist(sp[i],i,1,length(eclean_D)); + +eclean_sminw:hypsimp(eclean_invV.eclean_wm); +eclean_spinw:hypsimp(eclean_invV.eclean_wp); + +eclean_wmins:hypsimp(eclean_V.eclean_sm); +eclean_wpins:hypsimp(eclean_V.eclean_sp); + +eclean_Emins:makelist(eclean_wmins[i,1],i,1,3); +eclean_Hmins:makelist(eclean_wmins[i,1],i,4,6); +eclean_phimins:eclean_wmins[7,1]; diff --git a/contrib/maxima/ehclean.mac b/contrib/maxima/ehclean.mac new file mode 100644 index 00000000..0d5b608e --- /dev/null +++ b/contrib/maxima/ehclean.mac @@ -0,0 +1,29 @@ +kill(all); +load("ehcleanbase.mac"); + +/* ------------------------------------------------------------------------- */ +/* flux */ +/* ------------------------------------------------------------------------- */ + +ehclean_Dinc:ratsubst(c,1/(sqrt(%epsilon)*sqrt(%mu)), ehclean_D); +print(ehclean_Dinc); +ehclean_sflux:hyp_upwind_flux([-%chi*c,-c,c,%chi*c], ehclean_Dinc); + +print("eh-clean system flux in terms of characteristic variables:"); +print(covect(ehclean_sflux)); + +/* FIXME: ehclean_V should not depend on epsilon and mu, but it does + For now, make cp and cm equal. */ +/* +ehclean_sflux:subst( + [cp=1/(sqrt(%epsilon)*sqrt(%mu)), + cm=1/(sqrt(%epsilon)*sqrt(%mu)), + %chip=%chi, + %chim=%chi], + ehclean_sflux); + +print("e-clean system flux in terms of physical variables:"); +ehclean_wflux:fullhypsimp(ehclean_V.ev(ehclean_sflux, + [sm=ehclean_sminw,sp=ehclean_spinw])); +print(ehclean_wflux); +*/ diff --git a/contrib/maxima/ehcleanbase.mac b/contrib/maxima/ehcleanbase.mac new file mode 100644 index 00000000..70dd8410 --- /dev/null +++ b/contrib/maxima/ehcleanbase.mac @@ -0,0 +1,19 @@ +load("maxwellbase.mac"); + +assume(%chi>0); + +ehclean_Asimp:blockmat( + max_Asimp, + hstack( + vstack(epsinv*muinv*covect(n),zeromatrix(3,1)), + vstack(zeromatrix(3,1),epsinv*muinv*covect(n)) + ), + + vstack( + hstack(%chi^2*n,zeromatrix(1,3)), + hstack(zeromatrix(1,3),%chi^2*n) + ), + zeromatrix(2,2) +); +ehclean_A:max_invsubst(ehclean_Asimp); +[ehclean_V, ehclean_D, ehclean_invV]:max_invsubst(hypdiagonalize(ehclean_Asimp)); diff --git a/contrib/maxima/em_units.mac b/contrib/maxima/em_units.mac new file mode 100644 index 00000000..dabce501 --- /dev/null +++ b/contrib/maxima/em_units.mac @@ -0,0 +1,37 @@ +kill(all); +load("myhelpers.mac"); +assume(m>0); +assume(N>0); +assume(s>0); + +N:kg*m/s^2; + +C : A*s; +J : N*m; +V : J/C; +F : C/V; +T : V*s/m^2; +N : kg*m/s^2; + +u_c0 : m/s; +u_mu0 : N / A^2; +u_epsilon0 : 1/(u_c0^2*u_mu0); + + +c0 : 299792458 * u_c0; +mu0 : 4*%pi * 10^(-7) * u_mu0; +epsilon0 : 1/(c0^2*mu0); + +u_E : V/m; +u_D : C/m^2; +u_B : T; +u_H : A/m; +u_J : A/m^2; +u_t : s; +u_x : m; + +u_sigma : u_J/u_E; + +u_ndE : u_E; +u_ndH : sqrt(u_mu0/u_epsilon0) * u_H; + diff --git a/contrib/maxima/euler.mac b/contrib/maxima/euler.mac new file mode 100644 index 00000000..77799424 --- /dev/null +++ b/contrib/maxima/euler.mac @@ -0,0 +1,172 @@ +kill(all); + +load("myhelpers.mac"); + +assume(%gamma>1); + +d:2; +usyms:[u,v,w]; +coords:[x,y,z]; +uvec:makelist(usyms[i], i, 1, d); +rhouvec:rho*uvec; +E_expr:p/(%gamma-1)+rho/2*uvec.uvec; + +p_solved:rhs(solve(E=E_expr,p)[1]); +p_used: p; /* p or p_solved */ + +/* fluxes ------------------------------------------------------------------- */ +vars:append([rho, E], rhouvec); +depends(append([rho, E, p], uvec), [x,y,z,t]); + +rho_flux:rhouvec; +E_flux:uvec*(E+p_used); +rhou_flux:makelist(makelist( + rhouvec[i]*uvec[j] + if i = j then p_used else 0, + j, 1, d), i, 1, d); + +all_fluxes:makelist( + vstack([rho_flux[i]], [E_flux[i]], rhou_flux[i]), + i, 1, d); + +euler_eqns:makelist( + -'diff(vars[nr],t)=sum('diff(all_fluxes[i][nr], coords[i]), i, 1, d), + nr, 1, length(vars)); + +/* linearization ------------------------------------------------------------ */ +u0vec:makelist(concat(usyms[i], 0), i, 1, d); +duvec:makelist(concat('d, usyms[i]), i, 1, d); +drhouvec:rho0*duvec+drho*u0vec; + +assume(%gamma>1); +assume(p0>0); +assume(rho0>0); + +zero_subst:append( + [rho=rho0, p=p0], + makelist(usyms[i] =concat(usyms[i], 0), i, 1, d)); +lin_subst:append( + [rho=rho0+drho, p=p0+dp], + makelist(usyms[i] =concat(usyms[i], 0) + duvec[i], i, 1, d)); + +E0_expr:subst(zero_subst,E_expr); +dE:subst(lin_subst,E_expr)-E0_expr; + +lin_subst:append(lin_subst, [E=E0+dE]); + +kill_second_order_terms(x):=block( + for lv1 in all_lin_vars do + for lv2 in all_lin_vars do + block( + x:ratsubst(0,lv1*lv2, x), + for co in coords do + x:ratsubst(0,diff(lv1, co)*lv2, x)), + x); + +lin_vars:append([drho, dE], drhouvec); +all_lin_vars:append(duvec, [drho, dp]); + +depends(all_lin_vars, [x,y,z,t]); + +lin_euler_eqns:kill_second_order_terms(makelist( + /* ref solution is time-invariant*/ + -diff(lin_vars[nr],t) + =sum(diff(subst(lin_subst, all_fluxes[i][nr]), coords[i]), i, 1, d), + nr, 1, length(vars))); + +lin_euler_eqns:ratsimp(lin_euler_eqns); + +/* convert to primitive variables */ +ddrho_dt: rhs(solve(lin_euler_eqns[1],diff(drho,t))[1]); + +lin_euler_eqns_p:makelist(lin_euler_eqns[i], i, 1, length(lin_euler_eqns)); + +for i:2 thru d+2 do + lin_euler_eqns_p[i]:subst(diff(drho,t)=ddrho_dt, lin_euler_eqns_p[i]), + +for i:1 thru d do + lin_euler_eqns_p[2+i]:-solve(lin_euler_eqns_p[2+i], diff(duvec[i], t))[1]; +dduvec_dt: makelist( + rhs(solve(lin_euler_eqns_p[i+2],diff(duvec[i],t))[1]), + i, 1, d); + +lin_euler_eqns_p[2]:subst(E0=E0_expr, lin_euler_eqns_p[2]); +for i:1 thru d do + lin_euler_eqns_p[2]:subst(diff(duvec[i],t)=dduvec_dt[i], lin_euler_eqns_p[2]); +lin_euler_eqns_p[2]:-solve(lin_euler_eqns_p[2], diff(dp, t))[1]; + +lin_euler_eqns_p:kill_second_order_terms(lin_euler_eqns_p); + +/* matrix building ---------------------------------------------------------- */ +prim_lin_vars:append([drho,dp],duvec); +n:makelist(concat(n, coords[k]), k, 1, d); +euler_mat:genmatrix( + lambda([i,j], + sum(n[k] + *diff(rhs(lin_euler_eqns_p[i]), diff(prim_lin_vars[j], coords[k])), + k, 1, d)), + length(lin_vars), length(lin_vars)); + +/* diagonalize system ------------------------------------------------------- */ +[euler_V, euler_D, euler_invV]:hypdiagonalize(euler_mat); + +rel_D:ratsimp(euler_D - n.u0vec*ident(d+2)); + +/* auxiliary variables, external and internal states ------------------------ */ +c0:sqrt(%gamma*p0/rho0); + +euler_wm:makelist(concat(prim_lin_vars[i], m), i, 1, length(prim_lin_vars)); +euler_wp:makelist(concat(prim_lin_vars[i], p), i, 1, length(prim_lin_vars)); + +euler_sminw:hypsimp(euler_invV.euler_wm); +euler_spinw:hypsimp(euler_invV.euler_wp); + +dumvec:makelist(euler_wm[i+2], i, 1, d); +dupvec:makelist(euler_wp[i+2], i, 1, d); + +/* subsonic outflow bc------------------------------------------------------- */ +euler_outflowbdryspinw:makelist( + /* outer state equals free-stream flow, about which we have linearized. + Hence the linearized state variable, which is the difference to free-stream + flow, is set to zero. */ + /* + = 0: convection: from interior + > 0: supersonic outflow: from interior + < 0: supersonic inflow: from exterior + */ + if rel_D[i,i] >= 0 then euler_sminw[i,1] else 0, + i, 1, d+2); + +euler_woutflowbdry:fullhypsimp(euler_V.euler_outflowbdryspinw); + +euler_known_woutflowbdry:vstack([ + drhom + n.dumvec*rho0/(2*c0) - dpm/(2*c0^2), + c0*rho0*(n.dumvec)/2 + dpm/2], + dumvec - n*(n.dumvec)/2 + dpm*n/(2*c0*rho0) + ); + +euler_diff_woutflowbdry:hypsimp(euler_woutflowbdry-euler_known_woutflowbdry); +assert(norm_2_squared(euler_diff_woutflowbdry)=0); + +/* subsonic inflow bc ------------------------------------------------------- */ +euler_inflowbdryspinw:makelist( + /* outer state equals free-stream flow, about which we have linearized. + Hence the linearized state variable, which is the difference to free-stream + flow, is set to zero. */ + /* + > 0: supersonic outflow: from interior + = 0: convection: from exterior + < 0: supersonic inflow: from exterior + */ + if rel_D[i,i] <= 0 then 0 else euler_sminw[i,1], + i, 1, d+2); + +euler_winflowbdry:fullhypsimp(euler_V.euler_inflowbdryspinw); + +euler_known_winflowbdry:vstack([ + n.dumvec*rho0/(2*c0) + dpm/(2*c0^2), + c0*rho0*(n.dumvec)/2 + dpm/2], + n*(n.dumvec)/2 + dpm*n/(2*c0*rho0) + ); + +euler_diff_winflowbdry:hypsimp(euler_winflowbdry-euler_known_winflowbdry); +assert(norm_2_squared(euler_diff_winflowbdry)=0); diff --git a/contrib/maxima/gedney-pml.mac b/contrib/maxima/gedney-pml.mac new file mode 100644 index 00000000..8b4a3ed1 --- /dev/null +++ b/contrib/maxima/gedney-pml.mac @@ -0,0 +1,115 @@ +kill(all); +load("maxwellbase.mac"); + +/* +See +S.D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for +the truncation of FDTD lattices,” IEEE Transactions on Antennas and +Propagation, vol. 44, 1996, S. 1630-1639. +*/ + +/* -------------------------------------------------------------------------- */ +/* Variable declarations */ +/* -------------------------------------------------------------------------- */ + +coords:[x,y,z]; +allvars:append([t],coords); + +max_E:makelist(concat(E,i),i,coords); +max_H:makelist(concat(H,i),i,coords); +max_w:append(max_E,max_H); + +max_D:makelist(concat(D,i),i,coords); +max_B:makelist(concat(B,i),i,coords); +aux_w:append(max_D, max_B); + +depends(max_w,allvars); +depends(aux_w,allvars); + +/* -------------------------------------------------------------------------- */ +/* Utilities */ +/* -------------------------------------------------------------------------- */ + +curl(x):=crossfunc(lambda([i,j], diff(x[j], coords[i]))); + +make_eqns(l, r):= makelist(l[i]=r[i], i, 1, length(l)); + +shift(l, amount):=makelist(l[mod(i-amount-1, length(l))+1], i, 1, length(l)); + +norm(x):=sum(x[i]^2, i, 1, length(x)); + +/* -------------------------------------------------------------------------- */ +/* Operator building */ +/* -------------------------------------------------------------------------- */ + +sigexpr(c):=(1+sigma[c]/(%i*omega*epsilon)); + +/* +This here is *NOT* in conservation form + u_t + A u_x = 0, +but in ODE form + u_t = (- A) u_x. + ---------- rhs +*/ + +max_rhs:append(curl(max_H),-curl(max_E)); + +s:makelist(sigexpr(c), c, coords); +sl:shift(s, -1); +sr:shift(s, 1); +sdiv:sr/s; + +max_D_values:epsilon*sdiv*max_E; +max_B_values:mu*sdiv*max_H; + +main_eqns:expand(make_eqns( + %i*omega * vstack( + sl*max_D, + sl*max_B + ), + max_rhs)); +aux_eqns:factor(ratsimp( + append( + make_eqns(max_D, max_D_values), + make_eqns(max_B, max_B_values) + ))); +aux_eqns:makelist( + lhs(e)*denom(rhs(e))=num(rhs(e)), e, aux_eqns); + +make_time_derivatives(vars, expr):= + block([e], + e:expr, + for v in vars do + e:ratsubst(diff(v,t),%i*omega*v,e), + e); + +aux_eqns:make_time_derivatives( + max_w, + make_time_derivatives(aux_w, expand(aux_eqns)/%i)); +main_eqns:make_time_derivatives(aux_w, main_eqns); + +soln:solve( + append(main_eqns,aux_eqns), + makelist(diff(v,t),v,append(max_w,aux_w))); +print(expand(covect(soln[1]))); + +/* -------------------------------------------------------------------------- */ +/* Compare vs 'nice' shorter expression */ +/* -------------------------------------------------------------------------- */ + +sig:makelist(sigma[c], c, coords); +sigl:shift(sig, -1); +sigr:shift(sig, 1); + +known:append( + max_D/epsilon*(sig - sigl)/epsilon - max_E*sigr/epsilon + + curl(max_H)/epsilon, + + max_B/mu*(sig - sigl)/epsilon - max_H*sigr/epsilon + -curl(max_E)/mu, + + -sigl/epsilon*max_D+curl(max_H), + -sigl/epsilon*max_B-curl(max_E) + ); +/* print(covect(expand(ratsimp(map(rhs, soln[1])-known))));*/ +assert(norm(ratsimp(map(rhs, soln[1])-known))=0); diff --git a/contrib/maxima/maxwell.mac b/contrib/maxima/maxwell.mac new file mode 100644 index 00000000..c336eb20 --- /dev/null +++ b/contrib/maxima/maxwell.mac @@ -0,0 +1,74 @@ +kill(all); + +load("maxwellbase.mac"); + +/* ------------------------------------------------------------------------- */ + +max_radbdryspinw:makelist( + if max_D[i,i] >= 0 then max_sminw[i,1] else 0, +i, 1, 6)$ + +max_wradbdry:fullhypsimp(max_V.max_radbdryspinw); +print("Radiation boundary condition for Maxwell's:"); +print(max_wradbdry); + +max_knownwradbdry:max_wm+1/2*vstack( + crossprod(n, crossprod(n, max_Em)) + - sqrt(mu/epsilon)*crossprod(n, max_Hm), + crossprod(n, crossprod(n, max_Hm)) + + sqrt(epsilon/mu)*crossprod(n, max_Em) + ); +assert(norm_2_squared(hypsimp(max_knownwradbdry - max_wradbdry))=0); + +/* ------------------------------------------------------------------------- */ +max_pecbdrywpins:ev(max_pecbdrywp,[Em=max_Emins, Hm=max_Hmins]); + +print("PEC bdry condition for Maxwell's in terms of characteristic variables (s)"); +max_pecbdryspinw:fullhypsimp(max_invV.max_pecbdrywpins); +print(max_pecbdryspinw); +/* +/* ------------------------------------------------------------------------- */ +/* Simple ("stoopid") upwinded flux */ +/* ------------------------------------------------------------------------- */ + +max_sfluxstoopid:max_D . makelist( + if max_D[i,i] >= 0 then + sm[i,1] + else + sp[i,1], + i, 1, length(max_D)); + +max_wfluxstoopid:fullhypsimp(max_V . ev(max_sfluxstoopid, [sm=max_sminw,sp=max_spinw])); + +print("Stoopid upwind flux for Maxwell's in terms of characteristic variables:"); +print(max_sfluxstoopid); +print("Stoopid upwind flux for Maxwell's in terms of physical variables:"); +print(max_wfluxstoopid); + +/* ------------------------------------------------------------------------- */ +/* Rankine-Hugoniot flux */ +/* ------------------------------------------------------------------------- */ + +print("Maxwell's flux in terms of characteristic variables:"); +print(max_sflux); + +print("Maxwell's flux in terms of physical variables:"); +print(max_wflux); + +/* ------------------------------------------------------------------------- */ +/* known flux value*/ +/* ------------------------------------------------------------------------- */ +max_Z:sqrt(mu/epsilon)$ +max_Y:sqrt(epsilon/mu)$ + +max_knownstrongwflux:ratsimp(vstack( + -1/(2*epsilon) + *(crossprod(n,(max_Hm-max_Hp)-1/max_Z*crossprod(n,max_Em-max_Ep))), + 1/(2*mu) + *(crossprod(n,(max_Em-max_Ep)+1/max_Y*crossprod(n,max_Hm-max_Hp))) + ))$ + +assert(norm_2_squared(hypsimp( + (max_strongwflux) + -max_knownstrongwflux))=0); + */ diff --git a/contrib/maxima/maxwellbase.mac b/contrib/maxima/maxwellbase.mac new file mode 100644 index 00000000..73b68faf --- /dev/null +++ b/contrib/maxima/maxwellbase.mac @@ -0,0 +1,79 @@ +load("myhelpers.mac"); + +assume(c>0); +assume(mu>0); +assume(epsilon>0); +assume(epsinv>0); +assume(muinv>0); + +/* A hyperbolic system matrix resulting from a curl */ +curlmat(coord):=genmatrix( + lambda ([i,j], levi_civita([coord,j,i])), + 3,3); + +max_submat(i):=blockmat( + zeromatrix(3,3), + -epsinv*curlmat(i), /* epsinv = 1/epsilon */ + muinv*curlmat(i), /* muinv = 1/mu */ + zeromatrix(3,3) + ); + +n:[nx,ny,nz]; + +max_Asimp:sum(n[i]*max_submat(i),i,1,3); +max_A:subst([epsinv=1/epsilon,muinv=1/mu], max_Asimp); + +max_invsubst(x):=subst([epsinv=1/epsilon, muinv=1/mu], x); + +[max_V, max_D, max_invV]:max_invsubst(hypdiagonalize(max_Asimp)); + +max_Em:makelist(Em[i],i,1,3); +max_Ep:makelist(Ep[i],i,1,3); +max_Hm:makelist(Hm[i],i,1,3); +max_Hp:makelist(Hp[i],i,1,3); +max_wm:vstack(max_Em,max_Hm); +max_wp:vstack(max_Ep,max_Hp); + +max_sm:makelist(sm[i],i,1,6); +max_sp:makelist(sp[i],i,1,6); + +max_sminw:hypsimp(max_invV.max_wm); +max_spinw:hypsimp(max_invV.max_wp); + +max_wmins:hypsimp(max_V.max_sm); +max_wpins:hypsimp(max_V.max_sp); + +max_Emins:makelist(max_wmins[i,1],i,1,3); +max_Hmins:makelist(max_wmins[i,1],i,4,6); + +/* ------------------------------------------------------------------------- */ +/* Boundary conditions */ +/* ------------------------------------------------------------------------- */ +max_pecbdrywp:append(-max_Em,max_Hm); + +/* ------------------------------------------------------------------------- */ +/* Rankine-Hugoniot flux */ +/* ------------------------------------------------------------------------- */ + +max_Dinc:subst([1/(sqrt(epsilon)*sqrt(mu))=c], max_D); +max_sflux:hyp_upwind_flux([-c,0,c], max_Dinc); + +/* FIXME: max_V should not depend on epsilon and mu, but it does + For now, make cp and cm equal. */ + +max_sflux:subst( + [cp=1/(sqrt(epsilon)*sqrt(mu)), cm=1/(sqrt(epsilon)*sqrt(mu))], + max_sflux); + +max_wflux:fullhypsimp(max_V.ev(max_sflux, [sm=max_sminw,sp=max_spinw])); + +max_stronglocalpart:max_A.max_wm; + +max_strongwflux:max_stronglocalpart-max_wflux; + +/* ------------------------------------------------------------------------- */ +/* vector components */ +/* ------------------------------------------------------------------------- */ + +normalcomp(v):=(v.n)*n; +tangentialcomp(v):=v-normalcomp(v); diff --git a/contrib/maxima/myhelpers.mac b/contrib/maxima/myhelpers.mac new file mode 100644 index 00000000..e61eb12a --- /dev/null +++ b/contrib/maxima/myhelpers.mac @@ -0,0 +1,142 @@ +load("eigen"); +load("itensor"); +load("diag"); + +assert(condition):=if not condition then error("Assertion violated") else true$ + +norm_2_squared(v):=v.v; + +/* A constant matrix of size n x m */ +constmatrix(n,m,c):=genmatrix(lambda ([i,j], c),n,m); + +vstack:append; + +hstack(a,b):=transpose(append(transpose(a),transpose(b))); + +blockmat(a11,a12,a21,a22):=vstack(hstack(a11,a12),hstack(a21,a22))$ + +crossfunc(f):=makelist( + sum(sum( + levi_civita([i,j,k])*f(j,k), + j,1,3),k,1,3),i,1,3)$ + +crossprod(a,b):=crossfunc(lambda([j,k], a[j]*b[k])); + +subrange(v, start, stop):=makelist(v[i,1],i,start,stop); + +extract_diagonal(A):=makelist(A[i][i],i,1,length(A)); +make_diagonal_matrix(v):=genmatrix( + lambda([i,j], if i=j then v[i] else 0), + length(v),length(v)); + +/* ------------------------------------------------------------------------- */ +/* Simplification for expressions stemming from hyperbolic systems */ +/* ------------------------------------------------------------------------- */ + +hypsimp(x):=ratsimp(ratsubst(1,n.n,x))$ + +fullhypsimp(x):=hypsimp( + ratsubst( + last(n)^2, + 1-sum(n[i]^2,i,1,length(n)-1), + x) + ); + +/* ------------------------------------------------------------------------- */ +/* diagonalize a given hyperbolic operator A */ +/* ------------------------------------------------------------------------- */ + +hypdiagonalize(A):=block([evA, V, invV,D], + evA:hypsimp(apply(append, eigenvectors(A)[2])), + V:transpose(apply(matrix, evA)), + invV:hypsimp(invert(V)), + assert(hypsimp(V.invV)=ident(length(A))), + D:hypsimp(invV.A.V), + [V, D, invV]); + +/* ------------------------------------------------------------------------- */ +/* compute upwind flux for a given operator with eigenvalues evs, sorted + * in ascending order. + * Sign assumptions for all variables occuring in evs must be in place. + */ +/* ------------------------------------------------------------------------- */ +hyp_upwind_flux(evs, D):=block([evvars, Dp, Dm, n, midstates, states, unknowns], + evvars:listofvars(evs), + + add_evvars_suffix(suffix, x):=subst(makelist(v=concat(''v, suffix), v, evvars), x), + + evsm:add_evvars_suffix(m, evs), + evsp:add_evvars_suffix(p, evs), + + Dm:add_evvars_suffix(m, D), + Dp:add_evvars_suffix(p, D), + + midstates:makelist(makelist(concat(s,state,i), i, 1, length(D)), + state, 1, length(evs)-1), + + states:append( + [makelist(concat(sm, i), i, 1, length(D))], + midstates, + [makelist(concat(sp,i), i, 1, length(D))]), + + unknowns:apply(append, midstates), + + result:if member(0, evs) then + block([biasedD, veceqns, eqns, soln], + biasedD:makelist( + if evs[i] = 0 then [Dp,Dm] + else if evs[i] > 0 then [Dp,Dp] + else [Dm,Dm], + i, 1, length(evs)), + + veceqns:apply(append, makelist( + -(if evs[i] > 0 then evsp[i] else evsm[i]) *(states[i+1]-states[i]) + +(biasedD[i][1].states[i+1]-biasedD[i][2].states[i]), + i,1,length(evs))), + + eqns:makelist(veceqns[i,1], i, 1, length(veceqns)), + + soln:solve(eqns, unknowns), + assert(length(soln)=1), + + for i: 1 thru length(evs) do + if evs[i] = 0 then return(Dp.subst(soln[1], midstates[i])) + ) + else + block([straddle_idx, Dstates, veceqns, eqns, soln], + straddle_idx:for i: 1 thru length(evs)-1 do + if (evs[i] < 0) and (evs[i+1] > 0) then return(i), + + flux:makelist(concat(flux,i),i,1,length(D)), + + unknowns:append(unknowns, flux), + + Dstates:append( + [Dm.first(states)], + makelist( + if i = straddle_idx then flux + else if evs[i] > 0 then Dp.midstates[i] + else Dm.midstates[i], + i, 1, length(midstates)), + [Dp.last(states)]), + + veceqns:apply(append, makelist( + -(if evs[i] > 0 then evsp[i] else evsm[i]) *(states[i+1]-states[i]) + +(Dstates[i+1]-Dstates[i]), + i,1,length(evs))), + + eqns:makelist(veceqns[i,1], i, 1, length(veceqns)), + + print(covect(eqns)), + soln:solve(eqns, unknowns), + assert(length(soln)=1), + + subst(soln[1], flux) + ), + subst( + append( + makelist(concat(sm, i)=sm[i,1], i, 1, length(D)), + makelist(concat(sp, i)=sp[i,1], i, 1, length(D)) + ), + result) + ); diff --git a/contrib/maxima/pml.mac b/contrib/maxima/pml.mac new file mode 100644 index 00000000..0068fac2 --- /dev/null +++ b/contrib/maxima/pml.mac @@ -0,0 +1,73 @@ +kill(all); +load("maxwellbase.mac"); + +/* +See +S.G. Johnson, “Notes on Perfectly Matched Layers,” 2008; +http://math.mit.edu/~stevenj/18.369/pml.pdf. +*/ + +coords:[x,y,z]; +allvars:append([t],coords); + +max_E:makelist(concat(E,i),i,coords); +max_H:makelist(concat(H,i),i,coords); +max_w:vstack(max_E,max_H); + +depends(max_w,allvars); + +/* +The minus sign takes the operator from conservation form + u_t + A u_x = 0 +into ODE form + u_t = - A u_x. +*/ +sigexpr(c):=(1+%i*sigma[c]/omega); + +pml_op:-subst(makelist(concat('n, i)=1/sigexpr(i)*concat(d, i), i, coords), max_A); +orig_op:-subst(makelist(concat('n, i)=concat(d, i), i, coords), max_A); + +op_rhs(op_mat, state):=block([rhss], + rhss:op_mat.state, + + for c in coords do + for var in state do + rhss:subst(diff(var, c), concat(d, c)*var, rhss), + rhss); + +pml_rhs:op_rhs(pml_op, max_w); +orig_rhs:op_rhs(orig_op, max_w); + +make_ode(my_rhs, state):= + covect( + makelist(-%i*%omega*state[i]=my_rhs[i,1], i, 1, length(state))); + +pml_eqns:make_ode(pml_rhs, max_w); +orig_eqns:make_ode(orig_rhs, max_w); + +kill_denom(eqns):= + covect(makelist( + block([eqn], + eqn:eqns[i,1], + for c in coords do + if not freeof(sigma[c], eqn) then + eqn:eqn*sigexpr(c), + lhs(eqn)=multthru(rhs(eqn))), + i,1,length(eqns))); + +simpler_eqns:kill_denom(pml_eqns); + +/* +auxvar_coeffs: +auxvars:makelist(concat(X,i),i,1,length(auxvar_coeffs)); +depends(auxvars,allvars); + +new_rhss:pml_rhs; +for i: 1 thru length(auxvar_coeffs) do + new_rhss:ratsubst(concat(X,i), %i/omega*auxvar_coeffs[i,1], new_rhss); + +new_nonaux_eqns:covect(makelist( + diff(max_w[i], t) = new_rhss[i,1], + i, 1, length(max_w))); +aux_eqns:covect(makelist(diff(concat(X,i), t)=-auxvar_coeffs[i,1], i, 1, length(auxvar_coeffs))); +*/ diff --git a/contrib/maxima/wave.mac b/contrib/maxima/wave.mac new file mode 100644 index 00000000..cc812f5b --- /dev/null +++ b/contrib/maxima/wave.mac @@ -0,0 +1,99 @@ +kill(all); + +load("myhelpers.mac"); + +/* redefine this to change dimensionality: */ +n:[nx,ny]; + +dims:length(n); + +assume(c>0); + +if dims = 1 then n:[1]; + +esymmatrix(n, v, i,j):=ematrix(n,n,v,i,j)+ematrix(n,n,v,j,i); +wave_A:sum(n[i]*esymmatrix(dims+1, -c, 1+i,1),i,1,dims); + +[wave_V, wave_D, wave_invV]:hypdiagonalize(wave_A); + +wave_vm:makelist(concat(vm,i),i,1,dims); +wave_vp:makelist(concat(vp,i),i,1,dims); +wave_wm:append([um],wave_vm); +wave_wp:append([up],wave_vp); + +wave_sm:makelist(concat(sm,i),i,1,length(wave_D)); +wave_sp:makelist(concat(sp,i),i,1,length(wave_D)); + +wave_sminw:wave_invV.wave_wm; +wave_spinw:wave_invV.wave_wp; + +wave_wmins:wave_V.wave_sm; +wave_wpins:wave_V.wave_sp; + +wave_radbdryspinw:makelist( + if wave_D[i,i] >= 0 then wave_sminw[i,1] else 0, + i, 1, length(wave_D)); +wave_radbdrywp:fullhypsimp(wave_V.wave_radbdryspinw); + +wave_dirbdryspinw:makelist( + if wave_D[i,i] >= 0 then wave_sminw[i,1] else ubc, + i, 1, length(wave_D)); +wave_dirbdrywp:fullhypsimp(wave_V.wave_dirbdryspinw); + +print("Radiation boundary condition for the wave equation:"); +print(wave_radbdrywp); + +print("Dirichlet boundary condition for the wave equation:"); +print(expand(wave_dirbdrywp)); + +wave_known_dirbdrywp:vstack([-n.wave_vm/2 + um/2 + ubc], + -n*um/2 + n*ubc - 1/2*n*(n.wave_vm) + wave_vm); + +assert(norm_2_squared(fullhypsimp(wave_known_dirbdrywp - wave_dirbdrywp))=0); + +print("Hom-dirichlet in characteristic:"); +print(fullhypsimp(wave_invV.vstack( + -columnvector(first(wave_wmins)), + rest(wave_wmins)) + )); + +print("Hom-Neumann in characteristic:"); +print(fullhypsimp(wave_invV.vstack( + columnvector(first(wave_wmins)), + -rest(wave_wmins)) + )); + +/* ------------------------------------------------------------------------- */ + +wave_eigenvalues:makelist(wave_D[i,i], i, 1, length(wave_D)); + +if member(0, wave_eigenvalues) then + wave_sflux:hyp_upwind_flux([-c,0,c], wave_D) +else + wave_sflux:hyp_upwind_flux([-c,c], wave_D); +wave_wflux:ratsimp(wave_V.ev(wave_sflux, [sm=wave_sminw,sp=wave_spinw])); +wave_strongwflux:fullhypsimp(subst([c=cm], wave_A).wave_wm - wave_wflux); + +/* +print("Wave equation flux in terms of characteristic variables:"); +print(wave_sflux); + + +print("Wave equation flux in terms of physical variables:"); +print(wave_wflux); +print("Weak flux divided by (-c), as implemented in StrongWaveOperator:"); +print(hypsimp(ev(wave_wflux, cp=c, cm=c)/(-c))); +*/ + +print("Strong-form wave equation flux in terms of physical variables:"); +/*print(wave_strongwflux);*/ + +print(fullhypsimp(ev(wave_strongwflux))); + +/* Closed-form expression for upwind flux */ +wave_knownstrongwflux:vstack( + [(cp*n.wave_vp - cm*n.wave_vm + cp*up-cm*um)/2], + (cp*n*(n.wave_vp)-cm*n*(n.wave_vm) + cp*n*up - cm*n*um)/2 + ); + +assert(norm_2_squared(fullhypsimp(wave_knownstrongwflux - wave_strongwflux))=0); diff --git a/contrib/notes/dumka3.cpp b/contrib/notes/dumka3.cpp new file mode 100644 index 00000000..1cf78d54 --- /dev/null +++ b/contrib/notes/dumka3.cpp @@ -0,0 +1,393 @@ +#include +#include +#include +#include +using namespace std; +/* VectorFunctionBase - interface with three functions, which have to be implemented + * void compute(Vector& y, double t, Vector& result) - calculates right hand side of the ODE y'=f(y,t) + * double cour(Vector& y, double t) - used if "isComputeEigenValue=false" to + * calculate h_{euler} - step which you would use in explicit Euler method + * Element GetNonZeroElement() - can be used if "isComputeEigenValue=true" in eigenvalue function, if + * iterations of power method have to be restarted, but y,f(y) equals to zero, + * in this case program need to initialize random non-zero vector, because + * vector's Element is unknown in advance + * for example if Vector contains "double"(s) you can return 0.01*(rand()%100) + * but if Vector contains "Point"(s) you can + * return Point(0.01*(rand()%100),0.01*(rand()%100),0.01*(rand()%100)) + */ +template +class VectorFunctionBase { + public: + virtual void compute(Vector& y, double t, Vector& result) = 0; + virtual double cour(Vector& y, double t) = 0; + virtual Element GetNonZeroElement() = 0; +}; + +/* + * dumka3 - solves non-linear mildly - stiff system of ordinary differential equations y'=f(y,t) + * in order to be able to use this function you should implement interface VectorFunctionBase (above) + * t - time + * tend - end of the integration interval + * h0 - initial step-size (this value can be override by the program) + * atol - tolerance of the computation (should be strictly positive) + * rtol - relative tolerance (should be strictly positive) + * f - Class which implement VectorFunctionBase interface, or class , which has 3 functions + * compute(Vector& y, double t, Vector& result) - calculates right hand side + * cour(Vector& y, double t) - used if "isComputeEigenValue=false" to + * calculate h_{euler} - step which you would use in explicit Euler method + * GetNonZeroElement() - can be used if "isComputeEigenValue=true" in eigenvalue function, if + * iterations of power method have to be restarted, but y,f(y) equals to zero, + * in this case program need to initialize random non-zero vector, because + * vector's Element can be complicated like Point etc, user should provide random non-zero element + * y - solution of the ODE y'=f(y,t) + * isComputeEigenValue - if "true": function "eigenvalue" will be used to determine maximum eigenvalue and h_{euler}=2/eigenvalue + * - if "false": user should provide H_{euler} via "cour" function of the VectorFunction f + * + */ +template +void dumka3(double& t, const double& tend, const double& h0,const double& atol, + const double& rtol, VectorFunction& f, Vector& y, bool isComputeEigenValue = false ) { + static const long _n_p=14; + static const long _sz=2400; + static const long _sz2=6; + long numOfRejectedSteps=0; + double minimumValue = 1.e-15; + double err_n=0.; + double h_n=0; + long numOfSteps=0; + double maxCou = 0.; + double meanCou = 0.; + double timeInterval = tend - t; + if (timeIntervaltend" << endl; + return; + }; + if ( atol < minimumValue ) { + cout << "atol must be positive" << endl; + return; + } + if ( rtol < minimumValue ) { + cout << "rtol must be positive" << endl; + return; + } + long numOfRHSEvaluations=0; + long numOfRHSEvaluations0; + double c_2,c_3,c_4,a_21,a_31,a_41,a_32,a_42,a_43,t2,t3,t4,tmp_1,tmp_2; + double r,cou,h_new,t_en,h; + int n_deg[_n_p]={3,6,9,15,21,27,36,48,63,81,135,189,243,324}; + int index_first[_n_p]={1,2,4,7,12,19,28,40,56,77,104,149,212,293}; + int index_last[_n_p]={1,3,6,11,18,27,39,55,76,103,148,211,292,400}; + const double stab_reg[_n_p]={ 2.5005127005e+0, 9.784756428464169e+0, + 23.818760475282560e+0, 68.932817124349140e+0, 136.648186730571300e+0, + 226.8897061613812e+0, 404.7232537379578e+0, 720.7401653373073e+0, + 1337.1268312643200e+0, 3266.0271134242240e+0, 9072.0255133253400e+0, + 17784.544093341440e+0, 29376.4540800858800e+0,50666.52463738415826e+0}; + + const double coef[_sz]; + + + long j; + long _size=y.size(); + Vector z0(_size), z1(_size), z2(_size), z3(_size), oldEigenVector(_size); + double lambdaOld = 0.; + //Vector z0, z1, z2, z3; + if ( h0 < 1.e-14 ) { + cout << "Initial step h0 is too small" << endl; + return; + } + if ( tend <= t ) { + cout << "End-time Tend is less than initial time" << endl; + return; + } + h_new=h0; + // z0=f.compute(y,t); + f.compute(y,t,z0); + if ( isComputeEigenValue ) { + f.compute(z0,t,oldEigenVector); + numOfRHSEvaluations++; + } + numOfRHSEvaluations++; + long stepId=0; + bool stepOK = true; + do { + if ( isComputeEigenValue ) {//calculate eigen value only every 20th step + if ( stepId%20 == 0 || !stepOK ) { + double ev = abs( eigenvalue(y, z0, t, f, oldEigenVector, z1, lambdaOld, + numOfRHSEvaluations,minimumValue) ); + if ( ev < minimumValue ) { + ev = minimumValue; + } + cou = 2./ (ev*1.2); + } + } else { + cou = f.cour (y, t); + } + if ( maxCou < cou ) { + maxCou=cou; + }; + meanCou += cou; + numOfSteps++; + h=min(h_new,(tend-t)); + // find degree of the polynomials to be used + int index=0; + while( (index!=_n_p-1) && ((2.e+0*h/cou)>stab_reg[index]) ) { + index++; + }; + + if(index>0) { + if(((stab_reg[index-1]/((double) n_deg[index-1]))*((double) n_deg[index]))>(2.*h/cou)) { + index=index-1; + }; + }; + int n_pol_degree=n_deg[index]; + cout << "Time= " << t << " Degree of the polynomial=" << n_pol_degree << endl; + + h=min((stab_reg[index]*cou/2.e+0),h); + + long n_pol=0; + numOfRHSEvaluations0 = numOfRHSEvaluations; + for (int k=index_first[index]-1; k < index_last[index]; k++ ) { + // save initial conditions (solution, right hand side, time) + // that will be needed if time step will be rejected + if(k==index_first[index]-1) { + for ( j=0; j<_size; j++ ) { + z2[j]=z0[j]; //z2=z_en + z3[j]=y[j]; //z3=y_en + }; + t_en=t; + }; + n_pol=n_pol+3; + int _idx=_sz2*k; + a_21=h*coef[_idx]; + c_2=a_21; + a_31=h*coef[_idx+1]; + a_32=h*coef[_idx+2]; + c_3=a_31+a_32; + a_41=h*coef[_idx+3]; + a_42=h*coef[_idx+4]; + a_43=h*coef[_idx+5]; + c_4=a_41+a_42+a_43; + //long _size=y.size(); + for ( j=0; j<_size; j++ ) { + y[j]=y[j]+a_21*z0[j]; + }; + t2=t+c_2; + + // marker W ****************************** + // z1=f.compute(y,t2); + f.compute(y,t2,z1); + numOfRHSEvaluations++; + if(n_pol==n_deg[index]) { + r=h*(coef[_idx+1]-coef[_idx]); + for ( j=0; j<_size; j++ ) { + y[j]=y[j]+r*z0[j]+a_32*z1[j]; + }; + } else { + for ( j=0; j<_size; j++ ) { + y[j]=y[j]+a_32*z1[j]; + }; + }; + + // marker X ****************************** + t3=t+c_3; + if(n_pol==n_deg[index]) { + tmp_1=c_3/2.e+0; + tmp_2=(c_4-c_2)/2.e+0; + for ( j=0; j<_size; j++ ) { + z1[j]=tmp_1*z1[j]-tmp_2*z0[j]; + }; + }; + + // marker Y ****************************** + // z0 = f.compute(y,t3); + f.compute(y,t3,z0); + numOfRHSEvaluations++; + if(n_pol==n_deg[index]) { + for ( j=0; j<_size; j++ ) { + z1[j]=tmp_2*z0[j]+z1[j]; + }; + }; + // marker Z ****************************** + for ( j=0; j<_size; j++ ) { + y[j]=y[j]+a_43*z0[j]; + }; + + t4=t+c_4; + t=t4; + // z0=f.compute(y,t); + f.compute(y,t,z0); + numOfRHSEvaluations++; + }; + + //long _size=z0.size(); + for ( j=0; j < _size; j++ ) { + z1[j]=(1./(rtol*max(sqrt(y[j]*y[j]),sqrt(z3[j]*z3[j]))+atol))*((-tmp_1)*z0[j]+z1[j]); + }; + h_new=h; + // check error and return new recommended step size in variable h_new + + if(!isStepAccepted(t,y,z1,h_new,n_pol_degree,err_n, h_n, stepId)) { + numOfRejectedSteps = numOfRejectedSteps + (numOfRHSEvaluations - numOfRHSEvaluations0); + stepOK = false; + f.isStepAccepted = false; + //long _size=y.size(); + /* code for paraller exceution via OpemMP - you may want to use it +#pragma parallel shared(y,z0,z2,z3) local(j) +#pragma pfor +*/ + for ( j=0; j < _size; j++ ) { + z0[j]=z2[j]; + y[j]=z3[j]; + }; + t=t_en; + cout << "Step is rejected. Used degree of polynomial: " << n_pol_degree << " Time= " << t << endl; + } else { + stepOK = true; + f.isStepAccepted = true; + } + stepId++; + } while ( t < tend ); + meanCou = meanCou/((double)numOfSteps); + cout << "Number of RHS evaluations = " << numOfRHSEvaluations << endl; + cout << "Mean step-size = " << timeInterval/((double)numOfRHSEvaluations) << endl; + cout << "Mean value of cou = " << meanCou << endl; + // cout << "Maximum value of cou = " << maxCou<< endl; + cout << "Last value of cou = " << cou << endl; + cout << "Number of rejected RHS evaluations = " << numOfRejectedSteps << endl; + //delete [] z; +}; + +/* norm1 - function isStepAccepted can use different norms, but we implemente "max=max|u_i|, i= 1..n" norm + * one can optimized this function by replacing qrt(z[i]*z[i]) to abs(z[i]), but you must define + * function abs for Element of the Vector + */ +template +double norm1 (const Vector& z) { + double eps1=0.e+0; + long _size=z.size(); + for ( int i=0; i < _size; i++ ) { + double absValue = sqrt(z[i]*z[i]); + if(eps1 +bool isStepAccepted(const double &t, const Vector &y, const Vector &z, + double &h_new, const int &n_pol_degree, double err_n, double h_n, long stepId) { + + bool isnotreject = false; + double eps1=0.e0; + long size = y.size(); + for ( int i=0;i < size; i++ ) { + eps1=eps1+z[i]*z[i]; + } + double eps=sqrt( eps1/((double)size) ); + double fracmin=0.1e0; + if(eps==0.e0) { + eps=1.e-14; + } + double frac=pow(1.e0/eps,1.e0/3.e0); + if(eps <= 1.e0) { + if (err_n>0.e0 && h_n > 0.e0) { + double frac2=pow(err_n,(1.e0/3.e0))*frac*frac*(h_new/h_n); + frac=min(frac,frac2); + } + isnotreject = true; + double fracmax=2.e0; + frac=min(fracmax,max(fracmin,0.8e0*frac)); + double h_old = h_new; + h_new=frac*h_new; + h_n=h_old; + err_n=eps; + } else { + isnotreject = false; + double fracmax=1.e0; + frac=0.8e0*min(fracmax,max(fracmin,0.8e0*frac)); + if(stepId==0) { + h_new=fracmin*h_new; + } else { + h_new=frac*h_new; + } + cout << eps << t << n_pol_degree << h_new << h_n << endl; + } + return isnotreject; +}; + +/* norm - calculates norm of the Vector, this function works only if Element of + * the Vector has operation multiplication "*" + */ +template +double norm ( const Vector &x ) { + int size = x.size(); + double norm = 0.; + for ( int i=0; i < size; i++ ) { + norm += x[i]*x[i]; + } + norm = sqrt(norm); + return norm; +} + +/* eigenvalue - calculates eigen values of the right hand side + * The used algorithm is a slight modification of the algorithm + * used in ROCK and RKC + */ + template +double eigenvalue(const Vector& y,Vector& fy, double t, + VectorFunction& f, Vector& v, Vector& result, + double& lambdaOld, long& numOfRHSEvaluations, double minimumValue) +{ + cout << "Compute eigenvalue ..." << endl; + int maxit = 30; + double tol = 0.01; + int idx = 0; + int size = y.size(); + double radius = 1.e-10; + double lambda = lambdaOld; + double nrmY = norm(y); + double nrmV = norm(v); + //f.compute(result,t,v); + //f.compute(fy,t,v); + do { + idx++; + // all vectors should be ||y-v|| < radius, because problem can be non-linear + if ( nrmV < minimumValue ) + { + // change initial vector and restart - if restarted check that + // initial vector has not been used already + for ( int i=0; i < size; i++ ) + { + v[i] = f.GetNonZeroElement(); + }; + nrmV = norm(v); + if ( nrmV < minimumValue ) { + cout << "Error occur: your GetNonZeroElement function should return a value with strictly positive norm, but it returns 0!" << endl; + cout << "Fix GetNonZeroElement function, please." << endl; + return lambdaOld; + } + } + double norminv= radius/nrmV; + for ( int i=0; i < size; i++ ) { + v[i] = y[i] + norminv * v[i]; + } + + f.compute(v,t,result); + numOfRHSEvaluations++; + // approximation of eigenvector + for ( int i=0; i < size; i++ ) { + v[i] = result[i] - fy[i]; + } + // approximation of eigenvalue + nrmV = norm(v); + lambdaOld = lambda; + lambda = nrmV/radius; + cout << "Iteration = " << idx << " eigenvalue = " << lambda << endl; + } + while ( abs(lambda - lambdaOld) > abs(lambda)*tol && idx <= maxit ); + lambdaOld = lambda; + return lambda; +} diff --git a/contrib/notes/em-cheat.tm b/contrib/notes/em-cheat.tm new file mode 100644 index 00000000..337fe0bb --- /dev/null +++ b/contrib/notes/em-cheat.tm @@ -0,0 +1,50 @@ + + + + +<\body> + > + + 1. Units: + + <\eqnarray*> + ||>>|||=s>>>|]>||>>>|]>||A|mN>>>||||As>>>|||s|N*m>>>>> + + + 2. Equations + + <\eqnarray*> + >>||\H+J>>|>>||\E>>|\D>||>>|\B>||>>> + + + <\equation*> + B=\H,D=\E + + + <\equation*> + Z=|\>>,Y=|\>>. + + + 3. Boundary conditions + + <\equation*> + n\E=0,n\H=0. + + + 4. Cross product, curl + + <\equation*> + >>|>>|>>>>>\>>|>>|>>>>>=b-ab>>|b-ab>>|b-ab>>>>>,curl + F=F-\F>>|F-\F>>|F-\F>>>>> + + + <\equation*> + a\(b\c)=b(a\c)-c(a\b) + + + +<\initial> + <\collection> + + + \ No newline at end of file diff --git a/contrib/notes/fluxfan.fig b/contrib/notes/fluxfan.fig new file mode 100644 index 00000000..6672b2ff --- /dev/null +++ b/contrib/notes/fluxfan.fig @@ -0,0 +1,40 @@ +#FIG 3.2 Produced by xfig version 3.2.5 +Landscape +Center +Metric +A4 +100.00 +Single +-2 +1200 2 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 1 1 1.00 60.00 120.00 + 2160 7290 8640 7290 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 1 1 1.00 60.00 120.00 + 5220 7290 5220 3960 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 5220 7290 7920 6120 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 5220 7290 6210 5220 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 5220 7290 4410 5130 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 5220 7290 3060 5580 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 5220 7290 3690 4950 +2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2 + 5220 7290 7110 5220 +4 1 0 50 -1 0 15 0.0000 6 195 360 8820 7380 $x$\001 +4 1 0 50 -1 0 15 0.0000 6 195 315 5220 3780 $t$\001 +4 1 0 50 -1 0 15 0.0000 6 210 1275 8010 6120 $\\lambda_n$\001 +4 1 0 50 -1 0 15 0.0000 6 210 1275 4410 5040 $\\lambda_k$\001 +4 1 0 50 -1 0 15 0.0000 6 195 585 6840 7020 $s^+$\001 +4 1 0 50 -1 0 15 0.0000 6 195 525 3420 7020 $s^-$\001 +4 1 0 50 -1 0 15 0.0000 6 210 1275 2880 5490 $\\lambda_1$\001 +4 1 0 50 -1 0 15 0.0000 6 225 1080 5580 5670 $s^{*(k)}$\001 +4 1 0 50 -1 0 15 0.0000 6 195 825 4320 5580 $\\cdots$\001 +4 1 0 50 -1 0 15 0.0000 6 195 825 6390 5760 $\\cdots$\001 +4 1 0 50 -1 0 15 0.0000 6 225 1770 6120 5040 $\\lambda_{k+1}$\001 +4 1 0 50 -1 0 15 0.0000 6 225 1275 6750 6210 $s^{*(n-1)}$\001 +4 1 0 50 -1 0 15 0.0000 6 195 570 4140 6210 $s^*$\001 diff --git a/contrib/notes/hedge-notes.tm b/contrib/notes/hedge-notes.tm new file mode 100644 index 00000000..5b714be9 --- /dev/null +++ b/contrib/notes/hedge-notes.tm @@ -0,0 +1,1530 @@ + + +> + +<\body> + + + + <\input|> >> + )clear all + + + <\output> + \ \ \ All user variables and function definitions have been cleared. + + + <\input|> >> + verteq:=[2/sqrt(3)*vector [sin(i*2*%pi/3),cos(i*2*%pi/3)] for i in + 0..2] + + + <\output> + 0,|3>,1,-|3>,-1,-|3>(1)> + + + + + <\input|> >> + vertun:= [vector[-1,1], vector[1,-1], vector[-1,-1]] + + + <\output> + -1,1,1,-1,-1,-1(2)> + + + + + <\input|> >> + A:=matrix [[a,b],[c,d]] + + + <\output> + |||>||>>>>(3)> + + + + + <\input|> >> + v:=vector [e,f] + + + <\output> + e,f(4)> + + + + + <\input|> >> + eqns:=concat( + + [(A*verteq.i+v).1=vertun.i.1 for i in 1..3], + + [(A*verteq.i+v).2=vertun.i.2 for i in 1..3]) + + + <\output> + +3e|3>=-1,+3e+3a|3>=1,+3e-3a|3>=-1,+3f|3>=1,+3f+3c|3>=-1,+3f-3c|3>=-1(5)> + + + + + <\input|> >> + solve(eqns,[a,b,c,d,e,f]) + + + <\output> + a=1,b=->,c=0,d=>,e=-,f=-(6)> + + + + + <\input|> >> + \; + + > + + + + + <\input|> >> + )clear all + + + <\output> + \ \ \ All user variables and function definitions have been cleared. + + + <\input|> >> + verteq:=[ + + vector [-1,-1/sqrt(3),-1/sqrt(6)], + + vector [ 1,-1/sqrt(3),-1/sqrt(6)], + + vector [ 0, 2/sqrt(3),-1/sqrt(6)], + + vector [ 0, \ \ \ \ \ \ \ \ 0, 3/sqrt(6)] + + ] + + + <\output> + -1,-|3>,-|6>,1,-|3>,-|6>,0,|3>,-|6>,0,0,|2>(1)> + + + + + <\input|> >> + vertun:= [ + + [-1,-1,-1], + + [ 1,-1,-1], + + [-1, 1,-1], + + [-1,-1, 1] + + ] + + + <\output> + -1,-1,-1,1,-1,-1,-1,1,-1,-1,-1,1(2)> + + + + + <\input|> >> + A:=matrix [[a,b,c],[d,e,f],[g,h,i]] + + + <\output> + |||||>|||>|||>>>>(3)> + + + + + <\input|> >> + v:=vector [j,k,l] + + + <\output> + j,k,l(4)> + + + + + <\input|> >> + eqns:=concat( + + [(A*verteq.i+v).1=vertun.i.1 for i in 1..4], + + concat( + + [(A*verteq.i+v).2=vertun.i.2 for i in 1..4], + + [(A*verteq.i+v).3=vertun.i.3 for i in 1..4] + + )) + + + <\output> + j-|6>c-|3>b-a=-1,j-|6>c-|3>b+a=1,j-|6>c+|3>b=-1,j+|2>c=-1,k-|6>f-|3>e-d=-1,k-|6>f-|3>e+d=-1,k-|6>f+|3>e=1,k+|2>f=-1,l-|6>i-|3>h-g=-1,l-|6>i-|3>h+g=-1,l-|6>i+|3>h=-1,l+|2>i=1(5)> + + + + + <\input|> >> + solve(eqns,[a,b,c,d,e,f,g,h,i,j,k,l]) + + + <\output> + a=1,b=-|3>,c=-|6>,d=0,e=|3>,f=-|6>,g=0,h=0,i=|2>,j=-,k=-,l=-(6)> + + + + + <\input|> >> + \; + + > + + + + + <\input|> >> + )clear all + + + <\output> + \ \ \ All user variables and function definitions have been cleared. + + + <\input|> >> + f:=operator 'f + + + <\output> + (1)> + + + + + <\input|> >> + g:=operator 'g + + + <\output> + (2)> + + + + + <\input|> >> + a:=2*(1+r)/(1-s)-1 + + + <\output> + (3)> + + + + + <\input|> >> + D(a,r) + + + <\output> + (4)> + + + + + <\input|> >> + D(a,s) + + + <\output> + -2s+1>(5)> + + + + + <\input|> >> + expr:=sqrt(2)*f(a)*g(s)*(1-s)^i + + + <\output> + fgs-s+1(6)> + + + + + <\input|> >> + D(expr,r) + + + <\output> + gs-s+1f|s-1>(7)> + + + + + <\input|> >> + D(expr,s) + + * + + (s-1)^2 + + + <\output> + s-2s+1f-s+1gs+2r+2gs-s+1f+-i*s+2i*s-ifgs-s+1i-1>(9)> + + + + + <\input|> >> + \; + + > + + <\eqnarray*> + ||>s-2s+1f-s+1gs>>|||>2r+2gs-s+1f>>|||>-i*s+2i*s-ifgs-s+1i-1>>>|||fa1-sgs>>|||2r+2gs1-sfa>>|||ags1-si-1>>>>> + + + + + + <\input|> >> + )clear all + + + <\output> + \ \ \ All user variables and function definitions have been cleared. + + + <\input|> >> + f:=operator 'f + + + <\output> + (1)> + + + + + <\input|> >> + g:=operator 'g + + + <\output> + (2)> + + + + + <\input|> >> + h:=operator 'h + + + <\output> + (3)> + + + + + <\input|> >> + a:= -2*(1+r)/(s+t) - 1 + + + <\output> + (4)> + + + + + <\input|> >> + b:=2*(1+s)/(1-t) - 1 + + + <\output> + (5)> + + + + + <\input|> >> + expr:=sqrt(8)*f(a)*g(b)*(1-b)^i*h(c)*(1-c)**(i+j) + + + <\output> + fghc-c+1j+i>(6)> + + + + + <\input|> >> + D(expr,r) + + + <\output> + ghc-c+1j+i>f|t+s>(7)> + + + + + <\input|> >> + dfds:=D(expr,s); + + + <\output> + + + + <\input|> >> + dfdsden:=denominator dfds + + + <\output> + +2s-1t+s-2st-s(20)> + + + + + <\input|> >> + dfds*dfdsden + + + <\output> + -4t-8s*t-4sfhc-c+1j+i>g+4r+4t-4r-4ghc-c+1j+i>f+4i*t+8i*s*t+4i*sfghc-c+1j+i>i-1>(18)> + + + + + <\input|> >> + dfdt:=D(expr,t); + + + <\output> + + + + <\input|> >> + dfdtden:=denominator dfdt + + + <\output> + +2s-2t+s-4s+1t+-2s+2st+s(23)> + + + + + <\input|> >> + dfdt*dfdtden + + + <\output> + 4s+4t+8s+8st+4s+4sfhc-c+1j+i>g+4r+4t+-8r-8t+4r+4ghc-c+1j+i>f+-4i*s-4it+-8i*s-8i*st-4i*s-4i*sfghc-c+1j+i>i-1>(24)> + + + + + <\input|> >> + \; + + > + + + + + + Let > be the th Lagrange + interpolation polynomial and =(u)> be the + vector of nodal coefficients of our approximated function + (x)=ul(x)>. ( + without subscript denotes a continuous function .) + + \T> represents the discrete + faces of the element >, and + > is the Lagrange interpolation polynomial + on the face corresponding to the node with global number + , or constant zero if there is no such function (such as + when is not the number of a node in ). + Further, + + <\eqnarray*> + >|>|>ll,>>|x>>>|>|>l\>>l,>>|>|>|ll.>>>> + + + Note that =(M)> and + =(M)>. + + Recall |^>=\> where + p=ul>, and + + <\equation*> + V\(x)=(x)>||)>>|>||>>|(x)>|>|(x)>>>>>\(x)=(x)p(x)+\+l(x)p(x)>>|>>|(x)p(x)+\+l(x)p(x)>>>>>=(x)>>|>>|(x)>>>>>, + + + simply because the sum of Lagrange polynomials uniquely interpolates + (x)>. We thus find + + <\eqnarray*> + >||>ll=J([V]p)([V]p)>>|||[V][V]pp>>|||[V][V]\>>|||[V][V]>>|||[V][V]>>|||[(V*V)].>>>> + + + To obtain the stiffness matrix x>>>, realize that we really only need a + differentiation matrix >>> that maps + the nodal values of a function to those of its derivative. If we have + >>>, we can simply recycle the mass + matrix: + + <\equation*> + Sx>>=M*D>>. + + + Before finding the global differentiation matrix + >>>, we first endeavor to find the + local differentiation matrix >>>, where + we denote local (unit) coordinates and global + coordinates . We define the derivative Vandermonde + matrices r>>\[\>p(r)]>. Recalling |^>=\>, we obtain + r>>|^>=Vr>>V\=D>>\> + and thus >>=Vr>>V>. + For clarity, and again using the summation convention, consider the + identity + + <\equation*> + (x)=ul(x)=u(x)|\>>l(x). + + + If we define a function :r\x>, then the + global differentiation matrix is then given by + + <\equation*> + [D>>]=\>>(l(F(x)))=(\l)(r)(F)(x),\>=[D>>](F)(x),\>. + + + + + + + We're discretizing +v\\u=0>. The + analytic solution is + + <\equation*> + u(x,t)=u(x-v*t). + + + Recall + + <\eqnarray*> + \(v*u*\)>||\(v*u)\+v*u\\\=(v\\*u)\+v*u\\\.>>>> + + + Using the summation convention, we find + + <\eqnarray*> + ||>u\+>(v\\u)\>>|||>u\->v*u\\\+>\\(v*u*\)>>|||>u\->v*u\\\+T>v*u\\n>>||>|>u\->v*u\\\+T>(v*u)>\\n>>||>|>\ull->ul>>|ul>>>>>\>l>>|>l>>>>>+\T>u)>ll>>|u)>ll>>>>>\n>>|||>|||>>> + + + + + <\eqnarray*> + ||>u\->a*u\\\+T>(v*u)>\\n>>|||>u\+>(v\\u)\-T>(v*u-{v*u})\\n>>|||>u\+>(v\\u)\-T>v>>|>>>>>\>>|>>>>>\\n>>>> + + + + + + + \; + + We have -\u=0>, so + + <\eqnarray*> + -\\v>||>|-\u>||>>> + + + Note that is a vector. \ Then, using the summation + convention, + + <\eqnarray*> + ||>u\->(\\v)\>>|||>u\+>v\\\->\\(v\)>>|||>u\+>v\\\-T>v\\n>>||>|>u\+>v\\\-T>v>\\n>>||>|>\ull+>]l>>|]l>>>>>\>l>>|>l>>>>>-\T>]>ll>>|]>ll>>>>>\n>>|||\u+Sx>[v]+Sx>[v]-\T>[v]>Mn+[v]>Mn>>|||)\\+(Sx>)\+(Sx>)\-\T>[v]>(M)n+[v]>(M)n>>>> + + + We set >\{\}>. Likewise, + we get + + <\eqnarray*> + ||>v\\->\u\\>>|||>v\\+>u\\\->\\(u>)>>|||>v\\+>u\\\-T>u\\n>>||>|>v\\+>u\\\-T>u>\\n>>||>|>[v]l>>|[v]l>>>>>\>>|>>>>+>ul\\>>|>>>>-\T>>ll>>|>l0>>>>>\n>>>|>[v]l>>|[v]l>>>>>\>|>>>>>+>ul\\>|>>>>>-\T>>l0>>|>ll>>>>>\n>>>>>>>>|||>\[v]ll+>ul\>l-\T>u>lln>>>|>\[v]ll+>ul\>l-\T>u>lln>>>>>>>>|||[v]M+uSx>-\T>u>Mn>>>|[v]M+uSx>-\T>u>Mn>>>>>>>>>> + + + and we set >\{u}>. + + + + First equation: + + <\eqnarray*> + ||>u\+>v\\\-T>v>\\n>>|||>u\->(\\v)\+T>(v-v>)\\n>>|||>u\->(\\v)\+T>[v]\\n>>>> + + + Second equation: + + <\eqnarray*> + ||>v\\+>u\\\-T>u>\\n>>|||>v\\->(\u)\+T>(u-u>)\\n>>|||>v\\->(\u)\+T>[u]\\n>>>> + + + + + Begin with + + <\eqnarray*> + -\\(q)>||>|\u>||>>> + + + so + + <\eqnarray*> + [u-\\(q)]\>||>|u\+*q\\\-\\(*q*\)>||>|u\+*q\\\-T>*q*\\n>||>|u\+*q\\\-T>(*q)>\\n>||>>> + + + and + + <\eqnarray*> + [q-\u]\\>||>|q\\-\\(\u)>||>|q\\+\\(*u)\-\\(*u\)>||>|q\\+\\(*u)\-T>n\(*u)>\>||>>> + + + + + Golub-Welsch recursion: + + <\equation*> + p=(\x+\)p-\p + + + Hesthaven-Warburton ``recursion'': + + <\equation*> + x*p=ap+bp+ap + + + Solve for > and > from + G-W: + + <\eqnarray*> + >||x*p+\p-\p>>|x*p>||-\p+\p>>|>||>p-|\>p+|\>p>>||||\>|\>>p|\>|\>>p+>|\>>p.>>>> + + + + + + + Start with a linear hyperbolic system: + + <\equation*> + \+A\\=0 + + + <\itemize> + Pick a unit normal > of an interface where the states + > and > meet. + + Compute + + <\equation*> + A>(\)\nA> + + + and diagonalize it to >=V*>A>(\)*(V>)>, + resulting in eigenvalues + + <\equation*> + \>\\\\>\0\\>\\\\>. + + + Call >\V>\>>. + + Consider the space-time diagram of Figure . + Each eigenvalue is imagined to be the speed of a emanating + from the point under consideration, and each shock separating one + from another. The left- and rightmost state are + > and >, of course. The next + state from the left is labeled >>, the one + after that \>>, up to + (n-1)>>. These states are a-priori unknown. In + addition, the fluxes )>,\,(D\)(n-1)>> + are also unknown: For the moment, there are unknowns. + + The Rankine-Hugoniot condition + + <\equation*> + \>=)(i)>-(D\)(i-1)>|\(i)>-\(i-1)>> + + + has to hold across each such interface. For convenience, set + (0)>=\> and + (n)>=\>. This will provide + equations. + + Naturally, eigenvalues larger than zero (``right-travelling'') use their + > values, whereas left-travelling values smaller + than zero use their > values. + + We label > that eigenvalue whose state + straddles zero the line of zero speed--that includes the case of + =0>. If there are several zero eigenvalues (can + that even happen?), pick as large as possible. It is assumed + that the number is the same on both sides of the interface, and + therefore it carries no > superscript. + + The numerical flux, i.e. the solution that we seek, is the flux + )(k)>> that straddles zero. + + For non-zero-straddling states, the unknowns + )(i)>> are found as follows: + + <\itemize> + If k>, )(i)>=D\(i)>>. + + If k>, )(i)>=D\(i)>>. + + + This takes unknowns out of the picture, leaving us with + unknowns. + + The case of a zero eigenvalue deserves special mention, even though + its handling is technically no different than before: The states + (k-1)>> and (k)>> + fall cleanly to the left and right of zero, so we can use + + <\eqnarray*> + )(k-1)>>||\(k-1)>,>>|)(k)>>||\(k)>.>>>> + + + Observe that we end up with only unknowns (namely all the + (i)>}>), but we still have + equations. It seems, however, that the resulting system of + equations has rank , and so a unique solution can be + determined. () + + Lastly, observe that the th Rankine-Hugoniot condition dictates + + <\equation*> + (D\)(k-1)>=(D\)(k)>, + + + ensuring consistency across the zero-speed shock. + + + I believe that the above can also be applied in the case of degenerate + eigenvalues, but haven't checked in detail. + + |Space-Time + Diagram of a flux fan.> + + + + + + I'm not entirely sure that the case of zero eigenvalues actually comes out + well. Let's consider the following: adopting your notation above, let's + assume we have eigenvalues for each element, + >>, n>. Let + , , and represent the number of negative, zero, + and positive eigenvalues, respectively. (Recall that we're assuming that + ) = sign(\)> + for all .) We have , and we assume at least two of + these numbers are nonzero. (If only one of them is nonzero, upwinding is + trivial.) + + \; + + We start with the unknowns + (i)>, + D*(i)>>.\ + + \; + + I'm adopting the convention that positive eigenvalues correspond to + right-traveling waves. We have the following Rankine-Hugoniot conditions: + + <\equation> + |>||)(i)>-(D\)(i-1)>|\(i)>-\(i-1)>>,>||,N>>|||||>|||)(i)>-(D\)(i-1)>,>||,N+Z>>|||||>|>||)(i)>-(D\)(i-1)>|\(i)>-\(i-1)>>,>||,n>>>>> \ n \ (linear) + equations + + + \; + + \; + + If I understand what we're doing, we consider the unknowns + D*(i)>, + , \, n>> as anisotropic + (directionally-biased) unknowns with explicit connections to the + (i)>> given by the + auxilliary equations + + <\equation> + )(i)> = + D\(i)>>||,N-1>>|||>|)(i)> + = D\(i)>>||, + n>>>>> \ N+P-2 \ \ (linear) equations + + + \; + + However, I don't see what we do for the unknowns which make contact with + the zero eigenvalue(s). You mention that if with + =0> we would have\ + + <\equation*> + )(k)>>|>>|\(k)>.>>>>> + + + But then I'd also claim that\ + + <\equation*> + )(k-1)>>|>>|\(k-1)>.>>>>> + + + I cannot think of a reason why we'd treat + )(k-1)>> any differently from + )(k)>> since I don't think it's the + case that )(k)>> `belongs' to + =0> in any way. If , then indeed + () gives us equations to add to the + equations in () giving equations + to couple the unknowns. But if , then + () gives us equations to couple with + the from (). That's equations + for unknowns. Conversely, if instead we adopted the convention + + <\equation*> + )(k)>>|>>|\(k)>>>|||>|)(k-1)>>|>>|\(k)>,>>>>> + + + then we'd have 2 equations plus from + () plus from () giving + us equations for unknowns.\ + + \; + + In general, if we don't adopt any directional bias for states adjacent to + the zero eigenvalue, then for 0> zero eigenvalues, we'd have + equations for unknowns, and if we do adopt + directional bias for states adjacent to the zero eigenvalue, then we'd have + equations for unknowns (If + =\=0>, I don't see how you can + adopt a directional bias for )(k)>>.) + + \; + + You make the comment above that ``We simply regard + )(k)>> as another unkown of the system, which, + in addition to (i)>}>, brings + us to unknowns.'' + + \; + + If I understand this correctly, this indeed does make equations + and unknowns, but you treat )(k)>> + differently than )(k+1)>> (non + directionally-biased vs. directionally biased) and, again, I don't see how + you can consistently make this choice given that, to me, + )(k+1)>> is no more directionally + unbiased than )(k)>>.\ + + + + \; + + Another comment: degenerate eigenvalues should be ok: let's suppose that + =\>. Then there is an + intermediate state (k)>> + separating (k-1)>> + and (k+1)>>. This + will give us\ + + <\equation*> + =)(k+1)>-(D\)(k)>|\(k+1)>-\(k)>>>>|>|=)(k)>-(D\)(k-1)>|\(k)>-\(k-1)>>>>>>> + + + Using these equations to eliminate the unknowns + )(k)>> and + (k)>> gives us the Rankine-Hugoniot + condition relating (k-1)>> to + (k+1)>>: + + <\equation*> + \=)(k+1)>-(D\)(k-1)>|\(k+1)>-\(k-1)>>. + + + This is the same condition one would get if we pretended that states + )(k)>> and + (k)>>, the states straddled by the + shocks of the same speed, did not exist and merrily went about our + upwinding.\ + + \; + + Note that this also means that a -fold degeneracy in the zero + eigenvalue isn't a problem either, as long as we can sort out the simple + zero eigenvalue case. + + + + We seem to be pretty confident that if there does not exist any 0 + eigenvalue, then everything is fine. Kittens meow, rainbows shine, and the + world is happy. But we've also established that we don't know what's going + on when \ \ \=0>. However, let's boil things down + to the simplest case we can. We know degenerate eigenvalues pose no problem + (see my comment last time), and multiple simple (+) or () + eigenvalues are likewise fine. So let's consider this model case: + + <\equation*> + A>(|^>):\, + \ \ \>\0, + \>=0, \ \ \ and + \ \ \ \>\0. + + + We assume the greatest generality possible here: + >>> can be different magnitudes, + but they have the same sign. In line with a strong hyperbolic system, we + assume that >> is diagonalizable as + + <\equation*> + A>=V>\>V>, + + + where >> are diagonal matrices containing the + >>. We introduce the eigenvectors + + <\equation*> + (V>) =>>|>>|>>>>>> + + + Here, all eigenvalues are simple, and >> + span > for each (+) or (). + there is no reason to assume that, e.g. + , + , + > span + >. Or we shall not assume so here, at + least. Although, to be fair, if this is not the case, then we can have the + case of a genuine shock forming, which is kind of ridiculous. Of course, we + have that >> is associated with + >>. We also assume + that\ + + <\equation*> + =\>||=\>>>>>, + + + which we can do since the >'s + span. And of course, I just realized that the + >> are the entries in + >>. Go math.\ + + \; + + We introduce the intermediate states >> + and \>>, + with their appropriate fluxes )>> + and )\>>. + All four of these, at the moment, are unknown.\ + + \; + + We note that at the end of the day, we want some flux\ + + <\equation*> + (A) + + + (yes that's horrible notation; leave me alone! wahhh!) + + \; + + All of our algorithms will attempt to find this quantity. Now, we consider + three possible algorithms for upwind flux computations.\ + + + + The motivation here is simple: the 255 view of upwinding is the following: + the evolution equation for > should not allow boundary + conditions that impose >; + similarly, the evolution equation for > should not allow + boundary conditions that impose >. + This is algorithmically consistent with vanilla upwinding + =a*u>. The algorithm then is the following: + + <\enumerate> + We simply set >=-\*-\=\>. + Also set \>=-\*-\=\>.\ + + Set )>=A>> + and )\>=A*\>>.\ + + Use the flux )=(A)>+(A*)\>=\+\>. + + + \; + + Advantages of this flux: + + <\itemize> + It's easy + + It's very intuitive + + It's consistent with usual upwinding if =A>. + + + Disadvantages: + + <\itemize> + There is no motivation for it other than, ``it's upwinding''. + + Where'r the R-H conditions? Actually, these are + exactly those conditions: see below.... + + + + + In this formulation, we impose the Rankine-Hugoniot conditions, which in + this case read + + <\equation*> + |(>-)>||A>-A>>|||>|>>||\>>>|||>|(\>-)>||A\>-A>>>>> + + + \; + + However, we also impose `directional bias', which means that intermediate + states which are adjacent to, but not overlapping, the zero-speed shock + inherit the operators from their respective sides. In math, this reads + + <\equation*> + >=A>>||\>=A\>.>>>>> + + + This is very interesting set of conditions: it seemingly over-determines + the system, and the R-H conditions become\ + + <\equation*> + |A-\I>->||>|||>|>>||\>>>|||>|A-\I\>->||>>>> + + + \; + + This means that the quantity >-> + *must* be a multiple of >. + I.e.,\ + + <\equation*> + >=\+\+\, + + + for some unknown scalar >. A similar argument leads us + to\ + + <\equation*> + \>=\+\+\. + + + It is worth stressing that our unknown vectors + >> and + \>>, + previously a full 6 unknowns, are now merely two unknowns, + > and >. This will be true for any + size of the matrices >>.\ + + \; + + =0> + and =0>. If we impose directional bias, `Stoopid + upwinding' automatically satisfies all but one of the R-H conditions: the + one at =0>.> + + \; + + Of course, we are still overdetermined: we have one R-H condition left: + three equations, two unknowns. What's left? The =0> + conditions boils down to\ + + <\equation*> + \\*-\*\=\\-\\, + + + where the unknowns are > and >. + Note that we *need* this equation to be satisifed independent of the + unknowns > and >, + which are dependent on the incoming waves + >>, over which we have no + control. This means that we require that + , + \span,>. + Note that since \>, + and \>, + then we must have that \> + and \>. + Then (and only then!) will we be able to solve for the fluxes we wish for + any values of >>>.\ + + \; + + In this simple case, our condition states that + =c*> + and =c*>. + I.e., that left traveling waves stay left-traveling waves and + right-traveling waves do the same. Then we have that\ + + <\equation*> + =c*\*|\>>||=\|c*\>>>>>>, + + + which uniquely determines the flux as\ + + <\equation*> + A=A>=A\>=\\+\\=\\+\\. + + + \; + + The above considerations lead us to the following theorem: + + <\theorem> + (Well-posedness of Rankine-Hugoniot directionally-biased upwinding) + + Let >(|^>)> + represent the (discontinuous) operator of a multidimensional strongly + hyperbolic linear wave equation at an element interface. Assume + diagonalizations of the form\ + + <\equation*> + >=V>\>V>,>>|>|>=>>|>>|>|>>>>>>>>>>> + + + \; + + Let the evolving \>. + The known values of > at the element + interfaces admit the decomposition\ + + <\equation*> + =\>||=\.>>>>> + + + Assume that each diagonal matrix >> contains + eigenvalues >> such that the + following quantities are well-defined: + + <\equation*> + =i: + \\0>||=i: + \\0>||=i: + \=0>>|||||>|=i: + \\0>||=i: + \\0>||=i: + \=0>>|||||>|=>||=>||=.>>>>> + + + Note that this gives us . Associate with each eigenvalue + >> the eigenvector + >>. Define the + following vector subspaces of >: + + <\equation*> + >=span:i\>||>=span:i\>||>=span:i\>>|||||>|>=span:i\>||>=span:i\>||>=span:i\>>>>> + + + \; + + Define the intermediate states + (i)>> + and the associated fluxes A(i)>> + such that (i)>> and + A(i)>> + `lies between' eigenvalues >> and + >> (see figure + ). Now define the unsuperscripted eigenvalues as\ + + <\equation*> + \=sort\: + i\>>|>|\=sort\: + i\>>|>|\=sort\: + i\>>>>> + + + \; + + Define the numerical flux as\ + + <\equation*> + A\A(q)>, + + + where \>.\ + + \; + + Assume the following conditions: + + <\enumerate> + 0> + + Imposition of the Rankine-Hugoniot condition across eigenvalues: + + <\equation*> + \(i)>-(i-1)> + = A(i)>-A(i-1)> + \ \ \ i = 1, 2, \n, + + + where we define (0)>=> + and (n)>=>.\ + + Directional biasing: + + <\equation*> + (A)(i)>=(i)>>||||,N>>|||||>|(i)>>||||,n-P>>>>> + + + >=>> + and >=>> + \ (>=>>) + + + \; + + Then the numerical flux A> + exists and is unique. It is given by\ + + <\equation*> + A>||\\+\\>>|||>|||\\+\\,>>>>> + + + \; + + where the unknowns \> + satisfy the rank system\ + + <\equation*> + >|>||>||>>>>>>>|>>|>>|>>>>>=\\, + + + and the unknowns \> + satisfy the rank system\ + + <\equation*> + >|>||>||>>>>>>>|>>|>>|>>>>>=\\. + + + + Note that this theorem even takes into account degeneracies in any + eigenvalue. + + + + We're dealing with + + <\equation*> + u=c(x)\u, + + + which we rewrite into conservation form: + + <\eqnarray*> + -c(x)\\\>||>|-c(x)\u>||>>> + + + and from there into matrix form: + + <\eqnarray*> + +c(x)||>|||>|||>>>>|\>\>\+||>|||>|||>>>>|\>\>\>||>>>> + + + where we set \(u,\)>.\ + + (see and + for Maxima scripts for code + to find the upwind flux.) + + We find: + + <\eqnarray*> + *u)>>||\>+-[c*u]|2>>>|\)>>||+(\\\)\[c\-c\].>>>> + + + + + We're dealing with + + <\equation*> + u=\\\u, + + + which we rewrite into matrix conservation form as + + <\equation*> + u=vu+vu. + + + We find )=\\\>, and therefore + + <\equation*> + (A(\)u)>=\\u>|\\\0> + (outflow)>,>>|\\u>|\\\0> (inflow)>.>>>>> + + + + + We're discretizing +v\\u=0>. Recall + + <\eqnarray*> + \(v*u*\)>||\(v*u)\+v*u\\\=(v\\*u)\+v*u\\\.>>>> + + + Using the summation convention, we find + + <\eqnarray*> + ||>u\+>(v\\u)\>>|||>u\->v*u\\\+>\\(v*u*\)>>|||>u\->v*u\\\+T>v*u\\n>>||>|>u\->v*u\\\+T>(v*u)>\\n>>||>|>\ull->ul>>|ul>>>>>\>l>>|>l>>>>>+\T>u)>ll>>|u)>ll>>>>>\n>>|||>|||>>> + + + <\eqnarray*> + ||>u\->a*u\\\+T>(v*u)>\\n>>|||>u\+>(v\\u)\-T>(v*u-{v*u})\\n>>|||>u\+>(v\\u)\-T>v>>|>>>>>\>>|>>>>>\\n>>>> + + + +<\initial> + <\collection> + + + + + + + +<\references> + <\collection> + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + > + + + +<\auxiliary> + <\collection> + <\associate|figure> + Space-Time Diagram of a flux + fan.|> + + <\associate|toc> + |math-font-series||1Mapping + from 2D equilateral to Unit Coordinates> + |.>>>>|> + + + |math-font-series||2Mapping + from 3D equilateral to Unit Coordinates> + |.>>>>|> + + + |math-font-series||3Derivatives + of the 2D Basis functions> |.>>>>|> + + + |math-font-series||4Derivatives + of the 3D Basis functions> |.>>>>|> + + + |math-font-series||5DG + Schemes> |.>>>>|> + + + |5.1Notation + |.>>>>|> + > + + |5.2Advection Equation + |.>>>>|> + > + + |5.2.1Weak DG + |.>>>>|> + > + + |5.2.2Strong DG + |.>>>>|> + > + + |5.3Wave Equation + |.>>>>|> + > + + |5.3.1Weak DG + |.>>>>|> + > + + |5.3.2Strong DG + |.>>>>|> + > + + |5.4Heat Equation + |.>>>>|> + > + + |math-font-series||6Quadrature + Rules> |.>>>>|> + + + |math-font-series||7Upwind + Fluxes> |.>>>>|> + + + |7.1DIY Upwind Fluxes/Riemann + Solvers |.>>>>|> + > + + |7.2Comments by Akil + |.>>>>|> + > + + |7.2.1Zero eigenvalues + |.>>>>|> + > + + |7.2.2Degenerate eigenvalues + |.>>>>|> + > + + |7.3More random comments by + Akil |.>>>>|> + > + + |7.3.1Stoopid upwinding + |.>>>>|> + > + + |7.3.2R-H Conditions with + Directional Bias |.>>>>|> + > + + |7.4Deriving an Upwind Flux for + the Wave Equation |.>>>>|> + > + + |7.5Deriving an Upwind Flux for + the 2D Advection Equation |.>>>>|> + > + + + \ No newline at end of file diff --git a/contrib/notes/ip-primal-nu.tm b/contrib/notes/ip-primal-nu.tm new file mode 100644 index 00000000..eafbd0dc --- /dev/null +++ b/contrib/notes/ip-primal-nu.tm @@ -0,0 +1,89 @@ + + + + +<\body> + Variational formulation of Poisson with fluxes > + and |^>>, solving for + ,\>: + + <\eqnarray*> + >\\\>||>u\\\+K>\\\,>>|>\\\v>||>vf+K>v|^>\\.>>>> + + + Use + + <\equation*> + K>v\\\=>[v]\{\}+>{v}[\], + + + where > includes the interior faces, to convert to: + + <\eqnarray*> + >\\\>||>u\\\+>[]\{\}+>{}[\],>>|>\\>\v>||>vf+>[v]\{|^>}+>{v}[|^>].>>>> + + + Now pick \>\v> + to facilitate conversion to a bilinear form. This turns the first equation + into. + + <\eqnarray*> + >\\\v>||>u\\(>\v)+>[]\{>\v}+>{}[>\v].>>>> + + + We can now equate the RHSs of this latter equation and the second equation + above: + + <\equation*> + >vf+>[v]\{|^>}+>{v}[|^>]=->u\\(>\v)|\>+>[]\{>\v}+>{}[>\v]. + + + Now use the integration-by-parts formula + + <\equation*> + >\v\\+>v\\\=>[v]\{\}+>{v}[\] + + + to eliminate the double derivative, with + =\\v> and >. + + <\equation*> + A=>u\\(>\v)=>(\u)\(>\v)+>[u]\{>\v}+>{u}[>\v]> + + + and obtain + + <\equation*> + >vf+>[v]\{|^>}+>{v}[|^>]=>>\u\\v->[u]\{>\v}->{u}[>\v]>+>[]\{>\v}+>{}[>\v] + + + Regroup and reorder: + + <\equation*> + >>\u\\v+>{-u}[>\v]+>[-u]\{>\v}->[v]\{|^>}->{v}[|^>]=>vf. + + + For IP, e.g. + + <\equation*> + \{u},|^>\{>\u}-|h>>[u]. + + + Substitute in: + + <\equation*> + >>\u\\v+>{}>-u}[>\v]|\>+>[}>-u]\{>\v}->[v]\{>\u}-|h>>[u]}>->{v}[>\u}-|h>>[u]>]=>vf + + + and obtain: + + <\equation*> + >>\u\\v->[u]\{>\v}->[v]\{>\u}->[v]|h>>[u]=>vf. + + + +<\initial> + <\collection> + + + \ No newline at end of file diff --git a/contrib/notes/ip-primal.tm b/contrib/notes/ip-primal.tm new file mode 100644 index 00000000..17fa4041 --- /dev/null +++ b/contrib/notes/ip-primal.tm @@ -0,0 +1,87 @@ + + + + +<\body> + Variational formulation of Poisson with fluxes > + and |^>>, solving for + ,\>: + + <\eqnarray*> + >\\\>||>u\\\+K>\\\,>>|>\\\v>||>vf+K>v|^>\\.>>>> + + + Use + + <\equation*> + K>v\\\=>[v]\{\}+>{v}[\], + + + where > includes the interior faces, to convert to: + + <\eqnarray*> + >\\\>||>u\\\+>[]\{\}+>{}[\],>>|>\\\v>||>vf+>[v]\{|^>}+>{v}[|^>].>>>> + + + Now pick \\v> to facilitate + conversion to a bilinear form. This turns the first equation into. + + <\eqnarray*> + >\\\v>||>u\\(\v)+>[]\{\v}+>{}[\v].>>>> + + + We can now equate the RHSs of this latter equation and the second equation + above: + + <\equation*> + >vf+>[v]\{|^>}+>{v}[|^>]=->u\\(\v)|\>+>[]\{\v}+>{}[\v]. + + + Now use the integration-by-parts formula + + <\equation*> + >\v\\+>v\\\=>[v]\{\}+>{v}[\] + + + to eliminate the double derivative: + + <\equation*> + A=>u\\(\v)=>\u\\v+>[u]\{\v}+>{u}[\v]> + + + and obtain + + <\equation*> + >vf+>[v]\{|^>}+>{v}[|^>]=>\u\\v->[u]\{\v}->{u}[\v]>+>[]\{\v}+>{}[\v] + + + Regroup and reorder: + + <\equation*> + >\u\\v+>{-u}[\v]+>[-u]\{\v}->[v]\{|^>}->{v}[|^>]=>vf. + + + For IP, e.g. + + <\equation*> + \,|^>\{\u}-|h>[u]. + + + Substitute in: + + <\equation*> + >\u\\v+>{}>-u}[\v]|\>+>[}>-u]\{\v}->[v]\{u}-|h>[u]}>->{v}[u}-|h>[u]>]=>vf + + + and obtain: + + <\equation*> + >\u\\v->[u]\{\v}->[v]\{\u}->[v]|h>[u]=>vf. + + + +<\initial> + <\collection> + + + \ No newline at end of file diff --git a/contrib/notes/mystyle.ts b/contrib/notes/mystyle.ts new file mode 100644 index 00000000..e908e301 --- /dev/null +++ b/contrib/notes/mystyle.ts @@ -0,0 +1,32 @@ + + + + +<\body> + <\with|mode|math> + ,>>>> + + >>>> + + >>>> + + >>>> + + >> + + >> + + | ->>>> + + > + + + >)>> + + +<\initial> + <\collection> + + + + \ No newline at end of file -- GitLab