注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

阿弥陀佛

街树飘影未见尘 潭月潜水了无声 般若观照心空静...

 
 
 

日志

 
 
关于我

一直从事气象预报、服务建模实践应用。 注重气象物理场、实况场、地理信息、本体知识库、分布式气象内容管理系统建立。 对Barnes客观分析, 小波,计算神经网络、信任传播、贝叶斯推理、专家系统、网络本体语言有一定体会。 一直使用Java、Delphi、Prolog、SQL编程。

网易考拉推荐

Burgers  

2016-09-27 17:19:24|  分类: Scala |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
package stat

import stddraw.StdDraw

/**
* Created by hxf on 16-9-26.
*/
case class
Burgers(smooth: Boolean, //toggles Gaussian or shock tube initial data
delta_t: Double, //length of time steps
Nx: Int, //number of cells
Ng: Int, //number of ghost cells at each boundary
xmin: Double, //domain
xmax: Double,
delta_x: Double, //distance between cells
x: Array[Double] //holds positions of cell centers
)

object Burgers {
var numsteps: Int = 0

//number of time steps
def apply(Nx: Int, delta_t: Double, smooth: Boolean): Burgers = {
val xmin = 0.0
val xmax = 5.0
this (xmin, xmax, Nx, delta_t, smooth)
}

def apply(xmin: Double, xmax: Double, Nx: Int, delta_t: Double, smooth: Boolean): Burgers = {
val Ng = 2
val delta_x = (xmax - xmin) / Nx
val x = Array.tabulate(Nx)(i => xmin + (i + 0.5) * delta_x)
Burgers(smooth, delta_t, Nx, Ng, xmin, xmax, delta_x, x)
}

import scala.math._

def init_fluid(b: Burgers): Array[Double] = {
val q = Array.ofDim[Double](b.Nx)
b.smooth match {
case true => Array.tabulate(b.Nx)(i =>
q(i) = exp(-(b.x(i) - (b.xmax + b.xmin) / 2.0) * (b.x(i) - (b.xmax + b.xmin) / 2.0))
)
case false =>
val q_L = 0.1
val q_R = 1.0
Array.tabulate(b.Nx) { i =>
q(i) = if (i < b.Nx / 2) q_L else q_R
}
}
q
}

def minmod(a: Double, b: Double): Double = {
if (a * b < 0) 0.0
else if (abs(a) < abs(b)) a
else b
}

//calculating Jacobian matrix at x
private def calc_A(q: Double) = {
//the Jacobian is just q
val A = q
A
}

//calculating eigenvectors of Jacobian
private def calc_eta(A: Double): Double = {
//just the same as A in the scalar case
val lambda = A
lambda
}

//calculating eigenvalues of Jacobian
private def calc_lambda(A: Double): Double = {
//just the same as A in the scalar case
val lambda: Double = A
lambda
}

//calculating jumps
private def calc_omega(qL: Double, qR: Double, eta: Double): Double = {
val omega = qR - qL
omega
}

//calculate Roe flux
def calc_FF(fR: Double, fL: Double, lambda: Double, eta: Double, omega: Double): Double = {
val FF = 0.5 * (fR + fL - lambda * eta * omega)
FF
}

//flux as a function of q
def calc_f(q: Double) = {
val f = 0.5 * q * q
f
}

def dqL_1(dq_imh: Double, dq_iph: Double): Double = minmod(dq_iph, dq_imh)

def dqR_1(dq_iph: Double, dq_ip3h: Double): Double = minmod(dq_iph, dq_ip3h)

def dqL_2(dq_imh: Double, dq_iph: Double): Double = {
if (abs(dq_iph) <= 1e-6) 0
else {
val theta = dq_imh / dq_iph
max(max(0, min(1, 2 * theta)), min(2, theta)) * dq_iph
}
}

def dqR_2(dq_iph: Double, dq_ip3h: Double): Double = {
if (abs(dq_iph) <= 1e-6) 0
else {
val theta = dq_ip3h / dq_iph
max(max(0, min(1, 2 * theta)), min(2, theta)) * dq_iph
}
}

def dqL_3(dq_imh: Double, dq_iph: Double): Double = {
if (abs(dq_iph) <= 1e-6) 0
else {
val theta = dq_imh / dq_iph
max(0, min(min((1 + theta) / 2, 2), 2 * theta)) * dq_iph
}
}

def dqR_3(dq_iph: Double, dq_ip3h: Double): Double = {
if (abs(dq_iph) <= 1e-6) 0
else {
val theta = dq_ip3h / dq_iph
max(0, min(min((1 + theta) / 2, 2), 2 * theta)) * dq_iph
}
}

def dqL_4(s_imh: Double, s_iph: Double): Double = {
if (abs(s_iph) <= 1e-6) 0
else {
val theta = s_imh / s_iph
s_iph * (theta + abs(theta)) / (1 + abs(theta))
}
}

private def dqR_4(s_iph: Double, s_ip3h: Double): Double = {
if (abs(s_iph) <= 1e-6) 0
else {
val theta = s_ip3h / s_iph
s_iph * (theta + abs(theta)) / (1 + abs(theta))
}
}

//reconstructs primitive variable q from the left
private def recos_qL(q_im1: Double, q_i: Double, q_ip1: Double, r_im1: Double, r_i: Double, r_ip1: Double): Double = {
val dq_iph: Double = q_ip1 - q_i
val dq_imh: Double = q_i - q_im1
val dq_i: Double = dqL_2(dq_imh, dq_iph)
val rqL_iph: Double = q_i + dq_i * 0.5
rqL_iph
}

//reconstructs primitive variable q from the right
private def recos_qR(q_i: Double, q_ip1: Double, q_ip2: Double, r_i: Double, r_ip1: Double, r_ip2: Double): Double = {
val dq_iph: Double = q_ip1 - q_i
val dq_ip3h: Double = q_ip2 - q_ip1
val dq_i: Double = dqR_2(dq_iph, dq_ip3h)
val rqR_iph: Double = q_ip1 - dq_i * 0.5
rqR_iph
}

//calculating Roe fluxes at each cell boundary (not including ghost cells)
private def calc_flux(q: Array[Double], x: Array[Double], Nx: Int, Ng: Int): Array[Double] = {
val FF = Array.ofDim[Double](Nx)
for (i <- Ng - 1 until Nx - Ng) {
val rqL_iph: Double = recos_qL(q(i - 1), q(i), q(i + 1), x(i - 1), x(i), x(i + 1))
val rqR_iph: Double = recos_qR(q(i), q(i + 1), q(i + 2), x(i), x(i + 1), x(i + 2))
val rq_iph: Double = 0.5 * (rqL_iph + rqR_iph)
val A: Double = calc_A(rq_iph)
val eta: Double = calc_eta(A)
val lambda: Double = calc_lambda(A)
val omega: Double = calc_omega(rqL_iph, rqR_iph, eta)
val fR_iph: Double = calc_f(rqR_iph)
val fL_iph: Double = calc_f(rqL_iph)
val FF_iph: Double = calc_FF(fR_iph, fL_iph, eta, lambda, omega)
FF(i) = FF_iph
}
FF
}

//updating
private def update_q(qn: Array[Double], x: Array[Double], FF: Array[Double], delta_t: Double, Nx: Int, Ng: Int): Array[Double] = {
val dx: Double = x(1) - x(0)
val qnp1 = Array.ofDim[Double](Nx)
for (i <- Ng until Nx - Ng) {
qnp1(i) = qn(i) - delta_t / dx * (FF(i) - FF(i - 1))
}
//constant extrapolation
for (i <- 0 until Ng) {
qnp1(i) = qnp1(Ng)
qnp1(Nx - i - 1) = qnp1(Nx - Ng - 1)
}
/*Either zero or first order extrapolation must be commented out
//first order extrapolation
for (int i = 0; i < Ng; i++) {
qnp1[i] = qnp1[Ng] - (Ng - i)*(qnp1[Ng + 1] - qnp1[Ng]);
qnp1[Nx - i - 1] = qnp1[Nx - Ng - 1] + (Ng - i)*(qnp1[Nx - Ng - 1] - qnp1[Nx - Ng - 2]);
}
*/
qnp1
}

//moving forward one time step
private def step(b: Burgers, q: Array[Double]): Array[Double] = {
//half advanced time step
var FF: Array[Double] = calc_flux(q, b.x, b.Nx, b.Ng)
val q_nph: Array[Double] = update_q(q, b.x, FF, b.delta_t / 2, b.Nx, b.Ng)
//full advanced time step
FF = calc_flux(q_nph, b.x, b.Nx, b.Ng)
val q_np1: Array[Double] = update_q(q_nph, b.x, FF, b.delta_t, b.Nx, b.Ng)
q_np1
}

//animates evolution until time tmax
def animate(b: Burgers, tmax: Double) {
StdDraw.setXscale(b.xmin, b.xmax)
val q = init_fluid(b)
numsteps = Math.floor(tmax / b.delta_t).toInt
var curr: Array[Double] = q
for (i <- 0 until numsteps) {
StdDraw.clear()
var j: Int = 0
for (j <- 0 until b.Nx - 1) {
StdDraw.line(b.x(j), curr(j), b.x(j + 1), curr(j + 1))
}
//drawing axes
StdDraw.line(b.xmin, 0, b.xmax, 0)
StdDraw.line(0, 0, 0, 1)
for (k <- b.xmin.toInt to b.xmax.toInt) {
if (k == 0) {
}
else {
StdDraw.line(k.toDouble, 0, k.toDouble, 0.01)
val I: String = "" + k
StdDraw.text(k.toDouble, 0.025, I)
}
}
for (k <- 1 to 5) {
StdDraw.line(0, 0.2 * k, 0.1, 0.2 * k)
val I: String = "%.1f".format(0.2f * k)
StdDraw.text(0.35, 0.2 * k, I)
}
StdDraw.show(30)
curr = step(b, curr)
}
}

//displays the solution at a time tmax
def time(b: Burgers, tmax: Double) {
val q = init_fluid(b)
numsteps = Math.floor(tmax / b.delta_t).toInt
var curr: Array[Double] = q
for (i <- 0 until numsteps) {
curr = step(b, curr)
}
StdDraw.setXscale(b.xmin, b.xmax)
var j: Int = 0
for (j <- 0 until b.Nx - 1) {
StdDraw.line(b.x(j), curr(j), b.x(j + 1), curr(j + 1))
}
//drawing axes
StdDraw.line(b.xmin, 0, b.xmax, 0)
StdDraw.line(0, 0, 0, 1)
var i: Int = b.xmin.toInt
for (i <- b.xmin.toInt to b.xmax.toInt) {
if (i == 0) {
}
else {
StdDraw.line(i.toDouble, 0, i.toDouble, 0.01)
val I: String = "" + i
StdDraw.text(i.toDouble, 0.025, I)
}
}
for (i <- 0 to 5) {
StdDraw.line(0, 0.2 * i, 0.1, 0.2 * i)
val I: String = "%.1f".format(0.2f * i)
StdDraw.text(0.35, 0.2 * i, I)
}
}

//returns values of q at each mesh point after time tmax
def returnq(b: Burgers, tmax: Double): Array[Double] = {
val q = init_fluid(b)
numsteps = Math.floor(tmax / b.delta_t).toInt
var curr: Array[Double] = q
for (i <- 0 until numsteps) {
curr = step(b, curr)
}
curr
}

//graphs the convergence factor as a function of x
def convergence(Nx: Int, time: Double) {
val dt: Double = 5.0 / (Nx * 20)
val one: Burgers = Burgers(Nx, dt, true)
val two: Burgers = Burgers(2 * Nx, dt, true)
val four: Burgers = Burgers(4 * Nx, dt, true)
val One: Array[Double] = returnq(one, time)
val Two: Array[Double] = returnq(two, time)
val Four: Array[Double] = returnq(four, time)
val Q = Array.ofDim[Double](Nx)
StdDraw.setXscale(one.xmin, one.xmax)
StdDraw.setYscale(0, 10)
//drawing axes
StdDraw.line(one.xmin, 0, one.xmax, 0)
StdDraw.line(0, 0, 0, 10)
for (i <- one.xmin.toInt to one.xmax.toInt) {
if (i == 0) {
}
else {
StdDraw.line(i.toDouble, 0, i.toDouble, 0.1)
val I: String = "" + i
StdDraw.text(i.toDouble, 0.25, I)
}
}
for (i <- 0 until 10) {
if (i == 0) {
}
else {
StdDraw.line(0, i.toDouble, 0.1, i.toDouble)
val I: String = "" + i
StdDraw.text(0.25, i.toDouble, I)
}
}
StdDraw.setPenRadius(0.005)
for (i <- 0 until Nx) {
val e1: Double = Math.abs(One(i) - ((Two(2 * i) + Two((2 * i) + 1)) / 2))
val e2: Double = Math.abs(((Two(2 * i) + Two((2 * i) + 1)) / 2) - (Four(4 * i) + Four(4 * i + 1) + Four(4 * i + 2) + Four(4 * i + 3)) / 4)
Q(i) = e1 / e2
StdDraw.point(one.x(i), e1 / e2)
}
}

def main(args: Array[String]): Unit = {
val hi = Burgers(-2, 2, 200, .01, true)
animate(hi,2)
}
}
  评论这张
 
阅读(39)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017