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

阿弥陀佛

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

 
 
 

日志

 
 
关于我

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

网易考拉推荐

边界Roe数值通量实现  

2016-09-30 14:48:19|  分类: Spark |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

Implementation

The quantities we then need to calculate in order to obtain the Roe numerical flux at the boundary, x_iph (i.e. x_iplushalf) are as follows:

  1. We need the values of the conservative variable
    q[x]
    to the left and right of the boundary. Methods: recos_qL, recos_qR, with the help of dqL_1, dqR_1, etc.
  2. The idea is to estimate the value of q at the cell boundary xiph based on the averages in the adjacent cells. The most intuitive way to do this is to estimate the values of the slope on either side and extrapolate, so as follows

    $\displaystyle q^L_{i+\frac{1}{2}} = q_i + \frac{q_i - q_{i-1}}{x_i - x_{i-1}}
$

    $\displaystyle q^R_{i+\frac{1}{2}} = q_{i+1} - \frac{q_{i+1} - q_i}{x_{i+1} - x_i}
$

    However, this can lead to problems when there are discontinuities. Since the slopes around discontinuities are very steep, estimated slopes can sometimes overshoot true cell averages. This can lead to too much total variation, defined as

    $\displaystyle \sum_i \vert q_{i+1} - q_i\vert
$

    A high total variation indicates large amounts of oscillation, which is not desirable in a solution and indicates instability. Various "slope-limiter" reconstruction methods aim to limit the variations so the solution method is total variation diminishing, or TVD. The idea is to restrict the ratio between adjacent jumps in cell averages where such jumps seem to indicate oscillatory behavior, while (depending on the method) steepening this jumps in other places to sharpen resolution around discontinuities that are actually part of the solution. This way, unwanted variation is reduced without sacrificing too much accuracy. There are options for four different types of slope limiters this code (called minmod, superbee, MC, and van Leer). The methods dqL_1, dqR_1, dqL_2, dqR_2, etc. calculate limited jumps for use in the recos_qL and recos_qR methods. By default the minmod limiter (dq1) is used.
    private double dqL_1(double dq_imh, double dq_iph) {
    		return minmod(dq_iph, dq_imh);
    	}
    
    //minmod function
    private double minmod(double a, double b){
    	if (a*b < 0) return 0;
    	else if(Math.abs(a) < Math.abs(b)) return a;
    	else return b;
    	}
    
    In practical terms, this means that if the slops on either side of the boundary have different signs, which could indicate spurious oscillation, the value will be reconstructed using a slope of zero. Similarly, if the slopes on either side have the same sign, the less steep one is always chosen. This clearly reduces the total variation of the solution, but it may slightly smooth out discontinuities even when they are desired in the solution.
  3. We need the jumps between right and left values, omega_i. Method: calc_omega.
  4. We need the right and left fluxes based on (1). Method: calc_f
  5. We need the Jacobian matrix at x_iph, along with its eigenvalues and eigenvectors lambda and eta. Methods: calc_A (for the Jacobian), calc_eta, calc_lambda.

    In the case of Burger's equation, the conservative variable and the flux functions are both 1D, so the Jacobian and single eigenvalue lambda are just the scalar q, and the eigenvector is just 1. This makes the above methods trivial, but they would be non-trivial in a multidimensional situation.

Time-stepping method

This program uses a 2nd order Runge-Kutta time stepping method. This involves estimating the derivative twice at each time step, rather than once as in the first-order Euler method

$\displaystyle \textbf q^{n+1} = \textbf q^n + \textbf v(\textbf q^n, t_n) \triangle t + O(t^2)$$\displaystyle \mbox { where } \textbf v(\textbf q, t) = \frac{\partial \textbf q}{\partial t} \Big \vert _t
$

The Runge-Kutta method uses the slope at the nth time step to estimate q at $ t_{n+\frac{1}{2}}$ , and based on the value $ q^{n+\frac{1}{2}}$ the "midpoint" slope at the $ \frac{\partial \textbf q}{\partial t} \Big \vert _{t_{n+\frac{1}{2}}}$ level is obtained. This slope is used in place of $ \frac{\partial \textbf q}{\partial t} \Big \vert _{t_n}$ , and leads to a higher order of accuracy

$\displaystyle \tilde{\textbf q}^{n+\frac{1}{2}} = \textbf q^n +\textbf v(\textbf q^n, t_n) \frac{\triangle t}{2}
$

$\displaystyle \textbf q^{n+1} = \textbf q^n + \textbf v(\tilde{\textbf q}^{n + \frac{1}{2}}, t_{n+\frac{1}{2}}) \triangle t + O(t^3)
$

In this case,

$\displaystyle v(\textbf q, t) = \textbf F(
$

This is implemented as follows in the step() method.
//moving forward one time step
private double[] step(double[] q) {

	//half advanced time step to estimate v at the midpoint
	double[] FF = calc_flux(q, x, Nx, Ng);
	double[] q_nph = update_q(q, x, FF, delta_t/2, Nx, Ng);

	//full advanced time step
	FF = calc_flux(q_nph, x, Nx, Ng);
	double[] q_np1 = update_q(q_nph, x, FF, delta_t, Nx, Ng);

	return q_np1;

	}
where in each case updating uses the derivative v based on the fluxes as described above (note that FF[i] contains the flux at the boundary x_iph).
for (int i = Ng; i < Nx - Ng; i++) {
	qnp1[i] = qn[i] - delta_t/dx*(FF[i] - FF[i-1]);
}

  评论这张
 
阅读(42)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

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

页脚

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