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

阿弥陀佛

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

 
 
 

日志

 
 
关于我

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

网易考拉推荐

UltraRel  

2016-10-12 17:25:20|  分类: Scala |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
package burger

import stddraw.StdDraw

/**
* Created by hxf on 16-9-27.
*/

object UltraRel {
//number of time steps
var numsteps: Int = 0
//adiabatic index
val k: Double = 1.5
//conservative variables must always be above floor
val floor: Double = 1e-11

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[Array[Double]] = {
val vars = Array.ofDim[Double](4, b.Nx)
//if smooth = true, Gaussian initial data, otherwise shock tube
if (b.smooth)
Array.tabulate(b.Nx)(i => vars(0)(i) = exp(-(b.x(i) - 2.5) * (b.x(i) - 2.5)))
else {
//shock tube
val rho_L = 1.0
val rho_R = 0.1
Array.tabulate(b.Nx) { i =>
if (i < b.Nx / 2) vars(0)(i) = rho_L else vars(0)(i) = rho_R
}
}
val v_L = 0.0
val v_R = 0.0
Array.tabulate(b.Nx) { i =>
if (i < b.Nx / 2.0) vars(1)(i) = v_L else vars(1)(i) = v_R
}
val P = Array.tabulate(b.Nx)(i => (k - 1) * vars(0)(i))
Array.tabulate(b.Nx) { i =>
vars(2)(i) = (vars(0)(i) + P(i)) / (1 - vars(1)(i) * vars(1)(i)) - P(i)
}
Array.tabulate(b.Nx) { i =>
vars(3)(i) = vars(1)(i) * (vars(2)(i) + P(i)) / (1 - vars(1)(i) * vars(1)(i))
}
vars
}

//minmod function
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(tau: Double, S: Double, rho: Double, v: Double): Array[Array[Double]] = {
val beta: Double = (2 - k) / 4
val P_tau: Double = -2 * beta + (4 * beta * beta + k - 1) * tau / Math.sqrt(4 * beta * beta * tau * tau + (k - 1) * (tau * tau - S * S))
val P_s: Double = -(k - 1) * S / Math.sqrt(4 * beta * beta * tau * tau + (k - 1) * (tau * tau - S * S))
val A = Array.ofDim[Double](2, 2)
A(0)(0) = 0
A(0)(1) = 1
A(1)(0) = -v * v + (1 - v * v) * P_tau
A(1)(1) = 2 * v + (1 - v * v) * P_s
A
}

//calculating the first eigenvector of Jacobian
private def calc_eta1(A: Array[Array[Double]]): Array[Double] = {
val eta1 = Array.ofDim[Double](2)
val lambda1: Double = (A(0)(0) + A(1)(1) + Math.sqrt((A(0)(0) - A(1)(1)) * (A(0)(0) - A(1)(1)) + 4 * A(0)(1) * A(1)(0))) / 2
eta1(0) = 1
eta1(1) = (lambda1 - A(0)(0)) / A(0)(1)
eta1
}

//calculating the other eigenvector
private def calc_eta2(A: Array[Array[Double]]): Array[Double] = {
val eta2 = Array.ofDim[Double](2)
val lambda2: Double = (A(0)(0) + A(1)(1) - Math.sqrt((A(0)(0) - A(1)(1)) * (A(0)(0) - A(1)(1)) + 4 * A(0)(1) * A(1)(0))) / 2
eta2(0) = 1
eta2(1) = (lambda2 - A(0)(0)) / A(0)(1)
eta2
}

//returns both eigenvalues
private def calc_lambda(A: Array[Array[Double]]): Array[Double] = {
val lambda = Array.ofDim[Double](2)
lambda(0) = (A(0)(0) + A(1)(1) + Math.sqrt((A(0)(0) - A(1)(1)) * (A(0)(0) - A(1)(1)) + 4 * A(0)(1) * A(1)(0))) / 2
lambda(1) = (A(0)(0) + A(1)(1) - Math.sqrt((A(0)(0) - A(1)(1)) * (A(0)(0) - A(1)(1)) + 4 * A(0)(1) * A(1)(0))) / 2
lambda
}

//calculates jumps in both conservative variables
private def calc_omega(tauL: Double, sL: Double, tauR: Double, sR: Double, eta1: Array[Double], eta2: Array[Double]): Array[Double] = {
val omega = Array.ofDim[Double](2)
omega(0) = (eta2(1) * (tauR - tauL) + eta2(0) * (sL - sR)) / (eta1(0) * eta2(1) - eta1(1) * eta2(0))
omega(1) = (eta1(0) * (sR - sL) + eta1(1) * (tauL - tauR)) / (eta1(0) * eta2(1) - eta1(1) * eta2(0))
omega
}

//calculate Roe flux
private def calc_FF(fL: Array[Double], fR: Array[Double], lambda: Array[Double], eta1: Array[Double], eta2: Array[Double], omega: Array[Double]): Array[Double] = {
val FF = Array.ofDim[Double](2)
FF(0) = 0.5 * (fR(0) + fL(0) - Math.abs(lambda(0)) * eta1(0) * omega(0) - Math.abs(lambda(1)) * eta2(0) * omega(1))
FF(1) = 0.5 * (fR(1) + fL(1) - Math.abs(lambda(0)) * eta1(1) * omega(0) - Math.abs(lambda(1)) * eta2(1) * omega(1))
FF
}

//calculates flux in terms of conservative variables
private def calc_f(tau: Double, S: Double, rho: Double, v: Double): Array[Double] = {
val f = Array.ofDim[Double](2)
val P: Double = (k - 1) * rho
f(0) = S
f(1) = S * v + P
f
}

//calculates conservative variables from primitive variables
private def calc_cons(rho: Double, v: Double): Array[Double] = {
val cons = Array.ofDim[Double](2)
val P: Double = (k - 1) * rho
cons(0) = (rho + P) / (1 - v * v) - P
cons(1) = v * (cons(0) + P) / (1 - v * v)
cons(0) = Math.max(cons(0), floor + Math.abs(cons(1)))
cons
}

//calculates primitive variables from conservative variables
private def calc_prim(tau: Double, S: Double): Array[Double] = {
val prim = Array.ofDim[Double](2)
val beta: Double = (2 - k) / 4
val P: Double = -2 * beta * tau + Math.sqrt(4 * beta * beta * tau * tau + (k - 1) * (tau * tau - S * S))
prim(0) = P / (k - 1)
prim(1) = S / (tau + P)
prim
}

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

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

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

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

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

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

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

private def dqR_4(s_iph: Double, s_ip3h: Double): Double = if (Math.abs(s_iph) <= 1e-6) 0
else {
val theta: Double = s_ip3h / s_iph
s_iph * (theta + Math.abs(theta)) / (1 + Math.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_1(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_1(dq_iph, dq_ip3h)
val rqR_iph: Double = q_ip1 - dq_i * 0.5
rqR_iph
}

//calculates Roe fluxes at each cell boundary (not including ghost cells)
private def calc_flux(vars: Array[Array[Double]], x: Array[Double], Nx: Int, Ng: Int): Array[Array[Double]] = {
val FF = Array.ofDim[Double](2, Nx)
var i: Int = Ng - 1
for (i <- Ng - 1 until Nx - Ng) {
//reconstructing primitive variables
val rrhoL_iph: Double = recos_qL(vars(0)(i - 1), vars(0)(i), vars(0)(i + 1), x(i - 1), x(i), x(i + 1))
val rrhoR_iph: Double = recos_qR(vars(0)(i), vars(0)(i + 1), vars(0)(i + 2), x(i), x(i + 1), x(i + 2))
val rrho_iph: Double = 0.5 * (rrhoL_iph + rrhoR_iph)
val rvL_iph: Double = recos_qL(vars(1)(i - 1), vars(1)(i), vars(1)(i + 1), x(i - 1), x(i), x(i + 1))
val rvR_iph: Double = recos_qR(vars(1)(i), vars(1)(i + 1), vars(1)(i + 2), x(i), x(i + 1), x(i + 2))
val rv_iph: Double = 0.5 * (rvL_iph + rvR_iph)
//calculating conservative variables from reconstructed values
val rqL_iph: Array[Double] = calc_cons(rrhoL_iph, rvL_iph)
val rqR_iph: Array[Double] = calc_cons(rrhoR_iph, rvR_iph)
val rq_iph: Array[Double] = calc_cons(rrho_iph, rv_iph)
//calculating Jacobian
val A: Array[Array[Double]] = calc_A(rq_iph(0), rq_iph(1), rrho_iph, rv_iph)
val eta1: Array[Double] = calc_eta1(A)
val eta2: Array[Double] = calc_eta2(A)
val lambda: Array[Double] = calc_lambda(A)
val omega: Array[Double] = calc_omega(rqL_iph(0), rqL_iph(1), rqR_iph(0), rqR_iph(1), eta1, eta2)
val fR_iph: Array[Double] = calc_f(rqR_iph(0), rqR_iph(1), rrhoR_iph, rvR_iph)
val fL_iph: Array[Double] = calc_f(rqL_iph(0), rqL_iph(1), rrhoL_iph, rvL_iph)
val FF_iph: Array[Double] = calc_FF(fR_iph, fL_iph, lambda, eta1, eta2, omega)
FF(0)(i) = FF_iph(0)
FF(1)(i) = FF_iph(1)
}
FF
}

//updating
private def update_q(q_n: Array[Array[Double]], x: Array[Double], FF: Array[Array[Double]], delta_t: Double, Nx: Int, Ng: Int): Array[Array[Double]] = {
val dx: Double = x(1) - x(0)
val q_np1 = Array.ofDim[Double](4, Nx)
for (i <- Ng until Nx - Ng) {
q_np1(2)(i) = q_n(2)(i) - delta_t * (FF(0)(i) - FF(0)(i - 1)) / dx
q_np1(3)(i) = q_n(3)(i) - delta_t * (FF(1)(i) - FF(1)(i - 1)) / dx
q_np1(2)(i) = Math.max(q_np1(2)(i), Math.abs(q_np1(3)(i)) + floor)
}
// //constant extrapolation
// for (int i = 0; i < Ng; i++) {
// q_np1[2][i] = q_np1[2][Ng];
// q_np1[2][Nx - i - 1] = q_np1[2][Nx - Ng - 1];
// q_np1[3][i] = q_np1[3][Ng];
// q_np1[3][Nx - i - 1] = q_np1[3][Nx - Ng - 1];
// }
//first order extrapolation
for (i <- 0 until Ng) {
q_np1(2)(i) = q_np1(2)(Ng) - (Ng - i) * (q_np1(2)(Ng + 1) - q_np1(2)(Ng))
q_np1(2)(Nx - i - 1) = q_np1(2)(Nx - Ng - 1) + (Ng - i) * (q_np1(2)(Nx - Ng - 1) - q_np1(2)(Nx - Ng - 2))
q_np1(3)(i) = q_np1(3)(Ng) - (Ng - i) * (q_np1(3)(Ng + 1) - q_np1(3)(Ng))
q_np1(3)(Nx - i - 1) = q_np1(3)(Nx - Ng - 1) + (Ng - i) * (q_np1(3)(Nx - Ng - 1) - q_np1(3)(Nx - Ng - 2))
}
for (i <- 0 until Nx) {
q_np1(0)(i) = calc_prim(q_np1(2)(i), q_np1(3)(i))(0)
q_np1(1)(i) = calc_prim(q_np1(2)(i), q_np1(3)(i))(1)
}
q_np1
}

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

//shows the solution evolving until time tmax
def animate(b: Burgers, tmax: Double) {
StdDraw.setXscale(b.xmin, b.xmax)
val vars = init_fluid(b)
numsteps = Math.floor(tmax / b.delta_t).toInt
var curr: Array[Array[Double]] = vars
for (i <- 0 until numsteps) {
StdDraw.clear()
for (j <- 0 until b.Nx - 1) {
StdDraw.line(b.x(j), curr(0)(j), b.x(j + 1), curr(0)(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.2 * k) //String.format("%.1g%n", 0.2 * k)
StdDraw.text(0.35, 0.2 * k, I)
}
// val filename: String = "UltraRel" + i + ".png"
// StdDraw.save(filename)
StdDraw.show(30)
curr = step(b, curr)
}
}

//displays the solution at a time tmax
def time(b: Burgers, tmax: Double) {
val vars = init_fluid(b)
numsteps = Math.floor(tmax / b.delta_t).toInt
var curr: Array[Array[Double]] = vars
for (i <- 0 until numsteps) {
curr = step(b, curr)
}
StdDraw.setXscale(b.xmin, b.xmax)
for (j <- 0 until b.Nx - 1) {
StdDraw.line(b.x(j), curr(0)(j), b.x(j + 1), curr(0)(j + 1))
}

//drawing axes
StdDraw.line(b.xmin, 0, b.xmax, 0)
StdDraw.line(0, 0, 0, 1)
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 <- 1 to 5) {
StdDraw.line(0, 0.2 * i, 0.1, 0.2 * i)
val I: String = "%.1f".format(0.2f * i) //String.format("%.1g%n", 0.2 * i)
StdDraw.text(0.35, 0.2 * i, I)
}
}

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

//takes the L2 measure of difference between two arrays a and b of the same length,
//w/o contributions from j cells at each end
def norm(a: Array[Double], b: Array[Double], j: Int): Double = {
var c: Double = 0
var i: Int = j
for (i <- j until a.length - j) {
c += (a(i) - b(i)) * (a(i) - b(i))
}
val norm: Double = Math.sqrt(c / a.length)
norm
}

//given a starting number of mesh points Nx and Gaussian initial data, calculates "convergence factor"
//norm(rho(Nx) - rho(2*Nx))/norm(rho(2*Nx) - rho(4*Nx)) - should be 2^n for an nth order method.
def test(Nx: Int, time: Double): Double = {
val dt: Double = 5.0 / (Nx * 20)
val one = UltraRel(Nx, dt, true)
val two = UltraRel(2 * Nx, dt, true)
val four = UltraRel(4 * Nx, dt, true)
val One: Array[Double] = returnrho(one, time)
val Two1: Array[Double] = Array.ofDim[Double](Nx)
val Two2: Array[Double] = Array.ofDim[Double](2 * Nx)
val Four: Array[Double] = Array.ofDim[Double](2 * Nx)
for (i <- 0 until Nx) {
Two1(i) = (returnrho(two, time)(2 * i) + returnrho(two, time)((2 * i) + 1)) / 2
Two2(i) = returnrho(two, time)(i)
Four(i) = (returnrho(four, time)(2 * i) + returnrho(four, time)((2 * i) + 1)) / 2
}
val Q: Double = norm(One, Two1, 0) / norm(Two2, Four, 0)
Q
}

def convergence(Nx: Int, time: Double) {
val dt: Double = 5.0 / (Nx * 20)
val one = UltraRel(Nx, dt, true)
val two = UltraRel(2 * Nx, dt, true)
val four = UltraRel(4 * Nx, dt, true)
val One: Array[Double] = returnrho(one, time)
val Two: Array[Double] = returnrho(two, time)
val Four: Array[Double] = returnrho(four, time)
val Q: Array[Double] = new Array[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 to 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]) {
val hi = UltraRel(200, .01, true)
animate(hi, 3)
}
}
  评论这张
 
阅读(14)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

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

页脚

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