Mercurial > repos > public > sbplib_julia

comparison Notes.md @ 1858:4a9be96f2569 feature/documenter_logo

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 12 Jan 2025 21:18:44 +0100
parents fe058a0ebd97
children
comparison
equal deleted inserted replaced
1857:ffde7dad9da5 1858:4a9be96f2569
1 # Notes 1 # Notes
2
3 ## How to dispatch for different operators
4 We have a problem in how dispatch for different operators work.
5 * We want to keep the types simple and flat (Awkward to forward `apply`)
6 * We want to dispatch SATs on the parameters of the continuous operator. (a * div for example)
7 * We want to allow keeping the same stencil_set across different calls. (maybe not so bad for the user to be responsible)
8
9 Could remove the current opset idea and introduce a description of continuous operators
10 ```julia
11 abstract type DifferentialOperator end
12
13 struct Laplace <: DifferentialOperator end
14 struct Advection <: DifferentialOperator
15 v
16 end
17
18 difference_operator(::Laplace, grid, stencil_set) = ... # Returns a plain LazyTensor. Replaces the current `laplace()` function.
19 sat_tensors(::Laplace, grid, stencil_set, bc) = ...
20
21 sat(::DifferentialOperator, grid, stencil_set, bc) = ...
22 ```
23
24
25 ### Update 2024-06-26
26 We will run into trouble if we start assuming things about the coupling
27 between the continuous and discrete setting. We could add representations of
28 continuous operators but we will also need representations of discrete
29 operators. Ideally it should be possible to ignore the continuous
30 representations and only work with the discrete operators without losing
31 functionality. The discrete representations does not have to be LazyTensors.
32 The could be used as inputs to methods for `sat`, `difference_operator` and so
33 on.
34
35 To see need for a fully functional discrete layer we can consider the
36 optimization of material parameters or something similar. In this case we do
37 not necessarily want to handle continuous objects.
2 38
3 ## Reading operators 39 ## Reading operators
4 40
5 Jonatan's suggestion is to add methods to `Laplace`, `SecondDerivative` and 41 Jonatan's suggestion is to add methods to `Laplace`, `SecondDerivative` and
6 similar functions that take in a filename from which to read stencils. These 42 similar functions that take in a filename from which to read stencils. These
69 If possible the goal should be for the parsing to get all the way to the 105 If possible the goal should be for the parsing to get all the way to the
70 stencils so that a user calls `read_stencil_set` and gets a 106 stencils so that a user calls `read_stencil_set` and gets a
71 dictionary-structure containing stencils, tuples, scalars and other types 107 dictionary-structure containing stencils, tuples, scalars and other types
72 ready for input to the methods creating the operators. 108 ready for input to the methods creating the operators.
73 109
74 ## Variable second derivative
75
76 2020-12-08 after discussion with Vidar:
77 We will have to handle the variable second derivative in a new variant of
78 VolumeOperator, "SecondDerivativeVariable?". Somehow it needs to know about
79 the coefficients. They should be provided as an AbstractVector. Where they are
80 provided is another question. It could be that you provide a reference to the
81 array to the constructor of SecondDerivativeVariable. If that array is mutable
82 you are free to change it whenever and the changes should propagate
83 accordingly. Another option is that the counter part to "Laplace" for this
84 variable second derivate returns a function or acts like a functions that
85 takes an Abstract array and returns a SecondDerivativeVariable with the
86 appropriate array. This would allow syntax like `D2(a)*v`. Can this be made
87 performant?
88
89 For the 1d case we can have a constructor
90 `SecondDerivativeVariable(D2::SecondDerivativeVariable, a)` that just creates
91 a copy with a different `a`.
92
93 Apart from just the second derivative in 1D we need operators for higher
94 dimensions. What happens if a=a(x,y)? Maybe this can be solved orthogonally to
95 the `D2(a)*v` issue, meaning that if a constant nD version of
96 SecondDerivativeVariable is available then maybe it can be wrapped to support
97 function like syntax. We might have to implement `SecondDerivativeVariable`
98 for N dimensions which takes a N dimensional a. If this could be easily
99 closured to allow D(a) syntax we would have come a long way.
100
101 For `Laplace` which might use a variable D2 if it is on a curvilinear grid we
102 might want to choose how to calculate the metric coefficients. They could be
103 known on closed form, they could be calculated from the grid coordinates or
104 they could be provided as a vector. Which way you want to do it might change
105 depending on for example if you are memory bound or compute bound. This choice
106 cannot be done on the grid since the grid shouldn't care about the computer
107 architecture. The most sensible option seems to be to have an argument to the
108 `Laplace` function which controls how the coefficients are gotten from the
109 grid. The argument could for example be a function which is to be applied to
110 the grid.
111
112 What happens if the grid or the varible coefficient is dependent on time?
113 Maybe it becomes important to support `D(a)` or even `D(t,a)` syntax in a more
114 general way.
115
116 ```
117 g = TimeDependentGrid()
118 L = Laplace(g)
119 function Laplace(g::TimeDependentGrid)
120 g_logical = logical(g) # g_logical is time independent
121 ... Build a L(a) assuming we can do that ...
122 a(t) = metric_coeffs(g,t)
123 return t->L(a(t))
124 end
125 ```
126
127 ## Known size of range and domain? 110 ## Known size of range and domain?
128 Is there any reason to use a trait to differentiate between fixed size and unknown size? 111 Is there any reason to use a trait to differentiate between fixed size and unknown size?
129 112
130 When do we need to know the size of the range and domain? 113 When do we need to know the size of the range and domain?
131 * When indexing to provide boundschecking? 114 * When indexing to provide boundschecking?
132 * When doing specialised computations for different parts of the range/domain? 115 * When doing specialised computations for different parts of the range/domain?
133 * More? 116 * More?
134 117
135 Maybe if we should have dynamic sizing it could be only for the range. `domain_size` would not be implemented. And the `range_size` would be a function of a vector that the TensorMapping is applied to. 118 Maybe if we should have dynamic sizing it could be only for the range. `domain_size` would not be implemented. And the `range_size` would be a function of a vector that the LazyTensor is applied to.
136 119
137 ## Reasearch and thinking 120 ## Reasearch and thinking
138 - [ ] Use a trait to indicate that a TensorMapping har the same range and domain?
139 - [ ] Rename all the Tensor stuff to just LazyOperator, LazyApplication and so on?
140 - [ ] Check how the native julia doc generator works 121 - [ ] Check how the native julia doc generator works
141 - [ ] Check if Vidars design docs fit in there 122 - [ ] Check if Vidars design docs fit in there
142 - [ ] Create a macro @lazy which replaces a binary op (+,-) by its lazy equivalent? Would be a neat way to indicate which evaluations are lazy without cluttering/confusing with special characters. 123 - [ ] Create a macro @lazy which replaces a binary op (+,-) by its lazy equivalent? Would be a neat way to indicate which evaluations are lazy without cluttering/confusing with special characters.
143 - [ ] Specificera operatorer i TOML eller något liknande? 124 - [ ] Can we have a trait to tell if a LazyTensor is transposable?
144 H.. H_gamma etc.)
145 - [ ] Dispatch on Lower() instead of the type Lower so `::Lower` instead of `::Type{Lower}` ???
146 Seems better unless there is some specific reason to use the type instead of the value.
147 - [ ] How do we handle mixes of periodic and non-periodic grids? Seems it should be supported on the grid level and on the 1d operator level. Between there it should be transparent.
148 - [ ] Can we have a trait to tell if a TensorMapping is transposable?
149 - [ ] Is it ok to have "Constructors" for abstract types which create subtypes? For example a Grids() functions that gives different kind of grids based on input?
150 125
151 ## Regions and tensormappings 126 ## Regions and tensormappings
152 - [ ] Use a trait to indicate if a TensorMapping uses indices with regions. 127 - [ ] Use a trait to indicate if a LazyTensor uses indices with regions.
153 The default should be that they do NOT. 128 The default should be that they do NOT.
154 - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle? 129 - [ ] What to name this trait? Can we call it IndexStyle but not export it to avoid conflicts with Base.IndexStyle?
155 - [ ] Figure out repeated application of regioned TensorMappings. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work? 130 - [ ] Figure out repeated application of regioned LazyTensors. Maybe an instance of a tensor mapping needs to know the exact size of the range and domain for this to work?
131
132 ### Ideas for information sharing functions
133 ```julia
134 using StaticArrays
135
136 function regions(op::SecondDerivativeVariable)
137 t = ntuple(i->(Interior(),),range_dim(op))
138 return Base.setindex(t, (Lower(), Interior(), Upper()), derivative_direction(op))
139 end
140
141 function regionsizes(op::SecondDerivativeVariable)
142 sz = tuple.(range_size(op))
143
144 cl = closuresize(op)
145 return Base.setindex(sz, (cl, n-2cl, cl), derivative_direction(op))
146 end
147
148
149 g = EquidistantGrid((11,9), (0.,0.), (10.,8.)) # h = 1
150 c = evalOn(g, (x,y)->x+y)
151
152 D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,1)
153 @test regions(D₂ᶜ) == (
154 (Lower(), Interior(), Upper()),
155 (Interior(),),
156 )
157 @test regionsizes(D₂ᶜ) == ((1,9,1),(9,))
158
159
160 D₂ᶜ = SecondDerivativeVariable(g, c, interior_stencil, closure_stencils,2)
161 @test regions(D₂ᶜ) == (
162 (Interior(),),
163 (Lower(), Interior(), Upper()),
164 )
165 @test regionsizes(D₂ᶜ) == ((11,),(1,7,1))
166 ```
167
156 168
157 ## Boundschecking and dimension checking 169 ## Boundschecking and dimension checking
158 Does it make sense to have boundschecking only in getindex methods? 170 Does it make sense to have boundschecking only in getindex methods?
159 This would mean no bounds checking in applys, however any indexing that they do would be boundschecked. The only loss would be readability of errors. But users aren't really supposed to call apply directly anyway. 171 This would mean no bounds checking in applys, however any indexing that they do would be boundschecked. The only loss would be readability of errors. But users aren't really supposed to call apply directly anyway.
160 172
161 Preferably dimensions and sizes should be checked when lazy objects are created, for example TensorApplication, TensorComposition and so on. If dimension checks decreases performance we can make them skippable later. 173 Preferably dimensions and sizes should be checked when lazy objects are created, for example TensorApplication, TensorComposition and so on. If dimension checks decreases performance we can make them skippable later.
162 174
175 ## Changes to `eval_on`
176 There are reasons to replace `eval_on` with regular `map` from Base, and
177 implement a kind of lazy map perhaps `lmap` that work on indexable
178 collections.
179
180 The benefit of doing this is that we can treat grids as gridfunctions for the
181 coordinate function, and get a more flexible tool. For example `map`/`lmap`
182 can then be used both to evaluate a function on the grid but also get a
183 component of a vector valued grid function or similar.
184
185 Below is a partial implementation of `lmap` with some ideas
186 ```julia
187 struct LazyMapping{T,IT,F}
188 f::F
189 indexable_iterator::IT # ___
190 end
191
192 function LazyMapping(f,I)
193 IT = eltype(I)
194 T = f(zero(T))
195 F = typeof(f)
196
197 return LazyMapping{T,IT,F}(f,I)
198 end
199
200 getindex(lm::LazyMapping, I...) = lm.f(lm.I[I...])
201 # indexabl interface
202 # iterable has shape
203
204 iterate(lm::LazyMapping) = _lazy_mapping_iterate(lm, iterate(lm.I))
205 iterate(lm::LazyMapping, state) = _lazy_mapping_iterate(lm, iterate(lm.I, state))
206
207 _lazy_mapping_iterate(lm, ::Nothing) = nothing
208 _lazy_mapping_iterate(lm, (next, state)) = lm.f(next), state
209
210 lmap(f, I) = LazyIndexableMap(f,I)
211 ```
212
213 The interaction of the map methods with the probable design of multiblock
214 functions involving nested indecies complicate the picture slightly. It's
215 unclear at the time of writing how this would work with `Base.map`. Perhaps we
216 want to implement our own versions of both eager and lazy map.
217
218
219 ### 2024-04
220 MappedArrays.jl provides a simple array type and function like the description
221 of LazyMapping above. One option is to remove `eval_on` completely and rely on
222 destructuring arguments if handling the function input as a vector is
223 undesirable.
224
225 If we can let multi-block grids be iterators over grid points we could even
226 handle those by specialized implementation of `map` and `mappedarray`.
227
228 ## Multiblock implementation
229 We want multiblock things to work very similarly to regular one block things.
230
231 ### Grid functions
232 Should probably support a nested indexing so that we first have an index for
233 subgrid and then an index for nodes on that grid. E.g `g[1,2][2,3]` or
234 `g[3][43,21]`.
235
236 We could also possibly provide a combined indexing style `g[1,2,3,4]` where
237 the first group of indices are for the subgrid and the remaining are for the
238 nodes.
239
240 We should make sure the underlying buffer for grid functions are continuously
241 stored and are easy to convert to, so that interaction with for example
242 DifferentialEquations is simple and without much boilerplate.
243
244 #### `map` and `collect` and nested indexing
245 We need to make sure `collect`, `map` and a potential lazy map work correctly
246 through the nested indexing. Also see notes on `eval_on` above.
247
248 Possibly this can be achieved by providing special nested indexing but not
249 adhering to an array interface at the top level, instead being implemented as
250 an iterator over the grid points. A custom trait can let map and other methods
251 know the shape (or structure) of the nesting so that they can efficiently
252 allocate result arrays.
253
254 ### Tensor applications
255 Should behave as grid functions
256
257 ### LazyTensors
258 Could be built as a tuple or array of LazyTensors for each grid with a simple apply function.
259
260 Nested indexing for these is problably not needed unless it simplifies their own implementation.
261
262 Possibly useful to provide a simple type that doesn't know about connections between the grids. Antother type can include knowledge of the.
263
264 We have at least two option for how to implement them:
265 * Matrix of LazyTensors
266 * Looking at the grid and determining what the apply should do.
267
268 ### Overall design implications of nested indices
269 If some grids accept nested indexing there might be a clash with how LazyArrays work. It would be nice if the grid functions and lazy arrays that actually are arrays can be AbstractArray and things can be relaxed for nested index types.
270
163 ## Vector valued grid functions 271 ## Vector valued grid functions
164 Från slack konversation:
165
166 Jonatan Werpers:
167 Med vektorvärda gridfunktioner vill vi ju fortfarande att grid funktionen ska vara till exempel AbstractArray{LitenVektor,2}
168 Och att man ska kunna göra allt man vill med LitenVektor
169 typ addera, jämföra osv
170 Och då borde points returnera AbstractArray{LitenVektor{Float,2},2} för ett 2d nät
171 Men det kanske bara ska vara Static arrays?
172
173 Vidar Stiernström:
174 Ja, jag vet inte riktigt vad som är en rimlig representation
175 Du menar en vektor av static arrays då?
176
177 Jonatan Werpers:
178 Ja, att LitenVektor är en StaticArray
179
180 Vidar Stiernström:
181 Tuplar känns typ rätt inuitivt för att representera värdet i en punkt
182 men
183 det suger att man inte har + och - för dem
184
185 Jonatan Werpers:
186 Ja precis
187
188 Vidar Stiernström:
189 så kanske är bra med static arrays i detta fall
190
191 Jonatan Werpers:
192 Man vill ju kunna köra en Operator rakt på och vara klar eller?
193
194 Vidar Stiernström:
195 Har inte alls tänkt på hur det vi gör funkar mot vektorvärda funktioner
196 men känns som staticarrays är hur man vill göra det
197 tuplar är ju immutable också
198 blir jobbigt om man bara agerar på en komponent då
199
200 Jonatan Werpers:
201 Hm…
202 Tål att tänkas på
203 Men det lär ju bli mer indirektion med mutables eller?
204 Hur fungerar det?
205 Det finns ju hur som helst både SVector och MVector i StaticArrays
206
207 Vidar Stiernström:
208 När vi jobbat i c/c++ och kollat runt lite hur man brukar göra så lagrar man i princip alla sina obekanta i en lång vektor och så får man specificera i funktioerna vilken komponent man agerar på och till vilken man skriver
209 så man lagrar grejer enl: w = [u1, v1, u2, v2, …] i 1D.
210 Men alltså har ingen aning hur julia hanterar detta
211
212 Jonatan Werpers:
213 Det vi är ute efter kanske är att en grid funcktion är en AbstractArray{T,2} where T<:NågotSomViKanRäknaMed
214 Och så får den typen var lite vad som helst.
215
216 Vidar Stiernström:
217 Tror det kan vara farligt att ha nåt som är AbstractArray{LitenArray{NDof},Dim}
218 Jag gissar att det kompilatorn vill ha är en stor array med doubles
219
220 Jonatan Werpers:
221 Och sen är det upp till den som använder grejerna att vara smart
222 Vill man vara trixig kan man väl då imlementera SuperHaxxorGridFunction <: AbstractArray{Array{…},2} som lagrar allt linjärt eller något sånt
223 Det kommer väl lösa sig när man börjar implementera vektorvärda saker
224 Euler nästa!
225 New
226 Vidar Stiernström:
227 Det vore skönt att inte behöva skriva såhär varje gång man testar mot en tupel :smile: @test [gp[i]...] ≈ [p[i]...] atol=5e-13
228
229 Jonatan Werpers:
230 https://github.com/JuliaArrays/ArraysOfArrays.jl
231 https://github.com/jw3126/Setfield.jl
232 272
233 ### Test-applikationer 273 ### Test-applikationer
234 div och grad operationer 274 div- och grad-operationer
235 275
236 Enligt Wikipedia verkar det som att `∇⋅` agerar på första dimensionen av ett tensor fält och `div()` på sista. 276 Enligt Wikipedia verkar det som att `∇⋅` agerar på första dimensionen av ett tensorfält och `div()` på sista.
237 Om man generaliserar kanske `∇` i så fall bara lägger till en dimension i början. 277 Om man generaliserar kanske `∇` i så fall bara lägger till en dimension i början.
238 278
239 Kan vi implementera `⋅`(\cdot) så att de fungerar som man vill för både tensor-fält och tensor-operatorer? 279 Kan vi implementera `⋅`(\cdot) så att de fungerar som man vill för både tensor-fält och tensor-operatorer?
240 280
241 Är `∇` ett tensor-fält av tensor-operatorer? Vad är ett tensor-fält i vår kod? Är det en special-fall av en tensor-mapping? 281 Är `∇` ett tensor-fält av tensor-operatorer? Vad är ett tensor-fält i vår kod? Är det en special-fall av en tensor-mapping?
242 282
243 ### Grid-funktionen 283 ### Grid-funktionen
244 Grid-funktionon har typen `AbstractArray{T,2} where T`. 284 Grid-funktioner har typen `AbstractArray{T,2} where T`.
245 `T` kan vara lite vad som helst, tillexemel en SVector eller Array, eller tuple. TensorOperatorerna bryr sig inte om exakt vad det är, mer än att typen måste stödja de operationer som operatorn använder. 285 `T` kan vara lite vad som helst, tillexemel en SVector eller Array, eller Tuple. Tensoroperatorerna bryr sig inte om exakt vad det är, mer än att typen måste stödja de operationer som operatorn använder.
246
247 En nackdel kan vara hur man ska få ut gridfunktionen för tex andra komponenten.
248
249 Syntax:
250 ```
251 f(x̄) = x̄
252 gf = evalOn(g, f)
253 gf[2,3] # x̄ för en viss gridpunkt
254 gf[2,3][2] # x̄[2] för en viss gridpunkt
255 ```
256
257 Note: Behöver bestämma om eval on skickar in `x̄` eller `x̄...` till `f`. Eller om man kan stödja båda.
258 286
259 ### Tensor operatorer 287 ### Tensor operatorer
260 Vi kan ha tensor-operatorer som agerar på ett skalärt fält och ger ett vektorfält eller tensorfält. 288 Vi kan ha tensor-operatorer som agerar på ett skalärt fält och ger ett vektorfält eller tensorfält.
261 Vi kan också ha tensor-operatorer som agerar på ett vektorfält eller tensorfält och ger ett skalärt fält. 289 Vi kan också ha tensor-operatorer som agerar på ett vektorfält eller tensorfält och ger ett skalärt fält.
262 290
263 TBD: Just nu gör `apply_transpose` antagandet att domän-typen är samma som range-typen. Det behöver vi på något sätt bryta. Ett alternativ är låta en TensorMapping ha `T_domain` och `T_range` istället för bara `T`. Känns dock lite grötigt. Ett annat alternativ skulle vara någon typ av trait för transpose? Den skulle kunna innehålla typen som transponatet agerar på? Vet inte om det fungerar dock. 291 TBD: Just nu gör `apply_transpose` antagandet att domän-typen är samma som range-typen. Det behöver vi på något sätt bryta. Ett alternativ är låta en LazyTensor ha `T_domain` och `T_range` istället för bara `T`. Känns dock lite grötigt. Ett annat alternativ skulle vara någon typ av trait för transpose? Den skulle kunna innehålla typen som transponatet agerar på? Vet inte om det fungerar dock.
264 292
265 TBD: Vad är målet med `T`-parametern för en TensorMapping? Om vi vill kunna applicera en difference operator på vad som helst kan man inte anta att en `TensorMapping{T}` bara agerar på instanser av `T`. 293 TBD: Vad är målet med `T`-parametern för en LazyTensor? Om vi vill kunna applicera en difference operator på vad som helst kan man inte anta att en `LazyTensor{T}` bara agerar på instanser av `T`.
266 294
267 Man kan implementera `∇` som en tensormapping som agerar på T och returnerar `StaticVector{N,T} where N`. 295 Man kan implementera `∇` som en tensormapping som agerar på T och returnerar `StaticVector{N,T} where N`.
268 (Man skulle eventuellt också kunna låta den agera på `StaticMatrix{N,T,D} where N` och returnera `StaticMatrix{M,T,D+1}`. Frågan är om man vinner något på det...) 296 (Man skulle eventuellt också kunna låta den agera på `StaticMatrix{N,T,D} where N` och returnera `StaticMatrix{M,T,D+1}`. Frågan är om man vinner något på det...)
269 297
270 Skulle kunna ha en funktion `range_type(::TensorMapping, ::Type{domain_type})` 298 Skulle kunna ha en funktion `range_type(::LazyTensor, ::Type{domain_type})`
271 299
272 Kanske kan man implementera `⋅(tm::TensorMapping{R,D}, v::AbstractArray{T,D})` där T är en AbstractArray, tm på något sätt har komponenter, lika många som T har element. 300 Kanske kan man implementera `⋅(tm::LazyTensor{R,D}, v::AbstractArray{T,D})` där T är en AbstractArray, tm på något sätt har komponenter, lika många som T har element.
273 301
274 ### Ratade alternativ: 302 ### Prestanda-aspekter
275 303 [Vidar, Discord, 2023-03-03]
276 304 Typiskt sett finns det två sätt att representera vektorvärda gridfunktioner AbstractArray{T,Dim} där T är en vektor över komponenterna. Man skulle alltså i 1D ha
277 #### 2.AbstractArray{T,2+1} where T (NOPE!) 305 u = [ [u1[x1], u2[x1]] , [u1[x2], u2[x2]], ... [u1[xN], u2[xN]]]. Detta brukar kallas array of structs (AoS). Alternativet är struct of arrays (SoA), där man har alla gridpunkter för en given komponent u = [[u1[x1], u1[x2]],... u1[xN]], [u2[x1], u2[x2], ... u2[xN]]].
278 Blir inte den här. Bryter mot alla tankar om hur grid funktioner ska fungera. Om de tillåts ha en annan dimension än nätet blir allt hemskt. 306
279 307 Personligen tycker jag att AoS känns som den mer naturliga representationen? Det skulle göra det enklarare att parallelisera en vektorvärd gridfunktion över gridpunkterna, och om man opererar på olika komponenter i samma funktion så är det också bra ur en minnesaccess-synpunkt då dessa kommer ligga nära vandra i minnet. Problemet är att AoS sabbar vektorisering på CPU då två gridpunkter i en komponent ligger långt bort från varandra. Efter lite eftersökningar (och efter att snackat lite med Ossian) så verkar det ändå som att AoS är dåligt på GPU, där man vill att trådar typiskt sett utföra samma operation på närliggande minne.
280 Man låter helt enkelt arrayen ha en extra dimension. En fördel är att man har en väldigt "native" typ. En nackdel kan vara att det eventuellt blir rörigt vilken dimension olika operatorer ska agera på. I värsta fall behöver vi "kroneckra in" de tillagda dimensionerna. Vektorfältets index kommer också att bli det första eftersom vi vill att de ska lagras kontinuerligt i minnet pga chachen. (Går kanske att lösa med en custom typ men då krånglar man till det för sig). En fördel skulle vara att man enkelt får ut olika komponenter. 308
281 309 Vad tänker du kring detta ur ett interface-perspektiv? Jag hittade paketet https://github.com/JuliaArrays/StructArrays.jl som verkar erbjuda AoS-interface men SoA-minneslayout så det kanske kan vara något vi kan använda? Inte native-stödd på samma sätt som SVector, men verkar iaf utvecklas aktivt.
282 Syntax: 310
283 ``` 311 [Efter telefonsamtal] För optimal prestanda behöver vi antagligen se till att man kan räkna ut varje komponent i en punkt individuellt. Detta så att man har frihet att till exempel låta den innersta loopen hålla komponentindexet konstant för att underlätta intruktionsvektorisering.
284 gf = eval_on_grid(g,f) 312
285 gf[:,2,3] # Hela vektorn för en gridpunkt 313
286 gf[2,2,3] # Andra komponenten av vektor fältet i en punkt. 314 [Vidare tankar]
287 gf[2,:,:] # 315 * Det borde bara vara output-gridfunktionen som behöver special-indexeras? Det viktiga på inputsidan är att den är lagrad på rätt sätt i minnet.
288 ``` 316 * Det borde inte vara några problem att behålla det "optimala" interfacet (gf[1,1,1][2]) till gridfunktionerna. Om man verkligen behöver kan skapa parallella indexeringsmetoder som gör det man behöver, i.e, "deep indexing".
289 317 * Det är inte säkert att vi behöver göra något speciellt på outputsidan överhuvudtaget. Det känns inte orimligt att kompilatorn skulle kunna optimera bort den koden som räknar ut onödiga komponenter.
290 ### Evaluering av funktioner på nät 318 * Om vi behöver special-indexering kommer till exempel LazyTensorApplication att behöva implementera det.
291 Hur ska man skriva funktioner som evalueras på nätet? `f(x,y) = ...` eller `f(x̄) = ...`? Eller båda? Kan eval_on_grid se skillnad eller får användaren specificera? 319 * För att komma vidare med något mer avancerat behöver vi antagligen implementera några operatorer som ger och agerar på vektorvärda funktioner. Tex grad, elastiska operatorn, andra?
292 320
293 ```
294 f(x,y) = [x^2, y^2]
295 f(x̄) = [x̄[1]^2, x̄[2]^2]
296 ```
297
298 Påverkas detta av hur vi förväntar oss kunna skapa lata gridfunktioner?
299
300 ### Komponenter som gridfunktioner
301 En viktig operation för vektor fält är att kunna få ut komponenter som grid-funktioner. Detta behöver antagligen kunna ske lazy.
302 Det finns ett par olika lösningar:
303 * Implementera en egen typ av view som tar hand om detta. Eller Accessors.jl?
304 * Använda en TensorMapping
305 * Någon typ av lazy-broadcast
306 * En lazy array som applicerar en funktion för varje element.
307
308 Skulle vara en fördel om det är hyffsat generiskt så att en eventuell användare kan utöka det enkelt om de har någon egen exotisk typ. Eller ska man vila helt på
309
310 Syntax:
311 ```
312 gf = eval(...)
313 component(gf,2) # Andra komponenten av en vektor
314 component(gf,2,3) # (2,3) elementet av en matris
315 component(gf,:,2) # Andra kolumnen av en matris
316 @ourview gf[:,:][2]
317 ```
318
319 ## Grids embedded in higher dimensions
320
321 For grids generated by asking for boundary grids for a regular grid, it would
322 make sense if these grids knew they were embedded in a higher dimension. They
323 would return coordinates in the full room. This would make sense when
324 drawing points for example, or when evaluating functions on the boundary.
325
326 Implementation of this is an issue that requires some thought. Adding an extra
327 "Embedded" type for each grid would make it easy to understand each type but
328 contribute to "type bloat". On the other hand adapting existing types to
329 handle embeddedness would complicate the now very simple grid types. Are there
330 other ways of doing the implentation?
331 321
332 ## Performance measuring 322 ## Performance measuring
333 We should be measuring performance early. How does our effective cpu and memory bandwidth utilization compare to peak performance? 323 We should be measuring performance early. How does our effective cpu and memory bandwidth utilization compare to peak performance?
334 324
335 We should make these test simple to run for any solver. 325 We should make these test simple to run for any solver.
346 336
347 A different approach would be to include it as a trait for operators so that you can specify what the adjoint for that operator is. 337 A different approach would be to include it as a trait for operators so that you can specify what the adjoint for that operator is.
348 338
349 339
350 ## Name of the `VolumeOperator` type for constant stencils 340 ## Name of the `VolumeOperator` type for constant stencils
351 It seems that the name is too general. The name of the method `volume_operator` makes sense. It should return different types of `TensorMapping` specialized for the grid. A suggetion for a better name is `ConstantStencilVolumeOperator` 341 It seems that the name is too general. The name of the method `volume_operator` makes sense. It should return different types of `LazyTensor` specialized for the grid. A suggetion for a better name is `ConstantStencilVolumeOperator`
342
343
344 ## Implementation of LazyOuterProduct
345 Could the implementation of LazyOuterProduct be simplified by making it a
346 struct containing two or more LazyTensors? (using split_tuple in a similar way
347 as TensorGrid)
348
349 ## Implementation of boundary_indices for more complex grids
350 To represent boundaries of for example tet-elements we can use a type `IndexCollection` to index a grid function directly.
351
352 ```julia
353 I = IndexCollection(...)
354 v[I]
355 ```
356
357 * This would impact how tensor grid works.
358 * To make things homogenous maybe these index collections should be used for the more simple grids too.
359 * The function `to_indices` from Base could be useful to implement for `IndexCollection`
360
361
362 ## Stencil application pipeline
363 We should make sure that `@inbounds` and `Base.@propagate_inbounds` are
364 applied correctly throughout the stack. When testing the performance of
365 stencil application on the bugfix/sbp_operators/stencil_return_type branch
366 there seemed to be some strange results where such errors could be the
367 culprit.
368
369
370 ## Tiled loops and regions in apply
371 There should be easy ways to use functionalty splitting the application of a lazy array into regions and using tiled iteration. This could make the application more efficient by reducing branching and improving cache usage in the tight loop. On commit f215ac2a5c66 and before there were some early tests regarding this in a DiffOp submodule.
372
373 The main ideas were:
374 ```julia
375 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
376 apply_region!(D, u, v, Lower, Lower)
377 apply_region!(D, u, v, Lower, Interior)
378 apply_region!(D, u, v, Lower, Upper)
379 apply_region!(D, u, v, Interior, Lower)
380 apply_region!(D, u, v, Interior, Interior)
381 apply_region!(D, u, v, Interior, Upper)
382 apply_region!(D, u, v, Upper, Lower)
383 apply_region!(D, u, v, Upper, Interior)
384 apply_region!(D, u, v, Upper, Upper)
385 return nothing
386 end
387 ```
388
389 ```julia
390 using TiledIteration
391 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
392 ri = regionindices(D.grid.size, closuresize(D.op), (r1,r2))
393 # TODO: Pass Tilesize to function
394 for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
395 for j ∈ tileaxs[2], i ∈ tileaxs[1]
396 I = ri[i,j]
397 u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
398 end
399 end
400 return nothing
401 end
402 ```