Skip to content

Instantly share code, notes, and snippets.

@leonid-shevtsov
Last active September 18, 2024 21:07
Show Gist options
  • Save leonid-shevtsov/5c99582c22f406a2f62b6290af5a366e to your computer and use it in GitHub Desktop.
Save leonid-shevtsov/5c99582c22f406a2f62b6290af5a366e to your computer and use it in GitHub Desktop.
// adapted from https://gitlab.com/manzhara157/nanobelts/-/blob/main/example.pas
package main
import (
"bufio"
"fmt"
"math"
"os"
"strconv"
"strings"
"github.com/samber/lo"
)
const (
n0 = 200000
tmax = 1000000.0
dt = 0.1
ra_kinet = 1.0
rb_kinet = 1.0
rc_kinet = 1.0
ra_therm = 1.0
rb_therm = 1.0
rc_therm = 1.0
j = 0.008
xeq = 0.01
d0 = 0.01
a0 = 200.0
b0 = 200.0
c0 = 200.0
vtot = n0 * a0 * b0 * c0 * 100
lich0 = 100
)
var (
ra_bal = 1.0 / 2
rb_bal = math.Sqrt(2)
rc_bal = math.Sqrt(2)
a = make([]float64, 0, n)
b = make([]float64, 0, n)
c = make([]float64, 0, n)
suma, sumb, sumc, sump, sumv float64
v0 float64
anew = make([]float64, n)
bnew = make([]float64, n)
cnew = make([]float64, n)
mean_a, mean_b, mean_c float64
mean_v float64
t = 0.0
d = d0
p float64
n = n0
lich = lich0
i int
)
func loadFloatFile(filename string, into *[]float64) {
f := lo.Must(os.Open(filename))
defer f.Close()
scanner := bufio.NewScanner(f)
for scanner.Scan() {
line := scanner.Text()
fields := strings.Fields(line)
for _, field := range fields {
val := lo.Must(strconv.ParseFloat(field, 64))
*into = append(*into, val)
}
}
lo.Must0(scanner.Err())
}
func initData() {
loadFloatFile("a.txt", &a)
loadFloatFile("b.txt", &b)
loadFloatFile("c.txt", &c)
for i := range a {
v0 += a[i] * b[i] * c[i]
}
}
func cleanFiles() {
os.Remove("mean_abc_go.txt")
os.Remove("fmean_v_go.txt")
os.Remove("fd_go.txt")
os.Remove("fn_go.txt")
os.Remove("fp_go.txt")
os.Remove("fv_go.txt")
}
func writeFmtToFile(filename string, format string, args ...interface{}) {
f := lo.Must(os.OpenFile(filename, os.O_APPEND|os.O_WRONLY|os.O_CREATE, 0644))
defer f.Close()
fmt.Fprintf(f, format, args...)
}
func saveData() {
lnT := math.Log(t)
writeFmtToFile("mean_abc_go.txt",
"%f %f %f %f\n",
lnT, math.Log(mean_a), math.Log(mean_b), math.Log(mean_c),
)
writeFmtToFile("fmean_v_go.txt", "%f %f\n", lnT, math.Log(mean_v))
writeFmtToFile("fd_go.txt", "%f %f\n", lnT, math.Log(d))
writeFmtToFile("fn_go.txt", "%f %f\n", lnT, math.Log(float64(n)))
writeFmtToFile("fp_go.txt", "%f %f\n", lnT, math.Log(p))
writeFmtToFile("fv_go.txt", "%f %f\n", lnT, math.Log(sumv))
}
func main() {
cleanFiles()
initData()
for t = dt; t < tmax; t += dt {
suma = 0
sumb = 0
sumc = 0
sumv = 0
sump = 0
for i = 0; i < n; {
anew[i] = a[i] + dt*ra_kinet*(d-(rc_therm/c[i])-(rb_therm/b[i])-(ra_bal/ra_kinet)*j)
bnew[i] = b[i] + dt*rb_kinet*(d-(ra_therm/a[i])-(rc_therm/c[i])-(rb_bal/rb_kinet)*j)
cnew[i] = c[i] + dt*rc_kinet*(d-(ra_therm/a[i])-(rb_therm/b[i])-(rc_bal/rc_kinet)*j)
if anew[i] < 0 || bnew[i] < 0 || cnew[i] < 0 {
a[i] = a[n-1]
b[i] = b[n-1]
c[i] = c[n-1]
n -= 1
} else {
suma += anew[i]
sumb += bnew[i]
sumc += cnew[i]
sumv += anew[i] * bnew[i] * cnew[i]
sump += anew[i]*bnew[i] + bnew[i]*cnew[i] + cnew[i]*anew[i]
i += 1
}
}
a, anew = anew, a
b, bnew = bnew, b
c, cnew = cnew, c
d = ((xeq+d0)*(1-v0/vtot)+v0/vtot-sumv/vtot)/(1-sumv/vtot) - xeq
p = 2 * sump
if lich == lich0 {
mean_a = suma / float64(n)
mean_b = sumb / float64(n)
mean_c = sumc / float64(n)
mean_v = sumv / float64(n)
saveData()
lich = 0
}
lich = lich + 1
}
}
func main() {
cleanFiles()
initData()
const_a := (ra_bal / ra_kinet) * j
const_b := (rb_bal / rb_kinet) * j
const_c := (rc_bal / rc_kinet) * j
log.Println(lo.Duration(func() {
for t = dt; t < tmax; t += dt {
suma = 0
sumb = 0
sumc = 0
sumv = 0
sump = 0
for i = 0; i < n; {
over_c := rc_therm / c[i]
over_b := rb_therm / b[i]
over_a := ra_therm / a[i]
anew[i] = a[i] + dt*ra_kinet*(d-over_c-over_b-const_a)
bnew[i] = b[i] + dt*rb_kinet*(d-over_a-over_c-const_b)
cnew[i] = c[i] + dt*rc_kinet*(d-over_a-over_b-const_c)
if anew[i] < 0 || bnew[i] < 0 || cnew[i] < 0 {
a[i] = a[n-1]
b[i] = b[n-1]
c[i] = c[n-1]
n -= 1
} else {
if lich == lich0 {
suma += anew[i]
sumb += bnew[i]
sumc += cnew[i]
sump += anew[i]*bnew[i] + bnew[i]*cnew[i] + cnew[i]*anew[i]
}
sumv += anew[i] * bnew[i] * cnew[i]
i += 1
}
}
a, anew = anew, a
b, bnew = bnew, b
c, cnew = cnew, c
d = ((xeq+d0)*(1-v0/vtot)+v0/vtot-sumv/vtot)/(1-sumv/vtot) - xeq
p = 2 * sump
if lich == lich0 {
mean_a = suma / float64(n)
mean_b = sumb / float64(n)
mean_c = sumc / float64(n)
mean_v = sumv / float64(n)
saveData()
lich = 0
}
lich = lich + 1
}
}))
}
func main() {
cleanFiles()
initData()
log.Println(lo.Duration(func() {
(func(a, b, c []float64) {
anew := make([]float64, n0)
bnew := make([]float64, n0)
cnew := make([]float64, n0)
const_a := (ra_bal / ra_kinet) * j
const_b := (rb_bal / rb_kinet) * j
const_c := (rc_bal / rc_kinet) * j
d := d0
n := n0
lich := lich0
for t := dt; t < tmax; t += dt {
var suma, sumb, sumc, sumv, sump float64
for i := 0; i < n; {
over_c := rc_therm / c[i]
over_b := rb_therm / b[i]
over_a := ra_therm / a[i]
anew[i] = a[i] + dt*ra_kinet*(d-over_c-over_b-const_a)
bnew[i] = b[i] + dt*rb_kinet*(d-over_a-over_c-const_b)
cnew[i] = c[i] + dt*rc_kinet*(d-over_a-over_b-const_c)
if anew[i] < 0 || bnew[i] < 0 || cnew[i] < 0 {
a[i] = a[n-1]
b[i] = b[n-1]
c[i] = c[n-1]
n -= 1
} else {
if lich == lich0 {
suma += anew[i]
sumb += bnew[i]
sumc += cnew[i]
sump += anew[i]*bnew[i] + bnew[i]*cnew[i] + cnew[i]*anew[i]
}
sumv += anew[i] * bnew[i] * cnew[i]
i += 1
}
}
a, anew = anew, a
b, bnew = bnew, b
c, cnew = cnew, c
d = ((xeq+d0)*(1-v0/vtot)+v0/vtot-sumv/vtot)/(1-sumv/vtot) - xeq
if lich == lich0 {
mean_a := suma / float64(n)
mean_b := sumb / float64(n)
mean_c := sumc / float64(n)
mean_v := sumv / float64(n)
p := 2 * sump
saveData(n, t, mean_a, mean_b, mean_c, mean_v, d, p, sumv)
lich = 0
}
lich = lich + 1
}
})(aa0, bb0, cc0)
}))
}
type workerState struct {
lich int
d float64
}
type workerOutput struct {
suma, sumb, sumc, sumv, sump float64
n int
}
func worker(a, b, c []float64, in <-chan workerState, out chan<- workerOutput) {
anew := make([]float64, n0)
bnew := make([]float64, n0)
cnew := make([]float64, n0)
n := len(a)
const_a := (ra_bal / ra_kinet) * j
const_b := (rb_bal / rb_kinet) * j
const_c := (rc_bal / rc_kinet) * j
for payload := range in {
var output workerOutput
for i := 0; i < n; {
over_c := rc_therm / c[i]
over_b := rb_therm / b[i]
over_a := ra_therm / a[i]
anew[i] = a[i] + dt*ra_kinet*(payload.d-over_c-over_b-const_a)
bnew[i] = b[i] + dt*rb_kinet*(payload.d-over_a-over_c-const_b)
cnew[i] = c[i] + dt*rc_kinet*(payload.d-over_a-over_b-const_c)
if anew[i] < 0 || bnew[i] < 0 || cnew[i] < 0 {
a[i] = a[n-1]
b[i] = b[n-1]
c[i] = c[n-1]
n -= 1
} else {
if payload.lich == lich0 {
output.suma += anew[i]
output.sumb += bnew[i]
output.sumc += cnew[i]
output.sump += anew[i]*bnew[i] + bnew[i]*cnew[i] + cnew[i]*anew[i]
}
output.sumv += anew[i] * bnew[i] * cnew[i]
i += 1
}
}
a, anew = anew, a
b, bnew = bnew, b
c, cnew = cnew, c
output.n = n
out <- output
}
close(out)
}
func main() {
numGoroutines := lo.Must(strconv.Atoi(os.Getenv("NUM_GOROUTINES")))
cleanFiles()
initData()
chunkSize := int(math.Ceil(float64(n0) / float64(numGoroutines)))
aChunks := lo.Chunk(aa0, chunkSize)
bChunks := lo.Chunk(bb0, chunkSize)
cChunks := lo.Chunk(cc0, chunkSize)
var inChans []chan workerState
var outChans []chan workerOutput
for i := range numGoroutines {
inChan := make(chan workerState, 1)
outChan := make(chan workerOutput, 1)
inChans = append(inChans, inChan)
outChans = append(outChans, outChan)
go worker(aChunks[i], bChunks[i], cChunks[i], inChan, outChan)
}
log.Println(lo.Duration(func() {
state := workerState{
lich: lich0,
d: d0,
}
for t := dt; t < tmax; t += dt {
for _, input := range inChans {
input <- state
}
var suma, sumb, sumc, sumv, sump float64
var n int
for _, outChan := range outChans {
output := <-outChan
n += output.n
suma += output.suma
sumb += output.sumb
sumc += output.sumc
sump += output.sump
sumv += output.sumv
}
state.d = ((xeq+d0)*(1-v0/vtot)+v0/vtot-sumv/vtot)/(1-sumv/vtot) - xeq
if state.lich == lich0 {
mean_a := suma / float64(n)
mean_b := sumb / float64(n)
mean_c := sumc / float64(n)
mean_v := sumv / float64(n)
p := 2 * sump
saveData(n, t, mean_a, mean_b, mean_c, mean_v, state.d, p, sumv)
state.lich = 0
}
state.lich += 1
}
}))
for _, inChan := range inChans {
close(inChan)
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment