Hyperelastic 2D Riemann problem solver

Define P(F), (Subscript[∇, F])P

SetDirectory[NotebookDirectory[]]

Import 2d flux definitions

NotebookEvaluate["/home/student/ShearWaveAIR/NumericalModels/bearclaw/FluxDefinitions2d.nb"];

Eigensystem

Velocity eigenvector block

Find squared speeds

Form matrix

M=Table[H[[k]].U[[k]],{k,d}];
For[k=1,k≤d,k++,
For[i=1,i≤d,i++,
For[j=i+1,j≤d,j++,
M[[k,i,j]]=M[[k,j,i]]]]];
Map[TableForm[#]&,M]

Find squared wave speeds S=(s^2), i.e., the eigenvalues of M

s2=Map[Reverse[Simplify[Eigenvalues[#]]]&,M];
s2[[2]]=Reverse[s2[[2]]];
Map[TableForm[#]&,s2]

Define eigenvectors of M

Define eigenvectors in a form that reduces to a non-singular diagonal matrix for (Subscript[P, 3,1])=0, and (Subscript[P, 1,1])>(Subscript[P, 3,3]) along direction 1, and for (Subscript[P, 4,2])=0, and (Subscript[P, 4,4])>(Subscript[P, 2,2]) along direction 2

wrow = Table[Simplify[{{s2[[k,1]]-M[[k,2,2]],M[[k,2,1]]},{M[[k,2,1]],s2[[k,2]]-M[[k,1,1]]}}],{k,d}];
Map[TableForm[#]&,wrow]

Check the behavior

Check that the above eigenvector form simplifies correctly for  (Subscript[P, 3,1])=0, and (Subscript[P, 1,1])>(Subscript[P, 3,3]) along direction 1, and for  (Subscript[P, 4,2])=0, and (Subscript[P, 4,4])>(Subscript[P, 2,2]) along direction 2

FullSimplify[wrow[[1]] /. M[[1,2,1]]->0,M[[1,1,1]]>M[[1,2,2]]]

FullSimplify[wrow[[2]] /. M[[2,2,1]]->0,M[[2,2,2]]>M[[2,1,1]]]

Note that a singular diagonal matrix is obtained for (Subscript[P, 1,1])<(Subscript[P, 3,3]) along direction 1, (Subscript[P, 4,4])<(Subscript[P, 2,2]) along direction 2

FullSimplify[wrow[[1]] /. M[[1,2,1]]->0,M[[1,1,1]]<M[[1,2,2]]]

FullSimplify[wrow[[2]] /. M[[2,2,1]]->0,M[[2,2,2]]<M[[2,1,1]]]

Verify velocity block eigenmode relation

Table[Simplify[M[[k]].wrow[[k,i]] - s2[[k,i]] wrow[[k,i]]],{k,d},{i,d}] // TableForm

Replacement rules

Replacement rules for common subexpressions

Δ=Table[Discriminant[CharacteristicPolynomial[M[[k]],S],S],{k,d}]

sdRules = Table[Sqrt[Δ[[k]]]->sd,{k,d}]

trM=Table[Tr[M[[k]]],{k,d}]

trMRules = Table[trM[[k]]->tr,{k,d}]

MtrRules = Map[Reverse[#]&,trMRules]

StoPRules =  Table[(Subscript[S, i])-> s2[[k,i]] /. sdRules[[k]] /. trMRules[[k]],{k,d},{i,d}]

sdTotrRules=Map[Flatten[#]&,Table[Solve[ (Subscript[S, i]) == ((Subscript[S, i]) /. StoPRules[[k,i]]),sd ],{k,d},{i,d}]]

Column representation of eigenvectors

wSrow=Simplify[Table[wrow[[k,i]] /. sdRules[[k]] /. trMRules[[k]] /. sdTotrRules[[k,i]] /. MtrRules[[k]],{k,d},{i,d}]];
w=Map[Transpose[#]&,wSrow];
Map[TableForm[#]&,w]

Propagating eigenmodes

Backward/forward propagating block

𝒱[k_,s_,w_]:=Transpose[{Join[-s w,U[[k]].w],Join[s w,U[[k]].w]}];

Along direction 1

𝒱[1,(Subscript[s, 1]),{(Subscript[v, 1]),(Subscript[v, 2])}] // TableForm

Along direction 2

𝒱[2,(Subscript[s, 1]),{(Subscript[v, 1]),(Subscript[v, 2])}] // TableForm

Full, propagating eigenmodes

Define some utility functions.

RepeatWave[w_]:=Module[{ww},
l=Length[w]; ww={};
For[i=1,i≤l,i++,
AppendTo[ww,w[[i]]];AppendTo[ww,w[[i]]]];
Return[ww];
]
FormWave[s_,w_,U_]:= Flatten[{s w,U. w}];

Form velocity block of propagating waves

ww=Table[RepeatWave[wSrow[[k]]],{k,d}];
Dimensions[ww]

Transpose[ww[[1]]] // TableForm

Form speeds

speeds =Map[Flatten[#]&,Table[{-Sqrt[s2[[k,i]]],Sqrt[s2[[k,i]]]},{k,d},{i,d}]];

Simplify[speeds[[1]] /. Prules,{λ>0,μ>0,λ>μ}]

Form propagating waves

Wrow =Table[FormWave[speeds[[k,i]],ww[[k,i]],U[[k]]],{k,d},{i,2d}]; (* Eigenvectors stored in rows *)
W=Map[Transpose[#]&,Wrow]; (* Eigenvectors stored in columns *)

TableForm[W[[1]]]

Verify eigenvector relationship

Table[Simplify[A[k].Wrow[[k,;;,i]]-speeds[[k,i]] Wrow[[k,;;,i]]],{k,d},{i,2d}];

Simplified form of eigenvectors with identification of common subexpressions

Define rules to reflect that (Subscript[∇, F])P=(Subscript[∇, F])(Subscript[∇, F])E is symmetric

symPrules = Flatten[Table[(Subscript[P, j,i])->(Subscript[P, i,j]),{i,2d},{j,i+1,2d}]]

Eigenvectors in which common subexpressions have been identified

WS =Map[Factor[#]&, Table[Simplify[W[[k]]/2 /. dRules[[k]] /. sdTosRules[[k]] /. sdTosRules2[[k]] /. sdRules[[k]] ] /. s12Rules[[k]] ,{k,d}]]  /. symPrules;

Verify simplified form of eigenvectors

Table[Simplify[A[k].WS[[k,;;,i]]-speeds[[k,i]] WS[[k,;;,i]] /. SRules[[k]] /. sRules[[k]] /. symPrules],{k,d},{i,2d}]

Verify orthogonality of simplified form of eigenvectors

Table[FullSimplify[WS[[k,;;,1]].WS[[k,;;,3]] /. SRules[[k]] /. sRules[[k]]] /. symPrules,{k,d}]

Table[FullSimplify[WS[[k,;;,1]].WS[[k,;;,4]] /. SRules[[k]] /. sRules[[k]]] /. symPrules,{k,d}]

Table[FullSimplify[WS[[k,;;,2]].WS[[k,;;,3]] /. SRules[[k]] /. sRules[[k]]] /. symPrules,{k,d}]

Table[FullSimplify[WS[[k,;;,2]].WS[[k,;;,4]] /. SRules[[k]] /. sRules[[k]]] /. symPrules,{k,d}]

Display final, simplified form of eigensystem

Display the simplified form of eigenvalues (wave speeds). First rules to form square root of discriminant

dRules

Next rules to find the positive speeds

sdTosRules

Inversion of the above rules

sTosdRules = Map[Reverse[#]&,sdTosRules,{2}]

Additional rules to obtain self-contained eigenmode blocks

EigRules = Map[{#[[1]],#[[1]],#[[2]],#[[2]]}&,Map[Reverse[Flatten[#]]&,Table[Solve[Simplify[((Subscript[s, 1])^2)+((Subscript[s, 2])^2) /. sTosdRules[[k]] ]== (Subscript[sqs, 1])+(Subscript[sqs, 2]),(Subscript[sqs, i])],{k,d},{i,d}] /. Table[(Subscript[sqs, k])->((Subscript[s, k])^2),{k,d}]]]

𝒲=Map[Transpose[#]&,Table[WS[[k,;;,i]] /. EigRules[[k,i]],{k,d},{i,2d}]];

Display the eigenvectors along direction 1

TableForm[𝒲[[1]],TableDirections->Column]

Compare with

TableForm[WS[[1]],TableDirections->Column]

Display the eigenvectors along direction 2

TableForm[𝒲[[2]],TableDirections->Column]

Compare with

TableForm[WS[[2]],TableDirections->Column]

Null space

Find the null space of (H^T) to complete the eigenbasis

NHT=Table[NullSpace[Transpose[U[[k]]]],{k,d}];
𝒩=Map[Transpose[#]&,Map[Join[{0,0},#]&,NHT,{2}]];
Map[TableForm[#,TableDirections->Column]&,𝒩]

Verify that 𝒲, 𝒩 are orthogonal

Table[𝒲[[k,;;,i]] . 𝒩[[k,;;,j]],{k,d},{i,2d},{j,(d^2)-d}]

Riemann problem

Eigenvalues in terms of deformation gradient

dRules

sdTosRules

Display the eigenvectors along direction 1

TableForm[𝒲[[1]],TableDirections->Column]

Display the eigenvectors along direction 2

TableForm[𝒲[[2]],TableDirections->Column]

Decompose jump onto propagating waves

VV=Table[Simplify[(Transpose[W[[k]]] . W[[k]])],{k,d}];
Dimensions[PW]

symRules = { (Subscript[P, 3,1])->(Subscript[P, 1,3]),(Subscript[P, 4,2])->(Subscript[P, 2,4])}

symdRules = dRules /. symRules

sdRules[[1]]

TableForm[Factor[PW[[1]] /. symRules /. symdRules[[1]] /. sdRules[[1]]],TableDirections->Column]

Simplify[PW[[1]] /.dRules[[1]] /. sRules[[1]] /. sRules2[[1]] /. sdRules[[1]]  /. c12Rules[[1]] /. symPrules]

PWS = Map[Factor[#]&,Table[ Simplify[PW[[k]] /. dRules[[k]] /. sRules[[k]] /. sRules2[[k]] /. sdRules[[k]] /. c12Rules[[k]] ],{k,d} ]];
Dimensions[PWS]

TableForm[PWS[[1]],TableDirections->Column]

PW=Table[Simplify[(Transpose[WS[[k]]] . WS[[k]]) /. c2Rules[[k]] /. cRules[[k]]  /. dRules[[k]] /. sRules[[k]] /. sRules2[[k]] /. sdRules[[k]]  /. c12Rules[[k]]  /. symPrules],{k,d}];

TableForm[PW[[1]],TableDirections->Column]

General form

Define a function to return the speeds and eigenmodes of the conservative system (Subscript[q, t])+(Subscript[A, 1])(Subscript[q, ,1])+(Subscript[A, 2])(Subscript[q, ,2])=0.  The flux Jacobian is of dimensions m×m, m=d+(d^2) and has 2d propagating modes, hence the number of non-propagating modes (with zero eigenvalue) is (d^2)-d

Modes[A_]:=Module[{m,eigs,vecs,cRules,rules,denom,k,R,lam,vec,modes,speeds},
m=Length[A]; modes={};speeds={};
eigs=Eigenvalues[A];
cRules =Table[ eigs[[(d^2)-d+2i]]->(Subscript[c, i]),{i,d}];
rules = {Simplify[((eigs[[(d^2)-d+4]])^2)- ((eigs[[(d^2)-d+2]])^2)] ->((Subscript[c, 2])^2)-((Subscript[c, 1])^2)};
R=Eigenvectors[A];
Monitor[
For[k=1,k≤m,k++,
lam=Simplify[eigs[[k]] /. cRules /. rules];
vec= R[[k]] /. cRules  /. rules; Print["vec=",vec];
 vec = Simplify[Denominator[Together[Total[vec]]] vec]; Print["vec=",vec];
AppendTo[modes,vec];
]
,k];
Return[{eigs,modes}];
]

cRules =Table[ eigs[[(d^2)-d+2i]]->(Subscript[c, i]),{i,d}]

rules = {Simplify[((eigs[[(d^2)-d+4]])^2)- ((eigs[[(d^2)-d+2]])^2)] ->((Subscript[c, 2])^2)-((Subscript[c, 1])^2)}

eigs=Eigenvalues[A[1]] /. cRules /. rules

tt=Eigenvectors[A[1]] /. cRules /. rules

s=sR[[;;,1]]; Rrow=sR[[;;,2]];

Verify that the first (d^2)-d eigenvectors indeed correspond to null eigenvalues

Table[Simplify[A[i].Rrow[[i,j]]],{i,d},{j,(d^2)-d}]

Mathematica stores eigenvectors row-wise. Switch to columns.

R=Transpose[Rrow,{1,3,2}];

Display the eigensystem along the first spatial direction. Eigenvalues (speeds)

s[[1]]

Display the eigenmode along the first spatial direction

R[[1]] // TableForm

Eigenvalues along the second spatial direction

s[[2]]

Eigenmodes along second spatial direction

R[[2]] // TableForm

Verification for linear elastic

(Subscript[P, EL])=(Subscript[P, EL])=Simplify[λ (Sum[(Subscript[F, idx2[l,l]]),{l,d}]-d)δ+μ Table[(Subscript[F, idx2[i,k]])+(Subscript[F, idx2[k,i]])-2 δ[[i,k]],{i,d},{k,d}]] /. LinearizeFTensor;
TableForm[(Subscript[P, EL])]Simplify[λ (Sum[(Subscript[F, idx2[l,l]]),{l,d}]-d)δ+μ Table[(Subscript[F, idx2[i,k]])+(Subscript[F, idx2[k,i]])-2 δ[[i,k]],{i,d},{k,d}]] /. LinearizeFTensor;
TableForm[(Subscript[P, EL])]

StressGradientRules=Flatten[Table[(Subscript[P, d(i-1)+j,k])->D[(Subscript[P, EL])[[i,j]],(Subscript[F, k])],{i,d},{j,d},{k,(d^2)}]]

Simplify[s[[1]] /. StressGradientRules, {λ>0,μ>0}] // TableForm

Simplify[R[[1]] /. StressGradientRules /. {(Subscript[c, 1])->(Sqrt[μ]),(Subscript[c, 2])->(Sqrt[λ+2μ])},{λ>0,μ>0}]// TableForm

Left/Right states

Define left/right eigenbases

Lrules= Flatten[{Table[(Subscript[P, i,j])->(Subscript[P, i,j]^L),{i,(d^2)},{j,(d^2)}],{(Subscript[c, 1])->(Subscript[c, 1]^L),(Subscript[c, 2])->(Subscript[c, 2]^L)}}];
Rrules= Flatten[{Table[(Subscript[P, i,j])->(Subscript[P, i,j]^R),{i,(d^2)},{j,(d^2)}],{(Subscript[c, 1])->(Subscript[c, 1]^R),(Subscript[c, 2])->(Subscript[c, 2]^R)}}];
RL = R /. Lrules;
RR = R /. Rrules;

Form interface Jacobian

J =(RL+RR)/2; (* Start from arithmetic average of L/R states *)
(* Left-going modes *)
J[[;;,;;,3]]=RL[[;;,;;,3]];  J[[;;,;;,5]]=RL[[;;,;;,5]]; 
(* Right-going modes *)
J[[;;,;;,4]]=RL[[;;,;;,4]];  J[[;;,;;,6]]=RL[[;;,;;,6]]; 
J[[1]] //TableForm

Compute Jacobian determinants

detJ=Table[Simplify[Det[J[[k]]]],{k,d}]

Compute Jacobian inverse

Jinv =Monitor[ Table[Simplify[Inverse[J[[k]]]],{k,d}],k];

Table[MatrixForm[Simplify[Jinv[[k]].J[[k]]]],{k,d}]