-
Notifications
You must be signed in to change notification settings - Fork 0
/
gradientCalculator.go
160 lines (137 loc) · 4.24 KB
/
gradientCalculator.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
package pf
import (
"fmt"
"math"
"github.com/davidkleiven/gopf/pfutil"
)
// GradientCalculator calculates the gradient of a field
type GradientCalculator struct {
FT FourierTransform
Comp int
KeepNyquist bool
}
// Calculate calculates the gradient of the data passed
// data contain the field in real-space
func (g *GradientCalculator) Calculate(indata []complex128, data []complex128) {
copy(data, indata)
g.FT.FFT(data)
for i := range data {
f := g.FT.Freq(i)[g.Comp]
if math.Abs(f-0.5) < 1e-10 && !g.KeepNyquist {
f = 0.0
}
data[i] *= complex(0.0, 2.0*math.Pi*f)
}
g.FT.IFFT(data)
pfutil.DivRealScalar(data, float64(len(data)))
}
// ToDerivedField constructs a derived field from the gradient calculator.
// N is the number of grid points and brick is the brick that should be
// differentiated
func (g *GradientCalculator) ToDerivedField(name string, N int, brick Brick) DerivedField {
return DerivedField{
Data: make([]complex128, N),
Name: name,
Calc: func(data []complex128) {
for i := range data {
data[i] = brick.Get(i)
}
g.Calculate(data, data)
},
}
}
// DivGrad is a type used to represent the term Div F(c)Grad <field>
type DivGrad struct {
Field string
F GenericFunction
}
// FuncName returns the name of the generic function F
func (dg *DivGrad) FuncName() string {
return fmt.Sprintf("DivGrad_%s_Func", dg.Field)
}
// GradName returns the name of the gradient
func (dg *DivGrad) GradName(comp int) string {
return fmt.Sprintf("GRAD_%s_%d", dg.Field, comp)
}
// PrepareModel impose nessecary changes in the model in order to use the
// DivGrad term. N is the number of grid points in the simulation domain,
// dim is the dimension of the simulation domain. FT is a fourier transformer
// required for gradient evaluations
func (dg *DivGrad) PrepareModel(N int, m *Model, FT FourierTransform) {
dim := len(FT.Freq(0))
for d := 0; d < dim; d++ {
grad := GradientCalculator{
FT: FT,
Comp: d,
KeepNyquist: false,
}
dField := grad.ToDerivedField(dg.GradName(d), N, m.Bricks[dg.Field])
m.RegisterDerivedField(dField)
dField2 := DerivedField{
Name: dg.FuncName() + dg.GradName(d),
Data: make([]complex128, N),
Calc: func(data []complex128) {
for i := range data {
data[i] = dg.F(i, m.Bricks) * m.Bricks[dg.GradName(grad.Comp)].Get(i)
}
},
}
m.RegisterDerivedField(dField2)
}
}
// Construct builds the right hand side term
func (dg *DivGrad) Construct(bricks map[string]Brick) Term {
return func(freq Frequency, t float64, field []complex128) {
pfutil.Clear(field)
dim := len(freq(0))
for d := 0; d < dim; d++ {
brick := bricks[dg.FuncName()+dg.GradName(d)]
for i := range field {
f := freq(i)[d]
field[i] += complex(0.0, 2.0*math.Pi*f) * brick.Get(i)
}
}
}
}
// OnStepFinished does nothing as we don't need any updates in between steps
func (dg *DivGrad) OnStepFinished(t float64) {}
// WeightedLaplacian is a type used to represent terms of the form
// F(c) LAP <field>, where F is a function of all the fields.
type WeightedLaplacian struct {
// Name of the field to the right of the laplacian (must be registered in the model)
Field string
// Name of the function in front of the laplacian (must be registered in the model)
PreFactor string
// FourierTransformer of the correct domain size. This is used to
// go back and fourth between the real domain and the fourier domain
// internally
FT FourierTransform
}
// Construct returns the fourier transformed representation of this term
func (wl *WeightedLaplacian) Construct(bricks map[string]Brick) Term {
return func(freq Frequency, t float64, field []complex128) {
fieldBrick := bricks[wl.Field]
// Calcualte real space laplacian of the field
lap := LaplacianN{
Power: 1,
}
for i := range field {
field[i] = fieldBrick.Get(i)
}
lap.Eval(freq, field)
wl.FT.IFFT(field)
pfutil.DivRealScalar(field, float64(len(field)))
// Calculate the realspace representation of the prefactor
work := make([]complex128, len(field))
funcBrick := bricks[wl.PreFactor]
for i := range work {
work[i] = funcBrick.Get(i)
}
wl.FT.IFFT(work)
pfutil.DivRealScalar(work, float64(len(work)))
for i := range field {
field[i] *= work[i]
}
wl.FT.FFT(field)
}
}