How to use nInt method of td Package

Best Go-testdeep code snippet using td.nInt

euler_test.go

Source:euler_test.go Github

copy

Full Screen

1package Euler2D2import (3 "fmt"4 "math"5 "strconv"6 "sync"7 "testing"8 "github.com/notargets/gocfd/DG2D"9 "github.com/notargets/gocfd/types"10 "github.com/stretchr/testify/assert"11 "github.com/notargets/gocfd/utils"12)13var ipDefault = &InputParameters{14 Title: "",15 CFL: 1,16 FluxType: "Lax",17 InitType: "freestream",18 PolynomialOrder: 0,19 FinalTime: 1,20 Minf: 0,21 Gamma: 1.4,22 Alpha: 0,23 BCs: nil,24 LocalTimeStepping: false,25 MaxIterations: 5000,26 ImplicitSolver: false,27 Limiter: "",28}29func TestFluidFunctions(t *testing.T) {30 N := 131 plotMesh := false32 ip := *ipDefault33 ip.Minf = 2.34 ip.PolynomialOrder = N35 c := NewEuler(&ip, "../../DG2D/test_tris_6.neu", 1, plotMesh, false, false)36 funcs := []FlowFunction{Density, XMomentum, YMomentum, Energy, Mach, StaticPressure}37 values := make([]float64, len(funcs))38 for i, plotField := range []FlowFunction{Density, XMomentum, YMomentum, Energy, Mach, StaticPressure} {39 values[i] = c.FSFar.GetFlowFunction(c.Q[0], 0, plotField)40 //fmt.Printf("%s[%d] = %8.5f\n", plotField.String(), ik, values[i])41 }42 assert.InDeltaSlicef(t, []float64{1, 2, 0, 3.78571, 2, 0.71429}, values, 0.00001, "err msg %s")43}44func TestEuler(t *testing.T) {45 var (46 msg = "err msg %s"47 tol = 0.00000148 )49 ip := *ipDefault50 ip.FluxType = "average"51 if true {52 { // Test interpolation of solution to edges for all supported orders53 Nmax := 754 for N := 0; N <= Nmax; N++ {55 ip.PolynomialOrder = N56 //c := NewEuler(1, N, "../../DG2D/test_tris_5.neu", 1, FLUX_Average, FREESTREAM, 1, 0, 1.4, 0, false, 5000, None, false, false, false)57 c := NewEuler(&ip, "../../DG2D/test_tris_5.neu", 1, false, false, false)58 Kmax := c.dfr.K59 Nint := c.dfr.FluxElement.NpInt60 Nedge := c.dfr.FluxElement.NpEdge61 var Q, Q_Face [4]utils.Matrix62 for n := 0; n < 4; n++ {63 Q[n] = utils.NewMatrix(Nint, Kmax)64 Q_Face[n] = utils.NewMatrix(3*Nedge, Kmax)65 }66 for n := 0; n < 4; n++ {67 for i := 0; i < Nint; i++ {68 for k := 0; k < Kmax; k++ {69 ind := k + i*Kmax70 Q[n].DataP[ind] = float64(k + 1)71 }72 }73 }74 // Interpolate from solution points to edges using precomputed interpolation matrix75 for n := 0; n < 4; n++ {76 Q_Face[n] = c.dfr.FluxEdgeInterp.Mul(Q[n])77 }78 for n := 0; n < 4; n++ {79 for i := 0; i < 3*Nedge; i++ {80 for k := 0; k < Kmax; k++ {81 ind := k + i*Kmax82 assert.InDeltaf(t, float64(k+1), Q_Face[n].DataP[ind], tol, msg)83 }84 }85 }86 }87 }88 }89 if true {90 { // Test solution process91 /*92 Solver approach:93 0) Solution is stored on sol points as Q94 0a) Flux is computed and stored in X, Y component projections in the 2*NpInt front of F_RT_DOF95 1) Solution is extrapolated to edge points in Q_Face from Q96 2) Edges are traversed, flux is calculated and projected onto edge face normals, scaled and placed into F_RT_DOF97 */98 Nmax := 799 for N := 0; N <= Nmax; N++ {100 //c := NewEuler(1, N, "../../DG2D/test_tris_5.neu", 1, FLUX_Average, FREESTREAM, 1, 0, 1.4, 0, false, 5000, None, false, false, false)101 ip.PolynomialOrder = N102 c := NewEuler(&ip, "../../DG2D/test_tris_5.neu", 1, false, false, false)103 Kmax := c.dfr.K104 Nint := c.dfr.FluxElement.NpInt105 Nedge := c.dfr.FluxElement.NpEdge106 NpFlux := c.dfr.FluxElement.Np // Np = 2*NpInt+3*NpEdge107 var Q_Face, F_RT_DOF [4]utils.Matrix108 for n := 0; n < 4; n++ {109 Q_Face[n] = utils.NewMatrix(3*Nedge, Kmax)110 F_RT_DOF[n] = utils.NewMatrix(NpFlux, Kmax)111 }112 Q := c.Q[0]113 // Mark the initial state with the element number114 qD := [4][]float64{Q[0].DataP, Q[1].DataP, Q[2].DataP, Q[3].DataP}115 for i := 0; i < Nint; i++ {116 for k := 0; k < Kmax; k++ {117 ind := k + i*Kmax118 qD[0][ind] = float64(k + 1)119 qD[1][ind] = 0.1 * float64(k+1)120 qD[2][ind] = 0.05 * float64(k+1)121 qD[3][ind] = 2.00 * float64(k+1)122 }123 }124 // Flux values for later checks are invariant with i (i=0)125 Fr_check, Fs_check := make([][4]float64, Kmax), make([][4]float64, Kmax)126 for k := 0; k < Kmax; k++ {127 Fr_check[k], Fs_check[k] = c.CalculateFluxTransformed(k, Kmax, 0, c.dfr.Jdet, c.dfr.Jinv, Q)128 }129 // Interpolate from solution points to edges using precomputed interpolation matrix130 for n := 0; n < 4; n++ {131 Q_Face[n] = c.dfr.FluxEdgeInterp.Mul(Q[n])132 }133 // Calculate flux and project into R and S (transformed) directions134 rtD := [4][]float64{F_RT_DOF[0].DataP, F_RT_DOF[1].DataP, F_RT_DOF[2].DataP, F_RT_DOF[3].DataP}135 for n := 0; n < 4; n++ {136 for i := 0; i < Nint; i++ {137 for k := 0; k < Kmax; k++ {138 ind := k + i*Kmax139 Fr, Fs := c.CalculateFluxTransformed(k, Kmax, i, c.dfr.Jdet, c.dfr.Jinv, Q)140 rtD[n][ind], rtD[n][ind+Nint*Kmax] = Fr[n], Fs[n]141 }142 }143 // Check to see that the expected values are in the right place (the internal locations)144 rtTD := F_RT_DOF[n].Transpose().DataP145 for k := 0; k < Kmax; k++ {146 val0, val1 := Fr_check[k][n], Fs_check[k][n]147 is := k * NpFlux148 assert.True(t, nearVecScalar(rtTD[is:is+Nint], val0, 0.000001))149 is += Nint150 assert.True(t, nearVecScalar(rtTD[is:is+Nint], val1, 0.000001))151 }152 // Set normal flux to a simple addition of the two sides to use as a check in assert()153 for k := 0; k < Kmax; k++ {154 for i := 0; i < 3*Nedge; i++ {155 ind := k + (2*Nint+i)*Kmax156 Fr, Fs := c.CalculateFluxTransformed(k, Kmax, i, c.dfr.Jdet, c.dfr.Jinv, Q_Face)157 rtD[n][ind] = Fr[n] + Fs[n]158 }159 }160 // Check to see that the expected values are in the right place (the edge locations)161 rtTD = F_RT_DOF[n].Transpose().DataP162 for k := 0; k < Kmax; k++ {163 val := Fr_check[k][n] + Fs_check[k][n]164 is := k * NpFlux165 ie := (k + 1) * NpFlux166 assert.True(t, nearVecScalar(rtTD[is+2*Nint:ie], val, 0.000001))167 }168 }169 }170 }171 }172 if true {173 { // Test solution process part 2 - Freestream divergence should be zero174 Nmax := 7175 for N := 0; N <= Nmax; N++ {176 //c := NewEuler(1, N, "../../DG2D/test_tris_5.neu", 1, FLUX_Average, FREESTREAM, 1, 0, 1.4, 0, false, 5000, None, false, false, false)177 ip.PolynomialOrder = N178 //ip.Minf = 0.301179 //ip.Alpha = 2180 plotMesh := false181 //c := NewEuler(&ip, "../../DG2D/test_tris_5.neu", 1, plotMesh, false, false)182 c := NewEuler(&ip, "../../DG2D/test_tris_6_nowall.neu", 1, plotMesh, false, false)183 c.FSIn = c.FSFar184 Kmax := c.dfr.K185 Nint := c.dfr.FluxElement.NpInt186 Nedge := c.dfr.FluxElement.NpEdge187 NpFlux := c.dfr.FluxElement.Np // Np = 2*NpInt+3*NpEdge188 // Mark the initial state with the element number189 var Q_Face, F_RT_DOF [4]utils.Matrix190 var Flux, Flux_Face [2][4]utils.Matrix191 for n := 0; n < 4; n++ {192 F_RT_DOF[n] = utils.NewMatrix(NpFlux, Kmax)193 Q_Face[n] = utils.NewMatrix(3*Nedge, Kmax)194 Flux_Face[0][n] = utils.NewMatrix(3*Nedge, Kmax)195 Flux_Face[1][n] = utils.NewMatrix(3*Nedge, Kmax)196 Flux[0][n] = utils.NewMatrix(Nint, Kmax)197 Flux[1][n] = utils.NewMatrix(Nint, Kmax)198 }199 Q := c.Q[0]200 c.SetRTFluxInternal(Kmax, c.dfr.Jdet, c.dfr.Jinv, F_RT_DOF, Q)201 c.InterpolateSolutionToEdges(Q, Q_Face, Flux, Flux_Face)202 EdgeQ1 := make([][4]float64, Nedge)203 EdgeQ2 := make([][4]float64, Nedge)204 c.CalculateEdgeFlux(0, false, nil, nil, [][4]utils.Matrix{Q_Face},205 [][2][4]utils.Matrix{Flux_Face}, c.SortedEdgeKeys[0], EdgeQ1, EdgeQ2)206 c.SetRTFluxOnEdges(0, Kmax, F_RT_DOF)207 // Check that freestream divergence on this mesh is zero208 for n := 0; n < 4; n++ {209 var div utils.Matrix210 div = c.dfr.FluxElement.DivInt.Mul(F_RT_DOF[n])211 for k := 0; k < Kmax; k++ {212 for i := 0; i < Nint; i++ {213 ind := k + i*Kmax214 div.DataP[ind] /= c.dfr.Jdet.DataP[k]215 }216 }217 fmt.Printf("checking variable[%d]...", n)218 assert.True(t, nearVecScalar(div.DataP, 0., 0.000001))219 fmt.Printf("done.\n")220 }221 }222 }223 if false { // Test divergence of polynomial initial condition against analytic values224 /*225 Note: the Polynomial flux is asymmetric around the X and Y axes - it uses abs(x) and abs(y)226 Elements should not straddle the axes if a perfect polynomial flux capture is needed227 */228 Nmax := 7229 for N := 0; N <= Nmax; N++ {230 plotMesh := false231 // Single triangle test case232 var c *Euler233 //c = NewEuler(1, N, "../../DG2D/test_tris_1tri.neu", 1, FLUX_Average, FREESTREAM, 1, 0, 1.4, 0, false, 5000, None, plotMesh, false, false)234 ip.PolynomialOrder = N235 c = NewEuler(&ip, "../../DG2D/test_tris_1tri.neu", 1, plotMesh, false, false)236 CheckFlux0(c, t)237 // Two widely separated triangles - no shared faces238 //c = NewEuler(1, N, "../../DG2D/test_tris_two.neu", 1, FLUX_Average, FREESTREAM, 1, 0, 1.4, 0, false, 5000, None, plotMesh, false, false)239 c = NewEuler(&ip, "../../DG2D/test_tris_two.neu", 1, plotMesh, false, false)240 CheckFlux0(c, t)241 // Two widely separated triangles - no shared faces - one tri listed in reverse order242 //c = NewEuler(1, N, "../../DG2D/test_tris_twoR.neu", 1, FLUX_Average, FREESTREAM, 1, 0, 1.4, 0, false, 5000, None, plotMesh, false, false)243 c = NewEuler(&ip, "../../DG2D/test_tris_twoR.neu", 1, plotMesh, false, false)244 CheckFlux0(c, t)245 // Connected tris, sharing one edge246 // plotMesh = true247 //c = NewEuler(1, N, "../../DG2D/test_tris_6_nowall.neu", 1, FLUX_Average, FREESTREAM, 1, 0, 1.4, 0, false, 5000, None, plotMesh, false, false)248 c = NewEuler(&ip, "../../DG2D/test_tris_6_nowall.neu", 1, plotMesh, false, false)249 CheckFlux0(c, t)250 }251 }252 }253 if true {254 { // Test divergence of Isentropic Vortex initial condition against analytic values - density equation only255 N := 1256 plotMesh := false257 ip.PolynomialOrder = N258 ip.InitType = "ivortex"259 //c := NewEuler(1, N, "../../DG2D/test_tris_6.neu", 1, FLUX_Average, IVORTEX, 1, 0, 1.4, 0, false, 5000, None, plotMesh, false, false)260 c := NewEuler(&ip, "../../DG2D/test_tris_6.neu", 1, plotMesh, false, false)261 for _, e := range c.dfr.Tris.Edges {262 if e.BCType == types.BC_IVortex {263 e.BCType = types.BC_None264 }265 }266 Kmax := c.dfr.K267 Nint := c.dfr.FluxElement.NpInt268 Nedge := c.dfr.FluxElement.NpEdge269 NpFlux := c.dfr.FluxElement.Np // Np = 2*NpInt+3*NpEdge270 // Mark the initial state with the element number271 var Q_Face, F_RT_DOF [4]utils.Matrix272 var Flux, Flux_Face [2][4]utils.Matrix273 for n := 0; n < 4; n++ {274 F_RT_DOF[n] = utils.NewMatrix(NpFlux, Kmax)275 Q_Face[n] = utils.NewMatrix(3*Nedge, Kmax)276 Flux_Face[0][n] = utils.NewMatrix(3*Nedge, Kmax)277 Flux_Face[1][n] = utils.NewMatrix(3*Nedge, Kmax)278 Flux[0][n] = utils.NewMatrix(Nint, Kmax)279 Flux[1][n] = utils.NewMatrix(Nint, Kmax)280 }281 Q := c.Q[0]282 X, Y := c.dfr.FluxX, c.dfr.FluxY283 c.SetRTFluxInternal(Kmax, c.dfr.Jdet, c.dfr.Jinv, F_RT_DOF, Q)284 c.InterpolateSolutionToEdges(Q, Q_Face, Flux, Flux_Face)285 EdgeQ1 := make([][4]float64, Nedge)286 EdgeQ2 := make([][4]float64, Nedge)287 c.CalculateEdgeFlux(0, false, nil, nil, [][4]utils.Matrix{Q_Face},288 [][2][4]utils.Matrix{Flux_Face}, c.SortedEdgeKeys[0], EdgeQ1, EdgeQ2)289 c.SetRTFluxOnEdges(0, Kmax, F_RT_DOF)290 var div utils.Matrix291 // Density is the easiest equation to match with a polynomial292 n := 0293 //fmt.Printf("component[%d]\n", n)294 div = c.dfr.FluxElement.DivInt.Mul(F_RT_DOF[n])295 //c.DivideByJacobian(Kmax, NpInt, c.dfr.Jdet, div.DataP, 1)296 for k := 0; k < Kmax; k++ {297 for i := 0; i < Nint; i++ {298 ind := k + i*Kmax299 div.DataP[ind] /= -c.dfr.Jdet.DataP[k]300 }301 }302 // Get the analytic values of divergence for comparison303 for k := 0; k < Kmax; k++ {304 for i := 0; i < Nint; i++ {305 ind := k + i*Kmax306 x, y := X.DataP[ind], Y.DataP[ind]307 qc1, qc2, qc3, qc4 := c.AnalyticSolution.GetStateC(0, x, y)308 q1, q2, q3, q4 := Q[0].DataP[ind], Q[1].DataP[ind], Q[2].DataP[ind], Q[3].DataP[ind]309 assert.InDeltaSlicef(t, []float64{q1, q2, q3, q4}, []float64{qc1, qc2, qc3, qc4}, tol, msg)310 divC := c.AnalyticSolution.GetDivergence(0, x, y)311 divCalc := div.DataP[ind]312 assert.InDeltaf(t, divCalc/qc1, divC[n]/qc1, 0.001, msg) // 0.1 percent match313 }314 }315 }316 }317}318func TestFluxInterpolation(t *testing.T) {319 N := 2320 plotMesh := false321 ip := *ipDefault322 ip.Minf = 1.323 ip.PolynomialOrder = N324 ip.FluxType = "Roe"325 c := NewEuler(&ip, "../../DG2D/test_tris_6.neu", 1, plotMesh, false, false)326 rk := c.NewRungeKuttaSSP()327 c.InterpolateSolutionToEdges(c.Q[0], rk.Q_Face[0], rk.Flux[0], rk.Flux_Face[0])328 el := c.dfr.SolutionElement329 /*330 el.JB2D.V.Print("V")331 el.JB2D.Vinv.Print("Vinv")332 el.MassMatrix.Print("M")333 */334 RR := [][2]float64{{-1 / 3, -1 / 3}, {-1, -1}, {-1, 1}, {1, -1}}335 NpInt := el.Np336 locations := utils.NewMatrix(len(RR), NpInt)337 for ii, rr := range RR {338 var sm int339 for i := 0; i <= N; i++ {340 for j := 0; j <= N-i; j++ {341 //fmt.Printf("i,j = %d,%d\n", i, j)342 r, s := rr[0], rr[1]343 //fmt.Printf("Basis value at [%3.1f,%3.1f][%d,%d] = %5.3f\n", r, s, i, j, el.JB2D.PolynomialTerm(r, s, i, j))344 locations.Set(ii, sm, el.JB2D.PolynomialTerm(r, s, i, j))345 sm++346 }347 }348 }349 filter := make([]float64, 10)350 filter[0] = 1351 filter[1] = 1.0352 filter[2] = 1.0353 filter[3] = 0.00354 filter[4] = 0.00355 filter[5] = 0.00356 filter[6] = 0.00357 filter[7] = 0.00358 filter[8] = 0.00359 filter[9] = 0.00360 truncate := utils.NewDiagMatrix(NpInt, filter[:NpInt])361 //solution := utils.NewMatrix(NpInt, 1, []float64{1, 1, 1, 1, 1, 1})362 //modes := el.JB2D.Vinv.Mul(solution)363 //solution.Transpose().Print("solution-const")364 //modes.Transpose().Print("modes1")365 //locations.Mul(modes).Transpose().Print("Locations1")366 //solution = utils.NewMatrix(NpInt, 1, []float64{1, 1, 1, 0.5, 0.5, 1.5})367 solution := utils.NewMatrix(NpInt, 1, []float64{1.0, 0.5, 1.5, 1.5, 1.5, 1.5, 1, 1, 1, 1})368 modes := el.JB2D.Vinv.Mul(solution)369 solution.Transpose().Print("solution-variable")370 modes.Transpose().Print("modes2")371 locations.Mul(modes).Transpose().Print("Locations2")372 VinvFiltered := truncate.Mul(el.JB2D.Vinv)373 modes = VinvFiltered.Mul(solution)374 //modes.Transpose().Print("modes3")375 locations.Mul(modes).Transpose().Print("Locations3")376 R, S := utils.NewVector(4, []float64{-1 / 3, -1, -1, 1}), utils.NewVector(4, []float64{-1 / 3, -1, 1, -1})377 linterp := el.JB2D.GetInterpMatrix(R, S)378 linterp.Mul(solution).Transpose().Print("linterp-locations")379}380func TestFluxJacobian(t *testing.T) {381 var (382 tol = 0.000001383 msg = "err msg %s"384 )385 { // Flux Jacobian calculation386 c := Euler{}387 c.FSFar = NewFreeStream(0.1, 1.4, 0)388 Qinf := c.FSFar.Qinf389 Fx, Gy := c.FluxJacobianCalc(Qinf[0], Qinf[1], Qinf[2], Qinf[3])390 // Matlab: using the FSFar Q = [1,.1,0,1.79071]391 assert.InDeltaSlicef(t, []float64{392 0, 1.0000, 0, 0,393 -0.0080, 0.1600, 0, 0.4000,394 0, 0, 0.1000, 0,395 -0.2503, 2.5010, 0, 0.1400,396 }, Fx[:], tol, msg)397 assert.InDeltaSlicef(t, []float64{398 0, 0, 1.0000, 0,399 0, 0, 0.1000, 0,400 0.0020, -0.0400, 0, 0.4000,401 0, 0, 2.5050, 0,402 }, Gy[:], tol, msg)403 }404}405func TestEdges(t *testing.T) {406 dfr := DG2D.NewDFR2D(1, false, false, "../../DG2D/test_tris_9.neu")407 assert.Equal(t, len(dfr.Tris.Edges), 19)408 edges := make(EdgeKeySlice, len(dfr.Tris.Edges))409 var i int410 for key := range dfr.Tris.Edges {411 edges[i] = key412 i++413 }414 edges.Sort()415 //fmt.Printf("len(Edges) = %d, Edges = %v\n", len(edges), edges)416 l := make([]int, len(edges))417 for i, e := range edges {418 //fmt.Printf("vertex[edge[%d]]=%d\n", i, e.GetVertices(false)[1])419 l[i] = e.GetVertices(false)[1]420 }421 assert.Equal(t, []int{1, 2, 3, 4, 4, 4, 5, 5, 5, 6, 6, 7, 7, 8, 8, 8, 9, 9, 9}, l)422 edges2 := EdgeKeySliceSortLeft(edges)423 edges2.Sort()424 for i, e := range edges2 {425 //v := e.GetVertices(false)426 //fmt.Printf("vertex2[edge[%d]]=[%d,%d]\n", i, v[0], v[1])427 l[i] = e.GetVertices(false)[0]428 }429 assert.Equal(t, []int{0, 0, 0, 1, 1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 6, 7, 8}, l)430}431func TestDissipation(t *testing.T) {432 {433 dfr := DG2D.NewDFR2D(1, false, false, "../../DG2D/test_tris_9.neu")434 VtoE := NewVertexToElement(dfr.Tris.EToV)435 assert.Equal(t, VertexToElement{{0, 0, 0}, {0, 1, 0}, {1, 3, 0}, {1, 1, 0}, {1, 2, 0}, {2, 4, 0}, {2, 3, 0}, {3, 5, 0}, {3, 0, 0}, {4, 2, 0}, {4, 5, 0}, {4, 6, 0}, {4, 1, 0}, {4, 0, 0}, {4, 7, 0}, {5, 4, 0}, {5, 3, 0}, {5, 9, 0}, {5, 8, 0}, {5, 2, 0}, {5, 7, 0}, {6, 4, 0}, {6, 9, 0}, {7, 5, 0}, {7, 6, 0}, {8, 8, 0}, {8, 6, 0}, {8, 7, 0}, {9, 8, 0}, {9, 9, 0}},436 VtoE)437 vepFinal := [3]int32{9, 9, 0}438 for NPar := 1; NPar < 10; NPar += 2 {439 pm := NewPartitionMap(NPar, dfr.K)440 sd := NewScalarDissipation(0, dfr, pm)441 var vep [3]int32442 for np := 0; np < NPar; np++ {443 for _, val := range sd.VtoE[np] {444 vep = val445 }446 }447 assert.Equal(t, vepFinal[0], vep[0])448 }449 KMax := dfr.K450 assert.Equal(t, len(dfr.Jdet.DataP), KMax)451 for k := 0; k < KMax; k++ {452 area := 2. * dfr.Jdet.DataP[k] // Area of element is 2x Determinant, because the unit triangle is area=2453 assert.InDeltaf(t, area, 0.25, 0.000001, "err msg %s")454 }455 }456 if false { // Turn off value check tests while working on the constants in the artificial dissipation457 dfr := DG2D.NewDFR2D(2, false, false, "../../DG2D/test_tris_9.neu")458 Np, KMax := dfr.SolutionElement.Np, dfr.K459 pm := NewPartitionMap(1, KMax)460 Q := make([][4]utils.Matrix, 1)461 n := 0462 Q[0][n] = utils.NewMatrix(dfr.SolutionElement.Np, KMax)463 var val float64464 for k := 0; k < KMax; k++ {465 for i := 0; i < Np; i++ {466 if i < Np/2 {467 val = 2468 } else {469 val = 1470 }471 ind := k + KMax*i472 Q[0][n].DataP[ind] = val473 }474 }475 sd := NewScalarDissipation(0, dfr, pm)476 sd.Kappa = 4.477 sd.CalculateElementViscosity(0, Q)478 //assert.InDeltaSlicef(t, []float64{0.09903, 0.09903, 0.09903, 0.09903, 0.09903, 0.09903, 0.09903, 0.09903, 0.09903, 0.09903},479 assert.InDeltaSlicef(t, []float64{0.09903, 0.09903, 0.09903, 0.07003, 0.09903, 0.07003, 0.09903, 0.07003, 0.09903, 0.07003},480 sd.EpsilonScalar[0], 0.00001, "err msg %s")481 sd.Kappa = 0.75482 sd.CalculateElementViscosity(0, Q)483 //assert.InDeltaSlicef(t, []float64{0.01270, 0.01270, 0.01270, 0.01270, 0.01270, 0.01270, 0.01270, 0.01270, 0.01270, 0.01270},484 assert.InDeltaSlicef(t, []float64{0.01270, 0.01270, 0.01270, 0.00898, 0.01270, 0.00898, 0.01270, 0.00898, 0.01270, 0.00898},485 sd.EpsilonScalar[0], 0.00001, "err msg %s")486 }487 {488 dfr := DG2D.NewDFR2D(1, false, false, "../../DG2D/test_tris_9.neu")489 _, KMax := dfr.SolutionElement.Np, dfr.K490 for NP := 1; NP < 5; NP++ {491 pm := NewPartitionMap(NP, KMax)492 sd := NewScalarDissipation(0, dfr, pm)493 /*494 assert.InDeltaSlicef(t, []float64{0.666666, 0.166666, 0.166666, 0.166666, 0.666666, 0.166666, 0.166666, 0.166666,495 0.666666, 0.666666, 0.166666, 0.166666, 0.166666, 0.666666, 0.166666, 0.166666, 0.166666, 0.666666, 0.827326,496 0.172673, 0, 0.5, 0.5, 0, 0.172673, 0.827326, 0, 0, 0.827326, 0.172673, 0, 0.5, 0.5, 0, 0.172673, 0.827326,497 0.172673, 0, 0.827326, 0.5, 0, 0.5, 0.827326, 0, 0.172673},498 sd.BaryCentricCoords.DataP, 0.00001, "err msg %s")499 */500 // Set the epsilon scalar value to the element ID and check the vertex aggregation501 for np := 0; np < NP; np++ {502 KMaxLocal := pm.GetBucketDimension(np)503 for k := 0; k < KMaxLocal; k++ {504 kGlobal := pm.GetGlobalK(k, np)505 sd.EpsilonScalar[np][k] = float64(kGlobal)506 }507 }508 wg := &sync.WaitGroup{}509 for np := 0; np < NP; np++ {510 wg.Add(1)511 sd.propagateEpsilonMaxToVertices(np)512 }513 assert.Equal(t, []float64{1, 3, 4, 5, 7, 9, 9, 6, 8, 9}, sd.EpsVertex)514 }515 }516}517func TestDissipation2(t *testing.T) {518 // Test C0 continuity of Epsilon using element vertex aggregation519 {520 var (521 dfr = DG2D.NewDFR2D(1, false, false, "../../DG2D/test_tris_9.neu")522 )523 NP := 1524 _, KMax := dfr.SolutionElement.Np, dfr.K525 pm := NewPartitionMap(NP, KMax)526 sd := NewScalarDissipation(0, dfr, pm)527 for np := 0; np < NP; np++ {528 KMax = sd.PMap.GetBucketDimension(np)529 for k := 0; k < KMax; k++ {530 kGlobal := pm.GetGlobalK(k, np)531 sd.EpsilonScalar[np][k] = float64(kGlobal)532 }533 }534 wg := &sync.WaitGroup{}535 for np := 0; np < NP; np++ {536 wg.Add(1)537 sd.propagateEpsilonMaxToVertices(np)538 }539 assert.Equal(t, []float64{1, 3, 4, 5, 7, 9, 9, 6, 8, 9}, sd.EpsVertex)540 for np := 0; np < NP; np++ {541 sd.linearInterpolateEpsilon(np)542 //sd.baryCentricInterpolateEpsilon(np)543 }544 /*545 assert.InDeltaSlicef(t, []float64{546 5.66667, 3.33333, 7.66667, 7.16667, 8.16667, 6.00000, 7.50000, 8.00000, 8.83333, 9.00000, 4.66667, 5.33333,547 6.66667, 4.16667, 8.16667, 5.50000, 6.50000, 7.50000, 8.33333, 9.00000, 2.66667, 2.33333, 4.66667, 4.66667,548 5.66667, 6.50000, 7.00000, 8.50000, 8.83333, 9.00000, 5.66667, 3.33333, 7.66667, 7.16667, 8.16667, 6.00000,549 7.50000, 8.00000, 8.83333, 9.00000, 4.66667, 5.33333, 6.66667, 4.16667, 8.16667, 5.50000, 6.50000, 7.50000,550 8.33333, 9.00000, 2.66667, 2.33333, 4.66667, 4.66667, 5.66667, 6.50000, 7.00000, 8.50000, 8.83333, 9.00000,551 6.65465, 3.69069, 8.65465, 7.96396, 9.00000, 5.82733, 7.65465, 7.82733, 8.82733, 9.00000, 6.00000, 5.00000,552 8.00000, 6.00000, 9.00000, 5.50000, 7.00000, 7.50000, 8.50000, 9.00000, 5.34535, 6.30931, 7.34535, 4.03604,553 9.00000, 5.17267, 6.34535, 7.17267, 8.17267, 9.00000, 4.30931, 5.96396, 6.30931, 3.17267, 8.13663, 5.34535,554 6.17267, 7.34535, 8.17267, 9.00000, 3.00000, 4.00000, 5.00000, 3.50000, 6.50000, 6.00000, 6.50000, 8.00000,555 8.50000, 9.00000, 1.69069, 2.03604, 3.69069, 3.82733, 4.86337, 6.65465, 6.82733, 8.65465, 8.82733, 9.00000,556 2.03604, 1.34535, 4.03604, 4.86337, 4.86337, 6.82733, 7.17267, 8.82733, 9.00000, 9.00000, 4.00000, 2.00000,557 6.00000, 6.50000, 6.50000, 6.50000, 7.50000, 8.50000, 9.00000, 9.00000, 5.96396, 2.65465, 7.96396, 8.13663,558 8.13663, 6.17267, 7.82733, 8.17267, 9.00000, 9.00000},559 sd.Epsilon[0].DataP, 0.00001, "err msg %s")560 */561 //fmt.Printf(sd.Epsilon[0].Print("Epsilon"))562 }563 // Gradient test using GetSolutionGradient()564 {565 ip := *ipDefault566 ip.FluxType = "average"567 // Testing to fourth order in X and Y568 ip.PolynomialOrder = 4569 c := NewEuler(&ip, "../../DG2D/test_tris_5.neu", 1, false, false, false)570 var (571 dfr = c.dfr572 Kmax = dfr.K573 fel = dfr.FluxElement574 NpInt = fel.NpInt575 Nedge = fel.NpEdge576 NpFlux = fel.Np577 Q, Q_Face [4]utils.Matrix578 QGradXCheck, QGradYCheck [4]utils.Matrix579 )580 for n := 0; n < 4; n++ {581 Q[n] = utils.NewMatrix(NpInt, Kmax)582 Q_Face[n] = utils.NewMatrix(3*Nedge, Kmax)583 QGradXCheck[n] = utils.NewMatrix(NpFlux, Kmax)584 QGradYCheck[n] = utils.NewMatrix(NpFlux, Kmax)585 }586 // Each variable [n] is an nth order polynomial in X and Y587 for i := 0; i < NpFlux; i++ {588 for k := 0; k < Kmax; k++ {589 ind := k + i*Kmax590 X, Y := dfr.FluxX.DataP[ind], dfr.FluxY.DataP[ind]591 if i < NpInt {592 Q[0].DataP[ind] = X + Y593 Q[1].DataP[ind] = X*X + Y*Y594 Q[2].DataP[ind] = X*X*X + Y*Y*Y595 Q[3].DataP[ind] = X*X*X*X + Y*Y*Y*Y596 }597 // First order in X and Y598 QGradXCheck[0].DataP[ind] = 1.599 QGradYCheck[0].DataP[ind] = 1.600 // Second order in X and Y601 QGradXCheck[1].DataP[ind] = 2 * X602 QGradYCheck[1].DataP[ind] = 2 * Y603 // Third order in X and Y604 QGradXCheck[2].DataP[ind] = 3 * X * X605 QGradYCheck[2].DataP[ind] = 3 * Y * Y606 // Fourth order in X and Y607 QGradXCheck[3].DataP[ind] = 4 * X * X * X608 QGradYCheck[3].DataP[ind] = 4 * Y * Y * Y609 }610 }611 // Interpolate from solution points to edges using precomputed interpolation matrix612 fI := dfr.FluxEdgeInterp.Mul613 for n := 0; n < 4; n++ {614 Q_Face[n] = fI(Q[n])615 }616 // Test hard coded loop GradX and GradY617 {618 for n := 0; n < 4; n++ {619 var (620 DXmd, DYmd = dfr.DXMetric.DataP, dfr.DYMetric.DataP621 DOFX, DOFY = utils.NewMatrix(NpFlux, Kmax), utils.NewMatrix(NpFlux, Kmax)622 DOFXd, DOFYd = DOFX.DataP, DOFY.DataP623 )624 var Un float64625 for k := 0; k < Kmax; k++ {626 for i := 0; i < NpFlux; i++ {627 ind := k + i*Kmax628 switch {629 case i < NpInt: // The first NpInt points are the solution element nodes630 Un = Q[n].DataP[ind]631 case i >= NpInt && i < 2*NpInt: // The second NpInt points are duplicates of the first NpInt values632 Un = Q[n].DataP[ind-NpInt*Kmax]633 case i >= 2*NpInt:634 Un = Q_Face[n].DataP[ind-2*NpInt*Kmax] // The last 3*NpEdge points are the edges in [0-1,1-2,2-0] order635 }636 DOFXd[ind] = DXmd[ind] * Un637 DOFYd[ind] = DYmd[ind] * Un638 }639 }640 DX := dfr.FluxElement.Div.Mul(DOFX) // X Derivative, Divergence x RT_DOF is X derivative for this DOF641 DY := dfr.FluxElement.Div.Mul(DOFY) // Y Derivative, Divergence x RT_DOF is Y derivative for this DOF642 fmt.Printf("Order[%d] check ...", n+1)643 assert.Equal(t, len(DX.DataP), len(QGradXCheck[n].DataP))644 assert.Equal(t, len(DY.DataP), len(QGradYCheck[n].DataP))645 assert.InDeltaSlicef(t, DX.DataP, QGradXCheck[n].DataP, 0.000001, "err msg %s")646 assert.InDeltaSlicef(t, DY.DataP, QGradYCheck[n].DataP, 0.000001, "err msg %s")647 fmt.Printf("... validates\n")648 }649 }650 // Test function for GradX and GradY651 {652 var (653 DX, DY = utils.NewMatrix(NpFlux, Kmax), utils.NewMatrix(NpFlux, Kmax)654 DOFX, DOFY = utils.NewMatrix(NpFlux, Kmax), utils.NewMatrix(NpFlux, Kmax)655 )656 // Before this call, we need to load edge data into the edge store657 edgeValues := make([][4]float64, Nedge)658 myThread := -1659 for k := 0; k < Kmax; k++ {660 kGlobal := c.Partitions.GetGlobalK(k, myThread)661 for edgeNum := 0; edgeNum < 3; edgeNum++ {662 en := dfr.EdgeNumber[kGlobal+Kmax*edgeNum]663 primeElement := c.dfr.Tris.Edges[en].ConnectedTris[0]664 if int(primeElement) == kGlobal {665 for i := 0; i < Nedge; i++ {666 ind := k + (i+edgeNum*Nedge)*Kmax667 for n := 0; n < 4; n++ {668 edgeValues[i][n] = Q_Face[n].DataP[ind]669 }670 }671 c.EdgeStore.PutEdgeValues(en, QFluxForGradient, edgeValues)672 }673 }674 }675 for n := 0; n < 4; n++ {676 fmt.Printf("Order[%d] check ...", n+1)677 c.GetSolutionGradient(-1, n, Q, DX, DY, DOFX, DOFY)678 assert.Equal(t, len(DX.DataP), len(QGradXCheck[n].DataP))679 assert.Equal(t, len(DY.DataP), len(QGradYCheck[n].DataP))680 assert.InDeltaSlicef(t, DX.DataP, QGradXCheck[n].DataP, 0.000001, "err msg %s")681 assert.InDeltaSlicef(t, DY.DataP, QGradYCheck[n].DataP, 0.000001, "err msg %s")682 fmt.Printf("... validates\n")683 }684 }685 }686}687func TestEuler_GetSolutionGradientUsingRTElement(t *testing.T) {688 {689 ip := *ipDefault690 ip.FluxType = "average"691 // Testing to fourth order in X and Y692 ip.PolynomialOrder = 4693 ip.Minf = 0.8694 ip.Alpha = 2.695 c := NewEuler(&ip, "../../DG2D/test_tris_9.neu", 1,696 false, false, false)697 rk := c.NewRungeKuttaSSP()698 myThread := 0699 var (700 dfr = c.dfr701 Kmax = dfr.K702 fel = dfr.FluxElement703 Nedge = fel.NpEdge704 NpFlux = fel.Np705 QGradXCheck, QGradYCheck [4]utils.Matrix706 Q0 = c.Q[myThread]707 SortedEdgeKeys = c.SortedEdgeKeys[myThread]708 EdgeQ1 = make([][4]float64, Nedge) // Local working memory709 EdgeQ2 = make([][4]float64, Nedge) // Local working memory710 )711 for n := 0; n < 4; n++ {712 QGradXCheck[n] = utils.NewMatrix(NpFlux, Kmax)713 QGradYCheck[n] = utils.NewMatrix(NpFlux, Kmax)714 }715 // Flow is freestream, graident should be 0 in all variables716 for n := 0; n < 4; n++ {717 for i := 0; i < NpFlux; i++ {718 for k := 0; k < Kmax; k++ {719 ind := k + i*Kmax720 QGradXCheck[n].DataP[ind] = 0.721 QGradYCheck[n].DataP[ind] = 0.722 }723 }724 }725 // Interpolate from solution points to edges using precomputed interpolation matrix726 // Test function for GradX and GradY727 {728 var (729 DX, DY = utils.NewMatrix(NpFlux, Kmax), utils.NewMatrix(NpFlux, Kmax)730 DOFX, DOFY = utils.NewMatrix(NpFlux, Kmax), utils.NewMatrix(NpFlux, Kmax)731 )732 rk.MaxWaveSpeed[0] =733 c.CalculateEdgeFlux(rk.Time, true, rk.Jdet, rk.DT, rk.Q_Face, rk.Flux_Face, SortedEdgeKeys,734 EdgeQ1, EdgeQ2) // Global735 for n := 0; n < 4; n++ {736 fmt.Printf("Variable[%d] check ...", n+1)737 c.GetSolutionGradient(-1, n, Q0, DX, DY, DOFX, DOFY)738 assert.Equal(t, len(DX.DataP), len(QGradXCheck[n].DataP))739 assert.Equal(t, len(DY.DataP), len(QGradYCheck[n].DataP))740 assert.InDeltaSlicef(t, DX.DataP, QGradXCheck[n].DataP, 0.000001, "err msg %s")741 assert.InDeltaSlicef(t, DY.DataP, QGradYCheck[n].DataP, 0.000001, "err msg %s")742 fmt.Printf("... validates\n")743 }744 }745 }746}747func TestInputParameters_Parse(t *testing.T) {748 var (749 err error750 )751 fileInput := []byte(`752Title: Test Case753CFL: 1.754InitType: Freestream # Can be IVortex or Freestream755FluxType: Roe756PolynomialOrder: 2757FinalTime: 4.758BCs: 759 Inflow:760 37:761 NPR: 4.0762 Outflow:763 22:764 P: 1.5765`)766 var input InputParameters767 if err = input.Parse(fileInput); err != nil {768 panic(err)769 }770 // Check Inflow BC number 37771 assert.Equal(t, input.BCs["Inflow"][37]["NPR"], 4.)772 // Check Outflow BC number 22773 assert.Equal(t, input.BCs["Outflow"][22]["P"], 1.5)774 input.Print()775 assert.Equal(t, input.FinalTime, 4.)776}777func PrintQ(Q [4]utils.Matrix, l string) {778 var (779 label string780 )781 for ii := 0; ii < 4; ii++ {782 switch ii {783 case 0:784 label = l + "_0"785 case 1:786 label = l + "_1"787 case 2:788 label = l + "_2"789 case 3:790 label = l + "_3"791 }792 fmt.Println(Q[ii].Transpose().Print(label))793 }794}795func PrintFlux(F []utils.Matrix) {796 for ii := 0; ii < len(F); ii++ {797 label := strconv.Itoa(ii)798 fmt.Println(F[ii].Print("F" + "[" + label + "]"))799 }800}801func nearVecScalar(a []float64, b float64, tol float64) (l bool) {802 near := func(a, b float64, tolI ...float64) (l bool) {803 var (804 tol float64805 )806 if len(tolI) == 0 {807 tol = 1.e-08808 } else {809 tol = tolI[0]810 }811 bound := math.Max(tol, tol*math.Abs(a))812 val := math.Abs(a - b)813 if val <= bound {814 l = true815 } else {816 fmt.Printf("Diff = %v, Left = %v, Right = %v\n", val, a, b)817 }818 return819 }820 for i, val := range a {821 if !near(b, val, tol) {822 fmt.Printf("Diff = %v, Left[%d] = %v, Right[%d] = %v\n", math.Abs(val-b), i, val, i, b)823 return false824 }825 }826 return true827}828func InitializePolynomial(X, Y utils.Matrix) (Q [4]utils.Matrix) {829 var (830 Np, Kmax = X.Dims()831 )832 for n := 0; n < 4; n++ {833 Q[n] = utils.NewMatrix(Np, Kmax)834 }835 for ii := 0; ii < Np*Kmax; ii++ {836 x, y := X.DataP[ii], Y.DataP[ii]837 rho, rhoU, rhoV, E := GetStatePoly(x, y)838 Q[0].DataP[ii] = rho839 Q[1].DataP[ii] = rhoU840 Q[2].DataP[ii] = rhoV841 Q[3].DataP[ii] = E842 }843 return844}845func GetStatePoly(x, y float64) (rho, rhoU, rhoV, E float64) {846 /*847 Matlab script:848 syms a b c d x y gamma849 %2D Polynomial field850 rho=a*abs(x)+b*abs(y);851 u = c*x; v = d*y;852 rhou=rho*u; rhov=rho*v;853 p=rho^gamma;854 q=0.5*rho*(u^2+v^2);855 E=p/(gamma-1)+q;856 U = [ rho, rhou, rhov, E];857 F = [ rhou, rho*u^2+p, rho*u*v, u*(E+p) ];858 G = [ rhov, rho*u*v, rho*v^2+p, v*(E+p) ];859 div = diff(F,x)+diff(G,y);860 fprintf('Code for Divergence of F and G FluxIndex\n%s\n',ccode(div));861 fprintf('Code for U \n%s\n%s\n%s\n%s\n',ccode(U));862 */863 var (864 a, b, c, d = 1., 1., 1., 1.865 pow = math.Pow866 fabs = math.Abs867 gamma = 1.4868 )869 rho = a*fabs(x) + b*fabs(y)870 rhoU = c * x * (a*fabs(x) + b*fabs(y))871 rhoV = d * y * (a*fabs(x) + b*fabs(y))872 E = ((c*c)*(x*x)+(d*d)*(y*y))*((a*fabs(x))/2.0+(b*fabs(y))/2.0) + pow(a*fabs(x)+b*fabs(y), gamma)/(gamma-1.0)873 return874}875func GetDivergencePoly(t, x, y float64) (div [4]float64) {876 var (877 gamma = 1.4878 pow = math.Pow879 fabs = math.Abs880 a, b, c, d = 1., 1., 1., 1.881 )882 div[0] = c*(a*fabs(x)+b*fabs(y)) + d*(a*fabs(x)+b*fabs(y)) + a*c*x*(x/fabs(x)) + b*d*y*(y/fabs(y))883 div[1] = (c*c)*x*(a*fabs(x)+b*fabs(y))*2.0 + c*d*x*(a*fabs(x)+b*fabs(y)) + a*(c*c)*(x*x)*(x/fabs(x)) + a*gamma*(x/fabs(x))*pow(a*fabs(x)+b*fabs(y), gamma-1.0) + b*c*d*x*y*(y/fabs(y))884 div[2] = (d*d)*y*(a*fabs(x)+b*fabs(y))*2.0 + c*d*y*(a*fabs(x)+b*fabs(y)) + b*(d*d)*(y*y)*(y/fabs(y)) + b*gamma*(y/fabs(y))*pow(a*fabs(x)+b*fabs(y), gamma-1.0) + a*c*d*x*y*(x/fabs(x))885 div[3] = c*(((c*c)*(x*x)+(d*d)*(y*y))*((a*fabs(x))/2.0+(b*fabs(y))/2.0)+pow(a*fabs(x)+b*fabs(y), gamma)+pow(a*fabs(x)+b*fabs(y), gamma)/(gamma-1.0)) + d*(((c*c)*(x*x)+(d*d)*(y*y))*((a*fabs(x))/2.0+(b*fabs(y))/2.0)+pow(a*fabs(x)+b*fabs(y), gamma)+pow(a*fabs(x)+b*fabs(y), gamma)/(gamma-1.0)) + c*x*((c*c)*x*((a*fabs(x))/2.0+(b*fabs(y))/2.0)*2.0+(a*(x/fabs(x))*((c*c)*(x*x)+(d*d)*(y*y)))/2.0+a*gamma*(x/fabs(x))*pow(a*fabs(x)+b*fabs(y), gamma-1.0)+(a*gamma*(x/fabs(x))*pow(a*fabs(x)+b*fabs(y), gamma-1.0))/(gamma-1.0)) + d*y*((b*(y/fabs(y))*((c*c)*(x*x)+(d*d)*(y*y)))/2.0+(d*d)*y*((a*fabs(x))/2.0+(b*fabs(y))/2.0)*2.0+b*gamma*(y/fabs(y))*pow(a*fabs(x)+b*fabs(y), gamma-1.0)+(b*gamma*(y/fabs(y))*pow(a*fabs(x)+b*fabs(y), gamma-1.0))/(gamma-1.0))886 return887}888func FluxCalcMomentumOnly(rho, rhoU, rhoV, E float64) (Fx, Fy [4]float64) {889 Fx, Fy =890 [4]float64{rhoU, rhoU, rhoU, rhoU},891 [4]float64{rhoV, rhoV, rhoV, rhoV}892 return893}894func CheckFlux0(c *Euler, t *testing.T) {895 /*896 Conditions of this test:897 - Two duplicated triangles, removes the question of transformation Jacobian making the results differ898 - Flux is calculated identically for each equation (only density components), removes the question of flux899 accuracy being different for the more complex equations900 - Flowfield is initialized to a freestream for a polynomial field, interpolation to edges is not done,901 instead, analytic initialization values are put into the edges902 Result:903 - No test of different triangle shapes and orientations904 - No test of accuracy of interpolation to edges905 - No accuracy test of the complex polynomial fluxes in Q[1-3]906 */907 if c.Partitions.ParallelDegree != 1 {908 panic("parallel degree should be 1 for this test")909 }910 c.FluxCalcMock = FluxCalcMomentumOnly // For testing, only consider the first component of flux for all [4]911 // Initialize912 X, Y := c.dfr.FluxX, c.dfr.FluxY913 QFlux := InitializePolynomial(X, Y)914 Kmax := c.dfr.K915 Nint := c.dfr.FluxElement.NpInt916 Nedge := c.dfr.FluxElement.NpEdge917 NpFlux := c.dfr.FluxElement.Np918 var Q, Q_Face, F_RT_DOF [4]utils.Matrix919 var Flux_Face [2][4]utils.Matrix920 for n := 0; n < 4; n++ {921 Q[n] = utils.NewMatrix(Nint, Kmax)922 Q_Face[n] = utils.NewMatrix(3*Nedge, Kmax)923 Flux_Face[0][n] = utils.NewMatrix(3*Nedge, Kmax)924 Flux_Face[1][n] = utils.NewMatrix(3*Nedge, Kmax)925 F_RT_DOF[n] = utils.NewMatrix(NpFlux, Kmax)926 for k := 0; k < Kmax; k++ {927 for i := 0; i < Nint; i++ {928 ind := k + i*Kmax929 Q[n].DataP[ind] = QFlux[n].DataP[ind]930 }931 for i := 0; i < 3*Nedge; i++ {932 ind := k + i*Kmax933 ind2 := k + (i+2*Nint)*Kmax934 Q_Face[n].DataP[ind] = QFlux[n].DataP[ind2]935 }936 }937 }938 c.SetRTFluxInternal(Kmax, c.dfr.Jdet, c.dfr.Jinv, F_RT_DOF, Q)939 EdgeQ1 := make([][4]float64, Nedge)940 EdgeQ2 := make([][4]float64, Nedge)941 // No need to interpolate to the edges, they are left at initialized state in Q_Face942 c.CalculateEdgeFlux(0, false, nil, nil, [][4]utils.Matrix{Q_Face},943 [][2][4]utils.Matrix{Flux_Face}, c.SortedEdgeKeys[0], EdgeQ1, EdgeQ2)944 c.SetRTFluxOnEdges(0, Kmax, F_RT_DOF)945 var div utils.Matrix946 for n := 0; n < 4; n++ {947 div = c.dfr.FluxElement.DivInt.Mul(F_RT_DOF[n])948 d1, d2 := div.Dims()949 assert.Equal(t, d1, Nint)950 assert.Equal(t, d2, Kmax)951 for k := 0; k < Kmax; k++ {952 Jdet := c.dfr.Jdet.At(k, 0)953 for i := 0; i < Nint; i++ {954 ind := k + i*Kmax955 div.DataP[ind] /= Jdet956 }957 }958 // Get the analytic values of divergence for comparison959 nn := 0 // Use only the density component of divergence to check960 for k := 0; k < Kmax; k++ {961 for i := 0; i < Nint; i++ {962 ind := k + i*Kmax963 x, y := X.DataP[ind], Y.DataP[ind]964 divC := GetDivergencePoly(0, x, y)965 divCalc := div.DataP[ind]966 normalizer := Q[nn].DataP[ind]967 //test := near(divCalc/normalizer, divC[nn]/normalizer, 0.0001) // 1% of field value968 assert.InDeltaf(t, divCalc/normalizer, divC[nn]/normalizer, 0.0001, "err msg %s") // 1% of field value969 }970 }971 }972}...

Full Screen

Full Screen

leds.go

Source:leds.go Github

copy

Full Screen

1package main2import (3 "fmt"4 "io"5 "math"6 "strings"7 "time"8 "github.com/ghetzel/go-stockutil/colorutil"9 "github.com/ghetzel/go-stockutil/log"10 "github.com/ghetzel/go-stockutil/stringutil"11 "github.com/ghetzel/go-stockutil/typeutil"12)13var DefaultTransitionShaderDuration = 3000 * time.Millisecond14var DefaultColor = colorutil.MustParse(`#00000000`)15type LEDSet map[int]colorutil.Color16func (self LEDSet) Has(i int) bool {17 if len(self) == 0 || i == math.MaxInt {18 return true19 }20 if _, ok := self[i]; ok {21 return true22 } else {23 return false24 }25}26func (self LEDSet) Get(i int) (colorutil.Color, bool) {27 if c, ok := self[i]; ok && !c.IsZero() {28 return c, true29 } else if c, ok := self[math.MaxInt]; ok && !c.IsZero() {30 return c, true31 } else {32 return colorutil.MustParse(`black`), false33 }34}35func ParseLEDRange(rangespec string) (leds LEDSet) {36 leds = make(LEDSet)37 rangespec = strings.TrimSpace(rangespec)38 for _, subrange := range strings.Split(rangespec, `,`) {39 var color colorutil.Color40 var index, colorspec = stringutil.SplitPairTrimSpace(subrange, `@`)41 if colorspec == `` {42 colorspec = subrange43 index = `*`44 }45 switch colorspec {46 case ``:47 continue48 case `-`:49 color = colorutil.MustParse(`rgba(0,0,0,1)`)50 default:51 color = colorutil.MustParse(colorspec)52 }53 if index == `*` {54 leds[math.MaxInt] = color55 return56 } else if a, b := stringutil.SplitPairTrimSpace(index, `:`); a != `` {57 var ai int = typeutil.NInt(a)58 if b != `` {59 var bi int = typeutil.NInt(b)60 for i := ai; i < bi; i++ {61 leds[i] = color62 }63 } else {64 leds[ai] = color65 }66 }67 }68 return69}70type Protocol byte71const (72 wled_WARLS Protocol = 0x173 wled_DRGB = 0x274 wled_DRGBW = 0x375 wled_DNRGB = 0x476 wled_NOTIFY = 0x077)78func (self LEDSet) Bytes(proto Protocol, timeout uint8) []byte {79 var payload = make([]byte, len(self)*4)80 for i := 0; i < len(self); i++ {81 if c, ok := self[i]; ok {82 var r, g, b, _ uint8 = c.RGBA255()83 payload[0+(i*4)] = byte(i)84 payload[1+(i*4)] = r85 payload[2+(i*4)] = g86 payload[3+(i*4)] = b87 }88 }89 return append([]byte{90 byte(proto),91 byte(timeout),92 }, payload...)93}94type Display struct {95 FrontBuffer LEDSet96 BackBuffer LEDSet97 FrameInterval time.Duration98 AutoclearTimeout uint899 TransitionShader ShaderFunc100 TransitionShaderDuration time.Duration101 TransitionArgs []string102 DefaultColor colorutil.Color103 ClearFirst bool104 ledcount int105 wledDest io.Writer106 proto Protocol107}108func NewDisplay(w io.Writer, ledcount int) *Display {109 var sset = &Display{110 FrontBuffer: make(LEDSet),111 BackBuffer: make(LEDSet),112 FrameInterval: 16 * time.Millisecond,113 AutoclearTimeout: 255,114 TransitionShader: nil,115 TransitionShaderDuration: DefaultTransitionShaderDuration,116 TransitionArgs: make([]string, 0),117 DefaultColor: DefaultColor,118 ledcount: ledcount,119 wledDest: w,120 proto: wled_WARLS,121 }122 sset.init()123 return sset124}125func (self *Display) SetTransitionEffect(effect string, args ...string) error {126 effect = strings.ToLower(effect)127 self.TransitionArgs = args128 switch effect {129 case ``, `fill`:130 self.TransitionShader = Fill131 self.TransitionShaderDuration = 0132 case `fade`:133 self.TransitionShader = Fade134 default:135 return fmt.Errorf("unknown effect %q", effect)136 }137 return nil138}139func (self *Display) init() {140 for i := 0; i < self.ledcount; i++ {141 self.FrontBuffer[i] = self.DefaultColor142 self.BackBuffer[i] = self.DefaultColor143 }144}145func (self *Display) flushBuffer(buffer LEDSet) error {146 var buf = buffer.Bytes(147 self.proto,148 self.AutoclearTimeout,149 )150 // log.Debugf("flush: %x", buf)151 var _, err = self.wledDest.Write(buf)152 return err153}154func (self *Display) Flush() error {155 var fi time.Duration = self.FrameInterval156 var td time.Duration = self.TransitionShaderDuration157 log.Debugf("flush requested: fn=%p transition=%v", self.TransitionShader, td)158 self.init()159 if tfn := self.TransitionShader; tfn != nil {160 var progressStep float64 = float64(fi) / float64(td)161 var frame int = 1162 var progress float64 = 0163 var start = time.Now()164 log.Debugf("anim start: step=%f duration=%v", progressStep, td)165 for progress <= 1.0 && !math.IsInf(progress, 0) {166 var transitionBuffer = make(LEDSet)167 for i := 0; i < self.ledcount; i++ {168 transitionBuffer[i] = self.BackBuffer[i]169 }170 var continueAnimation bool171 for i := 0; i < self.ledcount; i++ {172 var bc = self.BackBuffer[i]173 var fc = self.FrontBuffer[i]174 var lc = bc175 var tc = lc176 if i > 0 {177 lc = transitionBuffer[i-1]178 }179 tc, continueAnimation = tfn(&ShaderStep{180 Index: i,181 Frame: frame,182 AnimationProgress: progress,183 CurrentColor: fc,184 LastStepColor: lc,185 Args: self.TransitionArgs,186 })187 transitionBuffer[i] = tc188 self.BackBuffer[i] = tc189 }190 self.flushBuffer(transitionBuffer)191 if !continueAnimation {192 break193 }194 if fi > 0 {195 time.Sleep(fi)196 }197 frame += 1198 progress += progressStep199 }200 log.Debugf("anim done: frames=%d progress=%v took=%v", frame, progress, time.Since(start))201 }202 for i := 0; i < len(self.FrontBuffer); i++ {203 self.FrontBuffer[i] = self.BackBuffer[i]204 }205 return self.flushBuffer(self.FrontBuffer)206}...

Full Screen

Full Screen

serieslyinfo.go

Source:serieslyinfo.go Github

copy

Full Screen

1package main2import (3 "encoding/json"4 "flag"5 "fmt"6 "html/template"7 "io/ioutil"8 "log"9 "os"10 "text/tabwriter"11 "github.com/dustin/go-humanize"12 "github.com/dustin/seriesly/serieslyclient"13)14var (15 verbose = flag.Bool("v", false, "verbose")16 short = flag.Bool("short", false, "short form")17 templateSrc = flag.String("t", "", "Display template")18 templateFile = flag.String("T", "", "Display template filename")19)20const defaultTemplate = `{{range .}}21{{.DBName}}:22 Space Used: {{.SpaceUsed|bytes}}23 Last Sequence: {{.SpaceUsed|comma}}24 Header Position: {{.HeaderPos|comma}}25 Document Count: {{.DocCount|comma}}26 Deleted Count: {{.DeletedCount|comma}}27{{end}}28`29const shortTemplate = `dbname docs space used30------ ---- ----------31{{range .}}{{.DBName}} {{.DocCount|comma}} {{.SpaceUsed|bytes}}32{{end}}`33var funcMap = template.FuncMap{34 "comma": func(n json.Number) string {35 nint, err := n.Int64()36 maybeFatal(err, "Invalid int64: %v: %v", n, err)37 return humanize.Comma(nint)38 },39 "bytes": func(n json.Number) string {40 nint, err := n.Int64()41 maybeFatal(err, "Invalid int64: %v: %v", n, err)42 return humanize.Bytes(uint64(nint))43 }}44func init() {45 log.SetFlags(log.Lmicroseconds)46 flag.Usage = func() {47 fmt.Fprintf(os.Stderr, "Usage: %v [-v] http://serieslyhost:3133/ [dbnames...]\n",48 os.Args[0])49 flag.PrintDefaults()50 }51}52func maybeFatal(err error, fmt string, args ...interface{}) {53 if err != nil {54 log.Fatalf(fmt, args...)55 }56}57func vlog(s string, a ...interface{}) {58 if *verbose {59 log.Printf(s, a...)60 }61}62func describe(s *serieslyclient.Seriesly, dbs ...string) <-chan *serieslyclient.DBInfo {63 rv := make(chan *serieslyclient.DBInfo)64 go func() {65 defer close(rv)66 for _, db := range dbs {67 di, err := s.DB(db).Info()68 maybeFatal(err, "Couldn't fetch info for %v: %v", db, err)69 di.DBName = db70 rv <- di71 }72 }()73 return rv74}75func getTemplate(tdefault string) *template.Template {76 tmplstr := *templateSrc77 if tmplstr == "" {78 switch *templateFile {79 case "":80 tmplstr = tdefault81 case "-":82 td, err := ioutil.ReadAll(os.Stdin)83 maybeFatal(err, "Error reading template from stdin: %v", err)84 tmplstr = string(td)85 default:86 td, err := ioutil.ReadFile(*templateFile)87 maybeFatal(err, "Error reading template file: %v", err)88 tmplstr = string(td)89 }90 }91 tmpl, err := template.New("").Funcs(funcMap).Parse(tmplstr)92 maybeFatal(err, "Error parsing template: %v", err)93 return tmpl94}95func main() {96 flag.Parse()97 if flag.NArg() < 1 {98 flag.Usage()99 os.Exit(64)100 }101 s, err := serieslyclient.New(flag.Arg(0))102 maybeFatal(err, "Couldn't set up client: %v", err)103 dbs := flag.Args()[1:]104 if len(dbs) == 0 {105 dbs, err = s.List()106 maybeFatal(err, "Error listing DBs: %v", err)107 }108 tsrc := defaultTemplate109 if *short {110 tsrc = shortTemplate111 }112 tmpl := getTemplate(tsrc)113 tw := tabwriter.NewWriter(os.Stdout, 8, 8, 2, ' ', 0)114 tmpl.Execute(tw, describe(s, dbs...))115 defer tw.Flush()116}...

Full Screen

Full Screen

nInt

Using AI Code Generation

copy

Full Screen

1import (2func Walk(t *tree.Tree, ch chan int) {3 if t.Left != nil {4 Walk(t.Left, ch)5 }6 if t.Right != nil {7 Walk(t.Right, ch)8 }9}10func Same(t1, t2 *tree.Tree) bool {11 ch1 := make(chan int)12 ch2 := make(chan int)13 go Walk(t1, ch1)14 go Walk(t2, ch2)15 for i := 0; i < 10; i++ {16 if <-ch1 != <-ch2 {17 }18 }19}20func main() {21 ch := make(chan int)22 go Walk(tree.New(1), ch)23 for i := 0; i < 10; i++ {24 fmt.Println(<-ch)25 }26 fmt.Println(Same(tree.New(1), tree.New(1)))27 fmt.Println(Same(tree.New(1), tree.New(2)))28}

Full Screen

Full Screen

nInt

Using AI Code Generation

copy

Full Screen

1import (2type td struct {3}4func (t td) nInt() int {5}6func main() {7 t := td{a: 10}8 fmt.Println(t.nInt())9}

Full Screen

Full Screen

nInt

Using AI Code Generation

copy

Full Screen

1import "fmt"2func main() {3 fmt.Println(x.nInt())4}5import "fmt"6func main() {7 fmt.Println(x.nFloat())8}9import "fmt"10func main() {11 fmt.Println(x.nStr())12}13import "fmt"14func main() {15 fmt.Println(x.nBool())16}17import "fmt"18func main() {19 x = [3]int{1, 2, 3}20 fmt.Println(x.nArr())21}22import "fmt"23func main() {24 x = map[string]int{"a": 1, "b": 2, "c": 3}25 fmt.Println(x.nMap())26}27import "fmt"28func main() {29 fmt.Println(x.nPtr())30}31import "fmt"32func main() {33 x = make(chan int)34 fmt.Println(x.nChan())35}36import "fmt"37func main() {38 x = func() {}39 fmt.Println(x.nFunc())40}41import "fmt"42func main() {43 x = struct{}{}44 fmt.Println(x.nStruct())

Full Screen

Full Screen

nInt

Using AI Code Generation

copy

Full Screen

1import "fmt"2func main() {3 fmt.Println("Hello, playground")4 fmt.Println(td1.nInt())5}6import "fmt"7type td struct {8}9func main() {10 fmt.Println("Hello, playground")11}

Full Screen

Full Screen

nInt

Using AI Code Generation

copy

Full Screen

1import (2func main() {3 t = td.New(1, 2, 3, 4, 5)4 fmt.Println(t)5 fmt.Println(t.NInt(0))6 fmt.Println(t.NInt(1))7 fmt.Println(t.NInt(2))8 fmt.Println(t.NInt(3))9 fmt.Println(t.NInt(4))10}11import (12func main() {13 t := td.New(1, 2, 3, 4, 5)14 fmt.Println(t)15 fmt.Println(t.NInt(0))16 fmt.Println(t.NInt(1))17 fmt.Println(t.NInt(2))18 fmt.Println(t.NInt(3))19 fmt.Println(t.NInt(4))20}21func (t *TD) NInt(n int) int {22}23import (24func main() {25 t := td.New(1, 2, 3, 4, 5)26 fmt.Println(t)27 fmt.Println(t.NInt(0))28 fmt.Println(t.NInt(1))29 fmt.Println(t.NInt(2))

Full Screen

Full Screen

nInt

Using AI Code Generation

copy

Full Screen

1import (2func main() {3 td.set(3, 4, 5)4 fmt.Println(td.nInt())5}6import (7func main() {8 td.set(3, 4, 5)9 fmt.Println(td.nInt())10}11import (12func main() {13 td.set(3, 4, 5)14 fmt.Println(td.nInt())15}16import (17func main() {18 td.set(3, 4, 5)19 fmt.Println(td.nInt())20}21import (22func main() {23 td.set(3, 4, 5)24 fmt.Println(td.nInt())25}26import (27func main() {28 td.set(3, 4, 5)29 fmt.Println(td.nInt())30}31import (32func main() {33 td.set(3, 4, 5)34 fmt.Println(td.nInt())35}36import (37func main() {38 td.set(3, 4, 5)39 fmt.Println(td.nInt())40}41import (

Full Screen

Full Screen

nInt

Using AI Code Generation

copy

Full Screen

1import (2func main() {3 t.Init(10, 10)4 t.Set(1, 1, 1)5 t.Set(1, 2, 2)6 t.Set(1, 3, 3)7 t.Set(1, 4, 4)8 t.Set(1, 5, 5)9 t.Set(1, 6, 6)10 t.Set(1, 7, 7)11 t.Set(1, 8, 8)12 t.Set(1, 9, 9)13 t.Set(1, 10, 10)14 t.Set(2, 1, 11)15 t.Set(2, 2, 12)16 t.Set(2, 3, 13)17 t.Set(2, 4, 14)18 t.Set(2, 5, 15)19 t.Set(2, 6, 16)20 t.Set(2, 7, 17)21 t.Set(2, 8, 18)22 t.Set(2, 9, 19)23 t.Set(2, 10, 20)24 t.Set(3, 1, 21)25 t.Set(3, 2, 22)26 t.Set(3, 3, 23)27 t.Set(3, 4, 24)28 t.Set(3, 5, 25)29 t.Set(3, 6, 26)30 t.Set(3, 7, 27)31 t.Set(3, 8, 28)32 t.Set(3, 9, 29)33 t.Set(3, 10, 30)34 t.Set(4, 1, 31)35 t.Set(4, 2, 32)36 t.Set(4, 3, 33)37 t.Set(4, 4, 34)38 t.Set(4, 5, 35)39 t.Set(4, 6, 36)40 t.Set(4, 7, 37)41 t.Set(4, 8, 38)42 t.Set(4, 9

Full Screen

Full Screen

Automation Testing Tutorials

Learn to execute automation testing from scratch with LambdaTest Learning Hub. Right from setting up the prerequisites to run your first automation test, to following best practices and diving deeper into advanced test scenarios. LambdaTest Learning Hubs compile a list of step-by-step guides to help you be proficient with different test automation frameworks i.e. Selenium, Cypress, TestNG etc.

LambdaTest Learning Hubs:

YouTube

You could also refer to video tutorials over LambdaTest YouTube channel to get step by step demonstration from industry experts.

Run Go-testdeep automation tests on LambdaTest cloud grid

Perform automation testing on 3000+ real desktop and mobile devices online.

... View sample repo... View DocumentationSign up Free
›
_

Most used method in

Try LambdaTest Now !!

Get 100 minutes of automation test minutes FREE!!

...
Sign up with GoogleSign up with Email
Next-Gen App & Browser Testing Cloud

Was this article helpful?

Helpful

NotHelpful