forked from mallano/gofem
/
hydrost.go
135 lines (118 loc) · 3.23 KB
/
hydrost.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
// Copyright 2015 Dorival Pedroso and Raul Durand. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package fem
import (
"log"
"github.com/cpmech/gofem/inp"
"github.com/cpmech/gosl/chk"
"github.com/cpmech/gosl/la"
"github.com/cpmech/gosl/ode"
)
// HydroStatic computes water pressure (pl) and intrinsic liquid density (ρL)
// based on the following model
//
// ρL = ρL0 + Cl・pl thus dρL/dpl = Cl
//
// Z(z) = zmax + T・(z - zmax) with 0 ≤ T ≤ 1
// dZ = (z - zmax)・dT
// dpl = ρL(pl)・g・(-dZ)
// dpl = ρL(pl)・g・(zmax - z)・dT
// Δz = zmax - z
//
// / dpl/dT \ / ρL(pl)・g・Δz \
// dY/dT = | | = | |
// \ dρL/dT / \ Cl・dpl/dT /
//
type HydroStatic struct {
zwater float64
ρL0 float64
g float64
Cl float64
fcn ode.Cb_fcn
Jac ode.Cb_jac
sol ode.ODE
}
// Init initialises this structure
func (o *HydroStatic) Init() {
// basic data
o.zwater = Global.Sim.WaterLevel
o.ρL0 = Global.Sim.WaterRho0
o.g = Global.Sim.Gfcn.F(0, nil)
o.Cl = o.ρL0 / Global.Sim.WaterBulk
// x := {pl, ρL}
o.fcn = func(f []float64, x float64, y []float64, args ...interface{}) error {
Δz := args[0].(float64)
//ρL := o.ρL0
ρL := y[1]
f[0] = ρL * o.g * Δz // dpl/dT
f[1] = o.Cl * f[0] // dρL/dT
return nil
}
o.Jac = func(dfdy *la.Triplet, x float64, y []float64, args ...interface{}) error {
if dfdy.Max() == 0 {
dfdy.Init(2, 2, 4)
}
Δz := args[0].(float64)
dfdy.Start()
dfdy.Put(0, 0, 0)
dfdy.Put(0, 1, o.g*Δz)
dfdy.Put(1, 0, 0)
dfdy.Put(1, 1, o.Cl*o.g*Δz)
return nil
}
silent := true
o.sol.Init("Radau5", 2, o.fcn, o.Jac, nil, nil, silent)
o.sol.Distr = false // must be sure to disable this; otherwise it causes problems in parallel runs
}
// Calc computes pressure and density
func (o HydroStatic) Calc(z float64) (pl, ρL float64, err error) {
Δz := o.zwater - z
y := []float64{0, o.ρL0} // pl0, ρL0
err = o.sol.Solve(y, 0, 1, 1, false, Δz)
if err != nil {
err = chk.Err("HydroStatic failed when calculating pressure using ODE solver: %v", err)
return
}
return y[0], y[1], nil
}
// SetHydroSt sets the initial state to a hydrostatic condition
func (o *Domain) SetHydroSt(stg *inp.Stage) (ok bool) {
// set Sol
ndim := Global.Ndim
for _, n := range o.Nodes {
z := n.Vert.C[ndim-1]
dof := n.GetDof("pl")
if dof != nil {
pl, _, err := Global.HydroSt.Calc(z)
if LogErr(err, "hydrost: cannot compute pl") {
return
}
o.Sol.Y[dof.Eq] = pl
}
}
// set elements
var err error
for _, e := range o.ElemIntvars {
// build map with pressures @ ips
coords := e.Ipoints()
nip := len(coords)
pl := make([]float64, nip)
ρL := make([]float64, nip)
for i := 0; i < nip; i++ {
z := coords[i][ndim-1]
pl[i], ρL[i], err = Global.HydroSt.Calc(z)
if LogErr(err, "hydrost: cannot compute pl and ρL") {
return
}
}
ivs := map[string][]float64{"pl": pl, "ρL": ρL}
// set element's states
if LogErrCond(!e.SetIniIvs(o.Sol, ivs), "hydrost: element's internal values setting failed") {
return
}
}
// success
log.Printf("dom: initial hydrostatic state set")
return true
}