CF461B Appleman and Tree (树DP)

本文解析Codeforces Round #263 (Div. 2) D题——Appleman and Tree,采用树形DP方法解决如何切割树的问题,确保每个连通块恰好有一个黑色节点。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

CF462D

Codeforces Round #263 (Div. 2) D

Codeforces Round #263 (Div. 1) B

B. Appleman and Tree
time limit per test
2 seconds
memory limit per test
256 megabytes
input
standard input
output
standard output

Appleman has a tree with n vertices. Some of the vertices (at least one) are colored black and other vertices are colored white.

Consider a set consisting of k (0 ≤ k < n) edges of Appleman's tree. If Appleman deletes these edges from the tree, then it will split into (k + 1) parts. Note, that each part will be a tree with colored vertices.

Now Appleman wonders, what is the number of sets splitting the tree in such a way that each resulting part will have exactly one black vertex? Find this number modulo 1000000007 (109 + 7).

Input

The first line contains an integer n (2  ≤ n ≤ 105) — the number of tree vertices.

The second line contains the description of the tree: n - 1 integers p0, p1, ..., pn - 2 (0 ≤ pi ≤ i). Where pi means that there is an edge connecting vertex (i + 1) of the tree and vertex pi. Consider tree vertices are numbered from 0 to n - 1.

The third line contains the description of the colors of the vertices: n integers x0, x1, ..., xn - 1 (xi is either 0 or 1). If xi is equal to 1, vertex i is colored black. Otherwise, vertex i is colored white.

Output

Output a single integer — the number of ways to split the tree modulo 1000000007 (109 + 7).

Sample test(s)
Input
3
0 0
0 1 1
Output
2
Input
6
0 1 1 0 4
1 1 0 0 1 0
Output
1
Input
10
0 1 2 1 4 4 4 0 8
0 0 0 1 0 1 1 0 0 1
Output
27

题意:有n个结点的树,结点编号0~n-1,0为根,分别给出1~n-1的父亲,再给出0~n-1各个结点的颜色(0为白,1为黑),要将其中一些边切掉,使每个联通块有且只有1个黑点,求切法种类数。

题解:树形DP。

从根DFS,f[x][0]表示对{x点和它的子树、x点连接父亲结点的边}这一整坨,有多少种方案使得x这个联通块没黑点(x是黑点的时候这个也不为0,是把x的父边切掉的种类数)

f[x][1]是这个联通块有黑点的种类数。

 

太难了!怪不得大家都掉分飞起,虽然题解的代码看起来很短,我根本想不出来啊看了半天还是不懂啊!

具体还是看代码吧,写了点注释,这个统计方法太碉了,我也弄得不是很清楚,算了日后再说。

代码:

 1 //#pragma comment(linker, "/STACK:102400000,102400000")
 2 #include<cstdio>
 3 #include<cmath>
 4 #include<iostream>
 5 #include<cstring>
 6 #include<algorithm>
 7 #include<cmath>
 8 #include<map>
 9 #include<set>
