CF 461B Appleman and Tree 树形DP

本文介绍了一种利用树形动态规划方法解决特定树分割问题的方法。该问题要求在给定一棵带颜色标记的树后,计算出将树分割成若干部分,使得每部分恰好包含一个黑色节点的所有可能方案的数量。

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

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
每一个节点涂有黑色或者白色,1代表黑色,0代表白色
若在树上去掉k条边,就把树分成k+1部分,每一个部分也是一棵树,若每一部分都有且只有一个节点是黑色,
则这是一个合理的操作。
求合理操作的方案数%(1e9+7)

如果把黑色看成节点的值为1,白色看成节点的值为0
一棵树的值=树上所有节点的值之和
则这道题转化为:
把一棵树分成若干个部分,每一部分的值都为1的方案数。

树形DP

dp[i][0] : 以i为根的子树,i所在部分的值为0的方案数%mod
dp[i][1] : 以i为根的子树,i所在部分的值为1的方案数%mod

以root=0进行DFS
则输出:dp[0][1]


代码:



 1 #include<cstdio>
 2 #include<cstring>
 3 
 4 using namespace std;
 5 
 6 const int maxn=1e5+5;
 7 const int mod=1e9+7;
 8 #define ll long long
 9 
10 struct Edge
11 {
12     int to,next;
13 };
14 Edge edge[maxn<<1];
15 int head[maxn];
16 int tot=0;
17 ll dp[maxn][2];
18 int cost[maxn];
19 
20 void init()
21 {
22     memset(head,-1,sizeof head);
23     memset(dp,0,sizeof dp);
24 }
25 
26 void addedge(int u,int v)
27 {
28     edge[tot].to=v;
29     edge[tot].next=head[u];
30     head[u]=tot++;
31 }
32 
33 void solve(int );
34 void dfs(int ,int );
35 
36 int main()
37 {
38     init();
39     int n;
40     scanf("%d",&n);
41     for(int i=1;i<n;i++)
42     {
43         int p;
44         scanf("%d",&p);
45         addedge(i,p);
46         addedge(p,i);
47     }
48     for(int i=0;i<n;i++)
49     {
50         scanf("%d",&cost[i]);
51     }
52     solve(n);
53     return 0;
54 }
55 
56 void solve(int n)
57 {
58     dfs(0,-1);
59     printf("%d\n",(int)dp[0][1]);
60     return ;
61 }
62 
63 void dfs(int u,int pre)
64 {
65     if(cost[u])
66     {
67         dp[u][1]=1;
68         dp[u][0]=0;
69     }
70     else
71     {
72         dp[u][0]=1;
73         dp[u][1]=0;
74     }
75     for(int i=head[u];~i;i=edge[i].next)
76     {
77         int v=edge[i].to;
78         if(v==pre)
79             continue;
80         dfs(v,u);
81         if(cost[u])
82         {
83             dp[u][1]*=(dp[v][0]+dp[v][1]);
84             dp[u][1]%=mod;
85         }
86         else
87         {
88             dp[u][1]=dp[u][1]*(dp[v][0]+dp[v][1])+dp[u][0]*dp[v][1];
89             dp[u][1]%=mod;
90             dp[u][0]*=(dp[v][0]+dp[v][1]);
91             dp[u][0]%=mod;
92         }
93     }
94 }
View Code

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

























转载于:https://www.cnblogs.com/-maybe/p/4767145.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、付费专栏及课程。

余额充值