10 #include<stack>
11 #include<queue>
12 using namespace std;
13 #define ll long long
14 #define usll unsigned ll
15 #define mz(array) memset(array, 0, sizeof(array))
16 #define minf(array) memset(array, 0x3f, sizeof(array))
17 #define REP(i,n) for(i=0;i<(n);i++)
18 #define FOR(i,x,n) for(i=(x);i<=(n);i++)
19 #define RD(x) scanf("%d",&x)
20 #define RD2(x,y) scanf("%d%d",&x,&y)
21 #define RD3(x,y,z) scanf("%d%d%d",&x,&y,&z)
22 #define WN(x) printf("%d\n",x);
23 #define RE  freopen("D.in","r",stdin)
24 #define WE  freopen("1biao.out","w",stdout)
25 #define mp make_pair
26 #define pb push_back
27 
28 const int maxn=111111;
29 const int MOD=1e9+7;
30 
31 int n;
32 int a[maxn];
33 
34 struct Edge {
35     int next,v;
36 } e[2*maxn];
37 int en=0;
38 int head[maxn];
39 
40 void add(int x,int y) {
41     e[en].v=y;
42     e[en].next=head[x];
43     head[x]=en++;
44 }
45 
46 bool u[maxn];
47 ll f[maxn][2];///f[x][j] j=1表示x所在联通块有黑点,0表示无黑店 的种类数,包括x连接父亲的边和子树所有的边
48 void dfs(int x){
49     //printf("[in %d]",x);
50     int i;
51     u[x]=1;
52     f[x][0]=1;
53     f[x][1]=0;///先假设当前点是个白点
54     for(i=head[x]; i!=-1; i=e[i].next) {
55         if(!u[e[i].v]) {
56             dfs(e[i].v);
57             f[x][1]=(f[x][1]*f[e[i].v][0] + f[x][0]*f[e[i].v][1])%MOD;///有黑点的情况,先用已经统计的有黑点的情况乘一发儿子没黑点的情况,然后用已经统计的没黑点的情况乘一发儿子有黑点的情况
58             f[x][0]=f[x][0]*f[e[i].v][0]%MOD;///没黑点的情况直接乘儿子没黑点的情况
59         }
60     }
61     u[x]=0;
62     ///下面是对x点的父边的处理
63     if(a[x]==0)f[x][0]=(f[x][0]+f[x][1])%MOD;///x是白点,儿子要是有黑点,砍了x的父边就是没黑点,所以没黑点(f[x][0])的情况要加上有黑点的情况(f[x][1])
64     else f[x][1]=f[x][0];///x点是黑点,那不砍父边的情况(f[x][1])只有让x的儿子都不黑,砍父边的情况(f[x][0])也是x的儿子都不黑,因为x自己黑嘛,儿子再黑就连到一起了
65     //printf("[out %d,flag=%d,re=%I64d,a[x]=%d]\n",x,flag,re,a[x]);
66 }
67 
68 
69 ll farm() {
70     if(n==1)return 1;
71     mz(u);
72     dfs(0);
73     return f[0][1];
74 }
75 
76 int main() {
77     int i;
78     int x;
79     RD(n);
80     memset(head,-1,sizeof(head));
81     en=0;
82     REP(i,n-1) {
83         scanf("%d",&x);
84         add(i+1,x);
85         add(x,i+1);
86     }
87     for(i=0; i<n; i++)
88         scanf("%d",&a[i]);
89     printf("%I64d",farm());
90     return 0;
91 }
View Code

 

转载于:https://www.cnblogs.com/yuiffy/p/3940042.html

error #6649: This symbol has multiple INTENT statement/attribute declarations which is not allowed. [Y] Compilation Aborted (code 1) Compilation Aborted (code 1) error #6406: Conflicting attributes or multiple declaration of name. [Y] error #6406: Conflicting attributes or multiple declaration of name. [Y] error #6649: This symbol has multiple INTENT statement/attribute declarations which is not allowed. [Y]修改以下程序! 凝结尾迹形成模拟程序 ! 包含Schmidt-Appleman判据计算和Gyarmathy液滴生长模型 ! 添加了完整的控制方程组求解 !====================================================================== module SchmidtAppleman implicit none !!定义常量 real, parameter :: a_water = 7.63 !!a,b为经验常数,水面取7.63和241.9 冰面取9.5和265.5 real, parameter :: b_water = 241.9 real, parameter :: Q = 35e3 !!燃料的燃烧热 real, parameter :: b = 265.5 real, parameter :: EIH2O = 1.25 !!水蒸气排放指数(假设为1.25) real, parameter :: Cp = 1.004 !!空气的定压比热容 real, parameter :: epsilon = 0.6222 !!水干空气的分子量之比 real, parameter :: a = 9.5 real, parameter :: E0 = 6.11 !!温度为0时,E0的数值,hPa(马格努斯经验公式) real, parameter :: E0_water = 6.112 real, parameter :: EIT = 0.6 !!喷气发动机的平均推进效率 contains subroutine check_contrail_formation() real :: RHc !!温度T下的凝结尾形成的阈值湿度,计算可得 real :: RHi !!冰面饱和水汽,计算可得 real :: RHw !!相对于水面的饱和水汽压,kPa real :: G !!计算可得 real :: Tc !!计算可得 real :: RH real :: RH_0 = 80.0 !! 添加初始值80%相对湿度 real :: T=-50 !!大气环境温度 ℃ real :: P=26500 !!环境的大气压强 real :: e_T !!T温度下的饱和水汽压 real :: eTc !!温度T下的饱和水汽压,由马格努斯公式计算可得 integer :: have RH = RH_0 / 100.0 RHw = E0_water*10**((a_water*T)/(b_water + T)) RHi = RHw*(6.061*exp(18.102*T/(249.52+T))/6.1162/exp(22.577*T/(237.78+T))) G = (EIH2O * Cp * P) / (epsilon * Q * (1 - EIT)) Tc = -46.46 + 9.43*log(G-0.053) + 0.72*log(G-0.053)*log(G-0.053) eTc = E0*10**(a*T/(b+T)) RHc = (G*(T-Tc) + eTc)/eTc print *, "G, Tc, RHw =", G, Tc, RHw If (RHc <= RH) then If (RHc < 1) then If (RHi < 1) then have = 1 write (*,*) "emergence,能产生凝结尾迹" else have = 0 write (*,*) "will not emergence,不能产生凝结尾迹" end if end if end if end subroutine check_contrail_formation end module SchmidtAppleman !====================================================================== ! 凝结尾迹综合模拟模块 ! 包含气相方程(有限差分法)和液相方程(迎风格式) !====================================================================== module contrail_simulation use, intrinsic :: iso_fortran_env, only: dp => real64 implicit none !------------------ 物理常数定义 ------------------! real(dp), parameter :: pi = 3.1415926535897932_dp ! 圆周率 real(dp), parameter :: R_gas = 287.0_dp ! 空气气体常数 [J/(kg·K)] real(dp), parameter :: R_v = 461.5_dp ! 水蒸气气体常数 [J/(kg·K)] real(dp), parameter :: k_B = 1.380649e-23_dp ! 玻尔兹曼常数 [J/K] real(dp), parameter :: NA = 6.02214076e23_dp ! 阿伏伽德罗常数 [mol⁻¹] real(dp), parameter :: Cp_air = 1004.0_dp ! 空气定压比热容 [J/(kg·K)] real(dp), parameter :: lambda_v = 0.026_dp ! 蒸汽导热系数 [W/(m·K)] real(dp), parameter :: h_tg = 2.5e6_dp ! 汽化潜热 [J/kg] real(dp), parameter :: Pr = 0.71_dp ! 普朗特数(空气) real(dp), parameter :: lambda = 1.0_dp ! 分子自由程比 real(dp), parameter :: sigma = 0.072_dp ! 表面张力 [N/m] real(dp), parameter :: rho_l = 1000.0_dp ! 液态水密度 [kg/m³] real(dp), parameter :: mol_diameter = 3.0e-10_dp ! 水分子直径 [m] real(dp), parameter :: mu = 1.8e-5_dp ! 空气动力粘度 [Pa·s] real(dp), parameter :: k_thermal = 0.024_dp ! 空气热导率 [W/(m·K)] !------------------ 环境参数 ------------------! real(dp) :: Te = 223.15_dp ! 环境温度 [K] real(dp) :: P_env = 26500.0_dp ! 环境总压 [Pa] real(dp) :: Ue = 2000.0_dp ! 对流速度 [m/s] (修正为更合理的值) !------------------ 射流参数 ------------------! real(dp) :: T_jet = 1000.0_dp ! 射流温度 [K] real(dp) :: U_jet = 2000.0_dp ! 射流速度 [m/s] real(dp) :: P_jet = 50000.0_dp ! 射流压力 [Pa] !------------------ 计算网格参数 ------------------! integer, parameter :: nx = 100 ! 轴向网格数 integer, parameter :: nr = 50 ! 径向网格数 real(dp), parameter :: L = 50.0_dp ! 轴向计算域长度 [m] real(dp), parameter :: R_max = 10.0_dp ! 最大径向距离 [m] real(dp), parameter :: dx = L / nx ! 轴向步长 [m] real(dp), parameter :: dr = R_max / nr ! 径向步长 [m] real(dp), parameter :: dt = 1.0e-7_dp ! 时间步长 [s] (减小以适应高速射流) !------------------ 场变量定义 ------------------! ! 气相场变量 real(dp), dimension(nx, nr) :: rho_g = 0.3_dp ! 气相密度 [kg/m³] real(dp), dimension(nx, nr) :: u = 2000.0_dp ! 轴向速度 [m/s] (初始化为射流速度) real(dp), dimension(nx, nr) :: v = 0.0_dp ! 径向速度 [m/s] (轴对称假设) real(dp), dimension(nx, nr) :: p = 50000.0_dp ! 压力 [Pa] (初始化为射流压力) real(dp), dimension(nx, nr) :: T_g = 1000.0_dp ! 温度 [K] (初始化为射流温度) real(dp), dimension(nx, nr) :: h_g = 0.0_dp ! 焓 [J/kg] ! 液相场变量 real(dp), dimension(nx, nr) :: N_d = 1.0e-20_dp ! 液滴数目密度 [1/m³] real(dp), dimension(nx, nr) :: Y_d = 0.0_dp ! 液滴质量分数 [kg/kg] real(dp), dimension(nx, nr) :: r_d = 1.0e-9_dp ! 液滴半径 [m] ! 源项 real(dp), dimension(nx, nr) :: S_m = 0.0_dp ! 质量源项 real(dp), dimension(nx, nr) :: S_u = 0.0_dp ! 动量源项 real(dp), dimension(nx, nr) :: S_v = 0.0_dp ! 径向动量源项 real(dp), dimension(nx, nr) :: S_h = 0.0_dp ! 能量源项 real(dp) :: time = 0.0_dp ! 模拟时间 [s] integer :: step_count = 0 ! 时间步计数 contains !==================== 主计算程序 ====================! subroutine solve() real(dp) :: drdt, J ! 半径增长率 [m/s], 成核率 real(dp) :: mv_nucleation, mv_growth ! 成核生长质量源项 real(dp) :: vapor_consumption ! 蒸气消耗量 [kg/m³] real(dp) :: T_C, S, Kn, R_c, Delta_T ! 中间计算变量 real(dp) :: mean_free_path ! 分子平均自由程 [m] real(dp) :: x_pos, r_pos ! 位置坐标 real(dp) :: conv_term, diff_term ! 对流项和扩散项 real(dp) :: rho_air ! 空气密度 [kg/m³] real(dp) :: dudx, dudy, dvdx, dvdy ! 速度梯度 real(dp) :: dTdx, dTdy, dpdx, dpdy ! 梯度 real(dp) :: div_u, laplace_u ! 散度和拉普拉斯算子 integer :: i, j, step real(dp) :: area_factor ! 轴对称面积因子 !=============== 初始化阶段 ================! rho_air = P_env / (R_gas * Te) ! 计算空气密度 h_g = Cp_air * T_g ! 初始化焓值 ! 设置射流初始条件 (整个射流区域) do j = 1, nr do i = 1, nx rho_g(i,j) = P_jet / (R_gas * T_jet) ! 使用射流条件计算密度 u(i,j) = U_jet v(i,j) = 0.0_dp p(i,j) = P_jet T_g(i,j) = T_jet ! 仅在中心区域初始化液滴 if (i <= nx/5 .and. j <= nr/2) then N_d(i,j) = 1.0e15_dp Y_d(i,j) = 0.01_dp r_d(i,j) = 1.0e-6_dp else N_d(i,j) = 1.0e-20_dp Y_d(i,j) = 0.0_dp r_d(i,j) = 1.0e-9_dp end if end do end do ! 打开输出文件 open(unit=10, file="contrail_evolution.csv", status="replace") write(10, "(A)") "Step,Time[s],x[m],r[m],N[1/m3],Y[kg/kg],R[m],T[K],P[Pa],S" ! 打开云图输出文件 open(unit=20, file="temperature_field.csv", status="replace") open(unit=30, file="pressure_field.csv", status="replace") write(20, "(A)") "x[m],r[m],T[K]" write(30, "(A)") "x[m],r[m],P[Pa]" !=============== 时间推进循环 ================! do step = 1, 200000 ! 增加时间步数以适应更小的时间步长 step_count = step_count + 1 ! 0. 重置源项 S_m = 0.0_dp S_u = 0.0_dp S_v = 0.0_dp S_h = 0.0_dp ! 1. 液相求解 (迎风格式) do j = 2, nr-1 do i = 2, nx-1 ! 计算位置坐标 x_pos = (i - 0.5_dp) * dx r_pos = (j - 0.5_dp) * dr area_factor = 1.0_dp / r_pos ! 轴对称因子 ! 液滴数密度方程 (7) conv_term = 0.0_dp ! 轴向对流 (迎风) if (u(i,j) > 0.0_dp) then conv_term = conv_term + u(i,j) * (N_d(i,j) - N_d(i-1,j)) / dx else conv_term = conv_term + u(i,j) * (N_d(i+1,j) - N_d(i,j)) / dx end if ! 径向对流 (迎风) if (v(i,j) > 0.0_dp) then conv_term = conv_term + v(i,j) * (N_d(i,j) - N_d(i,j-1)) / dr else conv_term = conv_term + v(i,j) * (N_d(i,j+1) - N_d(i,j)) / dr end if ! 添加轴对称项 conv_term = conv_term + v(i,j) * N_d(i,j) * area_factor ! 更新液滴数密度 N_d(i,j) = N_d(i,j) - dt * (conv_term - J) ! 液滴质量分数方程 (8) conv_term = 0.0_dp ! 轴向对流 (迎风) if (u(i,j) > 0.0_dp) then conv_term = conv_term + u(i,j) * (Y_d(i,j) - Y_d(i-1,j)) / dx else conv_term = conv_term + u(i,j) * (Y_d(i+1,j) - Y_d(i,j)) / dx end if ! 径向对流 (迎风) if (v(i,j) > 0.0_dp) then conv_term = conv_term + v(i,j) * (Y_d(i,j) - Y_d(i,j-1)) / dr else conv_term = conv_term + v(i,j) * (Y_d(i,j+1) - Y_d(i,j)) / dr end if ! 添加轴对称项 conv_term = conv_term + v(i,j) * Y_d(i,j) * area_factor ! 更新质量分数 Y_d(i,j) = Y_d(i,j) - dt * (conv_term - mv_growth) ! 更新液滴半径 (9) if (N_d(i,j) > 1e-10_dp) then r_d(i,j) = (3.0_dp * Y_d(i,j) * rho_g(i,j) / (4.0_dp * pi * rho_l * N_d(i,j)))**(1.0_dp/3.0_dp) end if ! 计算成核和生长源项 mean_free_path = k_B * T_g(i,j) / (sqrt(2.0_dp)*pi*(mol_diameter**2)*p(i,j)) T_C = T_g(i,j) - 273.15_dp S = p(i,j) / (611.2_dp * exp(22.46_dp*T_C/(T_C + 272.62_dp))) R_c = (2.0_dp*sigma) / (rho_l * R_gas * T_g(i,j) * log(S)) J = 1e-3_dp * (p(i,j)/(R_gas*T_g(i,j)))**2 / rho_l * & sqrt(2.0_dp*sigma/(pi*(18e-3_dp/NA)**3)) * & exp(-4.0_dp*pi*R_c**2*sigma/(3.0_dp*k_B*T_g(i,j))) Delta_T = (T_g(i,j)**2 * R_v / h_tg) * log(S) Kn = mean_free_path / r_d(i,j) drdt = (lambda_v*Delta_T*(1.0_dp - R_c/r_d(i,j))) / & (p(i,j)*h_tg*r_d(i,j)*(1.0_dp + & (2.0_dp/(1.5_dp*Pr))*(r_d(i,j)/(1.0_dp+lambda))*Kn)) mv_nucleation = J * (4.0_dp/3.0_dp)*pi*R_c**3*rho_l mv_growth = N_d(i,j) * 4.0_dp*pi*r_d(i,j)**2 * drdt * rho_l ! 设置源项 (4-6) S_m(i,j) = -mv_growth S_u(i,j) = u(i,j) * S_m(i,j) S_v(i,j) = v(i,j) * S_m(i,j) S_h(i,j) = h_tg * S_m(i,j) end do end do ! 2. 气相求解 (有限差分法) do j = 2, nr-1 do i = 2, nx-1 ! 计算位置坐标 x_pos = (i - 0.5_dp) * dx r_pos = (j - 0.5_dp) * dr area_factor = 1.0_dp / r_pos ! 轴对称因子 ! 连续性方程 (1) dudx = (u(i+1,j) - u(i-1,j)) / (2.0_dp*dx) dvdy = (v(i,j+1) - v(i,j-1)) / (2.0_dp*dr) div_u = dudx + dvdy + v(i,j)*area_factor rho_g(i,j) = rho_g(i,j) - dt * (rho_g(i,j)*div_u + u(i,j)*(rho_g(i+1,j)-rho_g(i-1,j))/(2.0_dp*dx) & + v(i,j)*(rho_g(i,j+1)-rho_g(i,j-1))/(2.0_dp*dr)) + dt*S_m(i,j) ! x方向动量方程 (2) dudx = (u(i+1,j) - u(i-1,j)) / (2.0_dp*dx) dudy = (u(i,j+1) - u(i,j-1)) / (2.0_dp*dr) dpdx = (p(i+1,j) - p(i-1,j)) / (2.0_dp*dx) laplace_u = (u(i+1,j) - 2.0_dp*u(i,j) + u(i-1,j))/(dx*dx) + & (u(i,j+1) - 2.0_dp*u(i,j) + u(i,j-1))/(dr*dr) u(i,j) = u(i,j) - dt * (u(i,j)*dudx + v(i,j)*dudy + (1.0_dp/rho_g(i,j))*dpdx) & + dt * (mu/rho_g(i,j)) * (laplace_u + (1.0_dp/3.0_dp)*dudx) & + dt * S_u(i,j)/rho_g(i,j) ! r方向动量方程 (2) - 轴对称形式 dvdx = (v(i+1,j) - v(i-1,j)) / (2.0_dp*dx) dvdy = (v(i,j+1) - v(i,j-1)) / (2.0_dp*dr) dpdy = (p(i,j+1) - p(i,j-1)) / (2.0_dp*dr) laplace_v = (v(i+1,j) - 2.0_dp*v(i,j) + v(i-1,j))/(dx*dx) + & (v(i,j+1) - 2.0_dp*v(i,j) + v(i,j-1))/(dr*dr) v(i,j) = v(i,j) - dt * (u(i,j)*dvdx + v(i,j)*dvdy + (1.0_dp/rho_g(i,j))*dpdy) & + dt * (mu/rho_g(i,j)) * (laplace_v - v(i,j)/(r_pos*r_pos) + (1.0_dp/3.0_dp)*dvdy) & + dt * S_v(i,j)/rho_g(i,j) ! 能量方程 (3) dTdx = (T_g(i+1,j) - T_g(i-1,j)) / (2.0_dp*dx) dTdy = (T_g(i,j+1) - T_g(i,j-1)) / (2.0_dp*dr) h_g(i,j) = h_g(i,j) - dt * (u(i,j)*dTdx + v(i,j)*dTdy) * Cp_air & + dt * (k_thermal/(rho_g(i,j)*Cp_air)) * & ((T_g(i+1,j)-2.0_dp*T_g(i,j)+T_g(i-1,j))/(dx*dx) + & ((T_g(i,j+1)-2.0_dp*T_g(i,j)+T_g(i,j-1))/(dr*dr)) & + dt * S_h(i,j)/rho_g(i,j) T_g(i,j) = h_g(i,j) / Cp_air ! 状态方程 p(i,j) = rho_g(i,j) * R_gas * T_g(i,j) end do end do ! 3. 边界条件 ! 入口边界 (i=1) - 射流条件 i = 1 do j = 1, nr rho_g(i,j) = P_jet / (R_gas * T_jet) u(i,j) = U_jet v(i,j) = 0.0_dp p(i,j) = P_jet T_g(i,j) = T_jet if (j <= nr/2) then N_d(i,j) = 1.0e15_dp Y_d(i,j) = 0.01_dp r_d(i,j) = 1.0e-6_dp else N_d(i,j) = 1.0e-20_dp Y_d(i,j) = 0.0_dp r_d(i,j) = 1.0e-9_dp end if end do ! 出口边界 (i=nx) - 零梯度 i = nx do j = 1, nr rho_g(i,j) = rho_g(i-1,j) u(i,j) = u(i-1,j) v(i,j) = v(i-1,j) p(i,j) = P_env ! 出口压力设为环境压力 T_g(i,j) = T_g(i-1,j) N_d(i,j) = N_d(i-1,j) Y_d(i,j) = Y_d(i-1,j) r_d(i,j) = r_d(i-1,j) end do ! 对称轴 (j=1) j = 1 do i = 1, nx v(i,j) = 0.0_dp ! 对称条件 rho_g(i,j) = rho_g(i,j+1) u(i,j) = u(i,j+1) p(i,j) = p(i,j+1) T_g(i,j) = T_g(i,j+1) N_d(i,j) = N_d(i,j+1) Y_d(i,j) = Y_d(i,j+1) r_d(i,j) = r_d(i,j+1) end do ! 外边界 (j=nr) - 环境条件 j = nr do i = 1, nx rho_g(i,j) = rho_air u(i,j) = Ue v(i,j) = 0.0_dp p(i,j) = P_env T_g(i,j) = Te N_d(i,j) = 1.0e-20_dp Y_d(i,j) = 0.0_dp r_d(i,j) = 1.0e-9_dp end do ! 4. 输出结果 if (mod(step,100) == 0) then ! 输出中心线数据 i = nx/2 j = 1 x_pos = (i - 0.5_dp) * dx r_pos = (j - 0.5_dp) * dr T_C = T_g(i,j) - 273.15_dp S = p(i,j) / (611.2_dp * exp(22.46_dp*T_C/(T_C + 272.62_dp))) write(10,'(I0,9(",",ES12.5))') step, time, x_pos, r_pos, & N_d(i,j), Y_d(i,j), r_d(i,j), T_g(i,j), p(i,j), S print *, "步数:",step," 时间:",time," 温度:",T_g(i,j)," 压力:",p(i,j) ! 输出云图数据 (每500步) if (mod(step,500) == 0) then rewind(20) rewind(30) write(20, "(A)") "x[m],r[m],T[K]" write(30, "(A)") "x[m],r[m],P[Pa]" do j = 1, nr do i = 1, nx x_pos = (i - 0.5_dp) * dx r_pos = (j - 0.5_dp) * dr write(20, '(3(ES12.5,","))') x_pos, r_pos, T_g(i,j) write(30, '(3(ES12.5,","))') x_pos, r_pos, p(i,j) end do end do end if end if time = time + dt ! 终止条件:超过0.05秒 (高速射流需要更少的时间) if (time > 0.5_dp) exit end do !=============== 结束阶段 ================! close(10) close(20) close(30) ! 输出最终云图 open(unit=40, file="final_temperature.csv", status="replace") open(unit=50, file="final_pressure.csv", status="replace") write(40, "(A)") "x[m],r[m],T[K]" write(50, "(A)") "x[m],r[m],P[Pa]" do j = 1, nr do i = 1, nx x_pos = (i - 0.5_dp) * dx r_pos = (j - 0.5_dp) * dr write(40, '(3(ES12.5,","))') x_pos, r_pos, T_g(i,j) write(50, '(3(ES12.5,","))') x_pos, r_pos, p(i,j) end do end do close(40) close(50) print *, "=== 模拟完成 ===" print "(A,ES10.3)", "最终时间 [s]: ", time print "(A,ES10.3)", "最大液滴半径 [m]: ", maxval(r_d) print "(A,ES10.3)", "最高温度 [K]: ", maxval(T_g) print "(A,ES10.3)", "最低压力 [Pa]: ", minval(p) print "(A,ES10.3)", "最大速度 [m/s]: ", maxval(u) end subroutine solve end module contrail_simulation !======================== 主程序 ======================== program main use SchmidtAppleman use contrail_simulation implicit none ! 执行凝结尾迹形成判据 call check_contrail_formation() ! 执行综合模拟 call solve() end program main
最新发布
06-03
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值