TAO提供的initial reference

本文通过一个示例程序展示了如何使用 TAO (The ACE ORB) 获取初始化对象引用和服务列表。该程序运行于 Windows XP 系统下,通过 C++ 代码调用 TAO 相关 API 列出了所有可用的初始化服务。

应用程序需要可移值意味着它可以获取初始化对象引用。TAO提供了哪些初始化应用呢?我们可以通过程序来获取。

版本OCI TAO 1.6a

Windows XP

源代码:

#include "ace/OS_main.h"
#include "ace/Log_Msg.h"
#include "tao/tao/ObjectIdListC.h"
#include "tao/ORB.h"
#include "ace/streams.h"
int ACE_TMAIN(int argc, ACE_TCHAR *argv[])
{
	ACE_DEBUG((LM_DEBUG,ACE_TEXT("(%t) test start here./n")));
	try
	{
		CORBA::ORB_var orb = CORBA::ORB_init(argc,argv,"test");
	    CORBA::ORB_ObjectIdList_var lis = orb->list_initial_services();

		CORBA::ULong n = lis->length();
		for(CORBA::ULong index = 0UL; index < n; index ++)
		{
			 const char* name  = lis[index];
			 cout << index << "  " << name << endl;
		}
		orb->shutdown();
		orb->destroy();
	}	
	catch (CORBA::Exception& ex)
	{
		ex._tao_print_exception("main");
	}
	return 0;
}
 输出结果

0  NameService
1  TradingService
2  ImplRepoService
3  RootPOA
4  POACurrent
5  InterfaceRepository
6  ORBPolicyManager
7  PolicyCurrent
8  IORManipulation
9  IORTable
10  DynAnyFactory
11  TypeCodeFactory
12  CompressionManager
13  Monitor

这是我的代码,你看下电流密度为什么会从3900变到3400 原理是什么,什么公式 ! For modeling oxygen transport in Catalyst layer, details can be found in the following literature. ! Mesoscopic modeling impacts of liquid water saturation, and platinum distribution on gas transport resistances in a PEMFC catalyst layer ! Mesoscopic modeling of transport resistances in a polymer-electrolyte fuel-cell catalyst layer: Analysis of hydrogen limiting currents ! Please cite this reference: YT Mu, SR Yang, P He, WQ Tao, Electrochimica Acta 388, 138659 ! Please cite this reference: YT Mu, AZ Weber, ZL Gu, WQ Tao, Applied Energy 255, 113895-113908 ! Developed by Yutong Mu from Xi'an Jiaotong University, Email: yutongmu@mail.xjtu.edu.cn module start integer i,j,k,ij,mm,iter,nprocs,istart,iend,ntot,ncase parameter (nxmax=100,nprocsx=10,nx=nxmax/nprocsx) parameter (ny=100,nz=100) parameter (q = 19) ! D3Q19 model double precision,parameter::XL=2000.E-9,YL=300.E-9,ZL=300.E-9 integer,dimension(0:nx+1,0:ny+1,0:nz+1)::ls logical,dimension(0:nx+1,0:ny+1,0:nz+1)::wallg,walle,wallPt double precision,dimension(0:nx+1,0:ny+1,0:nz+1)::Por_o2,Naf_o2,diff double precision,dimension(0:nx+1,0:ny+1,0:nz+1)::chmcal,chPt ! D3Q19 discrete velocities integer,dimension(1:q)::fcx=(/1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 1,-1, 1,-1, 0, 0, 0, 0, 0/) integer,dimension(1:q)::fcy=(/0, 0, 1,-1, 0, 0, 1, 1,-1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0/) integer,dimension(1:q)::fcz=(/0, 0, 0, 0, 1,-1, 0, 0, 0, 0, 1, 1,-1,-1, 1, 1,-1,-1, 0/) double precision,dimension(0:q-1,0:nx+1,0:ny+1,0:nz+1)::fp,ffp,fn,ffn double precision,dimension(0:q-1)::w ! weight coefficients double precision::tau_p,tau_n ! BGK relaxation times double precision::chgin,Dg,Delec0,Delec,Dliq double precision::GNa_pn,Gup1_pn,Gup2_pn,Nup1_pn,Nup2_pn double precision::DX,DT,scale,sita,fx,EtaD,j0,kk1,kk2,kk3,PtO,Enerst double precision::farady,gascons,concref,excoeff,currentref,Etain,pressure,tmeaue double precision::one,two,react,RPt,RPt0,SumRate,TotRate,lamda_eq,scale_o2 double precision::RGDL,act,henry double precision::CurDensity,Rcl_delta,Rcl_old double precision::porosity double precision::Rcl,xiter0 double precision::RF,SumRF,SumRF0,ECSA double precision::kgdl double precision::totPot,totPor double precision::sumPot,sumPor double precision::K_GDL,K_CL,miu_liq,distance,theta_CL,theta_GDL double precision::Pc_CL,Pc_GDL,sigma,por_GDL,sat_GDL CHARACTER*20::name,fname,name1 integer idx,myid integer mxp,mxm integer nallgrp,ierr end module PROGRAM MAIN use start use mpi integer :: ramp_duration, iter_load double precision :: Etain_ramp_step, Etain_initial, Etain_final integer istatus(MPI_STATUS_SIZE) integer, parameter :: iter_total = 200000 integer, parameter :: output_freq = 200 integer :: transient_output_unit double precision :: sumNaf, totNaf ! =================== 1. initialize MPI =================== call MPI_INIT(ierr) call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr) nallgrp=MPI_COMM_WORLD idx=myid mxp=mod(idx+1,nprocsx) mxm=mod(idx+nprocsx-1,nprocsx) istart=1; iend=nx if(idx.eq.0) istart=3 if(idx.eq.nprocsx-1) iend=nx-2 !====================== initialize structure ===================== dx=3.0d-9 fx=0.5d0 pressure=101325.d0+110000.d0 tmeaue=273.15d0+80.d0 farady=96487.d0 gascons=8.314d0 concref=3.39d0 excoeff=0.5d0 currentref=60.d0 Dg =1.47e-6 ! *** MODIFICATION 1: Break the circular dependency *** ! First, set a fixed, stable value for tau_p. It MUST be > 0.5. tau_p = 0.6d0 ! Initialize D3Q19 weights call init_D3Q19 ! Second, calculate dt and other parameters based on the fixed tau_p. j0=1.d0/7.d0 sita=(1.d0-j0)*(tau_p-0.5d0)/3.d0 scale=Dg/sita dt=dx**2.d0/scale EtaD=(1.d0-j0)/3.d0 lamda_eq = 8.7d0 Delec0=(1.14698d0/10000000000.d0)*(lamda_eq**0.708d0) henry=101325.d0/(4.408d0-0.09712d0*lamda_eq) kp_pn=henry/gascons/tmeaue scale_o2=1.d0/5.d0 !scale_o2=10.0d0/5.d0 Delec=scale_o2*Delec0 ! Third, calculate tau_n based on the fixed tau_p and other parameters. tau_n=0.5d0+Delec/Dg*(tau_p-0.5d0) if(myid.eq.0) write(*,*) 'tau_p (fixed), tau_n (calculated), kp_pn:' if(myid.eq.0) write(*,*) tau_p,tau_n,kp_pn GNa_pn=-1.d0/kp_pn-1.d0 Gup1_pn=1.d0/kp_pn-1.d0 Gup2_pn=-2.d0 Nup1_pn=-1.d0/kp_pn+1.d0 Nup2_pn=-2.d0/kp_pn if(myid.eq.0) write(*,*) GNa_pn,Gup1_pn,Gup2_pn,Nup1_pn,Nup2_pn ! Now that dt is a valid non-zero number, these calculations will be correct. kk1=currentref*(dt/dx)/concref kk2=excoeff*farady/gascons/tmeaue kk3=3000.d0/gascons/tmeaue Rpt=7.3d0; Rpt0=Rpt*dx/dt ! Ramp loading parameters Etain_initial = 1.0d0 Etain_final = 0.6d0 iter_load = 20000 ramp_duration = 1 Etain_ramp_step = (Etain_final - Etain_initial) / dble(ramp_duration) Etain = Etain_initial Enerst=1.2d0-Etain PtO=1.d0/(1.d0+exp(22.4d0*(0.818d0-Enerst))) react=(1.d0-PtO)*kk1*exp(kk2*Etain-kk3*PtO) one=react*(EtaD*Rpt0+1.d0)/4.d0/farady two=one+EtaD if(myid.eq.0) then write(*,*) "Initial Etain: ", Etain write(*,*) "Initial PtO, one, two: ", PtO, one, two transient_output_unit = 20 open(transient_output_unit, file='transient_response.dat') write(transient_output_unit, '(A)') '# Time(s), CurrentDensity(A/m^2), Rcl(s/m), Avg_O2_Pore, Avg_O2_Nafion' endif if(myid.eq.0) write(*,*) "Initial PtO = ", PtO, kk2*Etain, kk3*PtO, one, two RGDL=20.d0 kgdl=dt/dx/RGDL if(myid.eq.0) write(*,*) kk1,kk2,kk3,RGDL call MYTSTRUCTURE chgtot =0.48952d0 chgin =chgtot do k=1,NZ do j=1,NY do i=1,NX if(.not.wallg(i,j,k)) then Por_o2(i,j,k)=chgin Naf_o2(i,j,k)=chgin/kp_pn call initialize_distributions(fp(:,i,j,k), Por_o2(i,j,k), w) call initialize_distributions(ffp(:,i,j,k), Por_o2(i,j,k), w) call initialize_distributions(fn(:,i,j,k), Naf_o2(i,j,k), w) call initialize_distributions(ffn(:,i,j,k), Naf_o2(i,j,k), w) else if(.not.walle(i,j,k)) then Por_o2(i,j,k)=chgin Naf_o2(i,j,k)=chgin/kp_pn call initialize_distributions(fp(:,i,j,k), Por_o2(i,j,k), w) call initialize_distributions(ffp(:,i,j,k), Por_o2(i,j,k), w) call initialize_distributions(fn(:,i,j,k), Naf_o2(i,j,k), w) call initialize_distributions(ffn(:,i,j,k), Naf_o2(i,j,k), w) else Por_o2(i,j,k)=chgin Naf_o2(i,j,k)=chgin/kp_pn call initialize_distributions(fp(:,i,j,k), Por_o2(i,j,k), w) call initialize_distributions(ffp(:,i,j,k), Por_o2(i,j,k), w) call initialize_distributions(fn(:,i,j,k), Naf_o2(i,j,k), w) call initialize_distributions(ffn(:,i,j,k), Naf_o2(i,j,k), w) endif enddo enddo enddo Rcl_old=0.0 do ITER=1,iter_total if (ITER >= iter_load .and. ITER < iter_load + ramp_duration) then Etain = Etain + Etain_ramp_step Enerst=1.2d0-Etain PtO=1.d0/(1.d0+exp(22.4d0*(0.818d0-Enerst))) react=(1.d0-PtO)*kk1*exp(kk2*Etain-kk3*PtO) one=react*(EtaD*Rpt0+1.d0)/4.d0/farady two=one+EtaD else if (ITER >= iter_load + ramp_duration) then if (Etain /= Etain_final) then Etain = Etain_final Enerst=1.2d0-Etain PtO=1.d0/(1.d0+exp(22.4d0*(0.818d0-Enerst))) react=(1.d0-PtO)*kk1*exp(kk2*Etain-kk3*PtO) one=react*(EtaD*Rpt0+1.d0)/4.d0/farady two=one+EtaD endif endif CALL COLLISION_BGK CALL STREAM_D3Q19 CALL LBM_BOUNDARY_D3Q19 CALL MACRO_BOUNDARY IF(mod(ITER,output_freq).eq.0) then CurDensity = TotRate*4.d0*farady/real(ny*nz) Rcl=chgin*4.d0*farady/CurDensity Rcl_delta=abs(Rcl-Rcl_old)/(Rcl+1.E-20) Rcl_old=Rcl call Cal_RoughFactor if(myid.eq.0) write(*,*) "RoughFactor= (m Pt 2 * m MEA -2)",SumRF chgin=0.85d0*chgin+0.15d0*(chgtot*Rcl/(Rcl+RGDL)) totPot=0.d0 totPor=0.d0 totNaf=0.d0 do k=1,nz do j=1,ny do i=istart,iend if((.not.wallg(i,j,k))) then totPot=totPot+Por_o2(i,j,k) totPor=totPor+1.d0 else if (.not.walle(i,j,k)) then totNaf=totNaf+Naf_o2(i,j,k) endif enddo enddo enddo call MPI_REDUCE(totPot,sumPot,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr) call MPI_REDUCE(totPor,sumPor,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr) call MPI_REDUCE(totNaf,sumNaf,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr) if(myid.eq.0) then write(*,*) "iter,CurDensity,Rcl,Rcl*sumRF :" write(*,*) iter,CurDensity,Rcl,Rcl*sumRF write(*,*) "chgin,Rcl_delta :" write(*,*) chgin,Rcl_delta write(*,*) "henry,tau_n,Delec :" write(*,*) henry,tau_n,Delec write(*,*) "sumPot/sumPor, lamda_eq " write(*,*) sumPot/sumPor,lamda_eq write(transient_output_unit, '(E12.5, 4(A, E12.5))') & ITER*dt, ',', CurDensity, ',', Rcl, ',', & sumPot/(sumPor+1.e-20), ',', sumNaf/((nxmax-4)*ny*nz*porosity + 1.e-20) endif if(iter.gt.20000 .and. Rcl_delta.le.1.d-5) then goto 110 endif endif enddo 110 continue if (myid.eq.0) then close(transient_output_unit) endif call output call MPI_FINALIZE(ierr) CONTAINS !************************************************************************ ! *** MODIFICATION 2: This subroutine now only initializes the weights. *** subroutine init_D3Q19 use start ! D3Q19 weight coefficients w(0) = 1.0d0/3.0d0 w(1:6) = 1.0d0/18.0d0 w(7:18) = 1.0d0/36.0d0 if(myid.eq.0) then write(*,*) "D3Q19 BGK model weights initialized" endif end subroutine subroutine initialize_distributions(f, rho, weight) use start double precision, intent(out) :: f(0:q-1) double precision, intent(in) :: rho, weight(0:q-1) integer :: m do m = 0, q-1 f(m) = rho * weight(m) enddo end subroutine subroutine equilibrium_D3Q19(feq, rho, u, v, w_vel, weight) use start double precision, intent(out) :: feq(0:q-1) double precision, intent(in) :: rho, u, v, w_vel, weight(0:q-1) integer :: m double precision :: cu, c_sq, u_sq u_sq = u*u + v*v + w_vel*w_vel do m = 0, q-1 cu = fcx(m+1)*u + fcy(m+1)*v + fcz(m+1)*w_vel feq(m) = rho * weight(m) * (1.d0 + 3.d0*cu + 4.5d0*cu*cu - 1.5d0*u_sq) enddo end subroutine subroutine MYTSTRUCTURE use start use mpi integer istatus(MPI_STATUS_SIZE) integer lentot,ii double precision,dimension(0:nxmax+1,0:ny+1,0:nz+1)::difftot integer,dimension(0:nxmax+1,0:ny+1,0:nz+1)::lstot if(myid.eq.0) then lstot=0 porosity=0 open(10, file="LSTOT.dat") do k=1,nz do j=1,ny do i=1,nxmax read(10,*) lstot(i,j,k) if(i.ge.3 .and. i.le.nxmax-2 .and. lstot(i,j,k).eq.0) porosity=porosity+1 enddo enddo enddo close(10) lstot(0:2,:,:)=0 lstot(nxmax-1:nxmax+1,:,:)=0 lstot(:,0,:)=lstot(:,ny,:) lstot(:,ny+1,:)=lstot(:,1,:) lstot(:,:,0)=lstot(:,:,nz) lstot(:,:,nz+1)=lstot(:,:,1) porosity=real(porosity)/real(ny*nz*(nxmax-4)) write(*,*) "Porosity = ", porosity difftot(:,:,:)=Dg open(10, file="DIFFTOT.dat") do k=1,nz do j=1,ny do i=1,nxmax read(10,*) difftot(i,j,k) enddo enddo enddo close(10) difftot(0:2,:,:)=Dg difftot(nxmax-1:nxmax+1,:,:)=Dg endif call MPI_BCAST (porosity,1,MPI_DOUBLE_PRECISION,0,nallgrp,ierr) lentot=(nxmax+2)*(ny+2)*(nz+2) call MPI_BCAST(lstot,lentot,MPI_integer,0,nallgrp,ierr) do k=0,nz+1 do j=0,ny+1 do i=0,nx+1 ii=nx*idx+i ls(i,j,k)=lstot(ii,j,k) enddo enddo enddo lentot=(nxmax+2)*(ny+2)*(nz+2) call MPI_BCAST(difftot,lentot,MPI_DOUBLE_PRECISION,0,nallgrp,ierr) do k=0,nz+1 do j=0,ny+1 do i=0,nx+1 ii=nx*idx+i diff(i,j,k)=difftot(ii,j,k) enddo enddo enddo wallg=.false. walle=.false. wallPt=.false. do k=0,nz+1 do j=0,ny+1 do i=0,nx+1 if(ls(i,j,k).eq.0) then diff(i,j,k)=0.5d0+diff(i,j,k)/Dg*(tau_p-0.5d0) endif if(ls(i,j,k).ne.0) wallg(i,j,k)=.true. if(ls(i,j,k).ne.1) walle(i,j,k)=.true. if(ls(i,j,k).ne.2) wallPt(i,j,k)=.true. enddo enddo enddo return end subroutine Subroutine Cal_RoughFactor use start use mpi integer istatus(MPI_STATUS_SIZE) double precision,dimension(0:nx+1,0:ny+1,0:nz+1)::ttemp !===================== Obtain the Real Roughness Factor ========================= ttemp=0.d0 do k=1,nz do j=1,ny do i=1,nx mm=0 if(Naf_o2(i,j,k).gt.1.d-6) then if(ls(i+1,j,k).eq.2) mm=mm+1 if(ls(i-1,j,k).eq.2) mm=mm+1 if(ls(i,j+1,k).eq.2) mm=mm+1 if(ls(i,j-1,k).eq.2) mm=mm+1 if(ls(i,j,k+1).eq.2) mm=mm+1 if(ls(i,j,k-1).eq.2) mm=mm+1 endif ttemp(i,j,k)=mm enddo enddo enddo RF=sum(ttemp(:,:,:)) call MPI_REDUCE(RF,SumRF,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr) call MPI_BCAST (SumRF,1,MPI_DOUBLE_PRECISION,0,nallgrp,ierr) if(myid.eq.0) then open(11,file="UpdateRoughFactor_1.txt") SumRF=SumRF/real(ny*nz) write(11,*) "RoughFactor= (m Pt 2 * m MEA -2)",SumRF ECSA=SumRF/real(nxmax-4)/21.45d0/1.0e6/dx write(11,*) "Specific surface area over volumePt (m 2 * g -1)=",ECSA close(11) endif return End Subroutine subroutine COLLISION_BGK use start use mpi integer istatus(MPI_STATUS_SIZE) double precision,dimension(0:q-1)::feq_p, feq_n double precision::rho_p, rho_n integer :: m_i, j_i, k_i do k_i=1,nz do j_i=1,ny do i=1,nx if(.not.wallg(i,j_i,k_i)) then rho_p = 0.d0 do m_i=0,q-1 rho_p = rho_p + fp(m_i,i,j_i,k_i) enddo Por_o2(i,j_i,k_i) = rho_p call equilibrium_D3Q19(feq_p, rho_p, 0.d0, 0.d0, 0.d0, w) do m_i=0,q-1 ffp(m_i,i,j_i,k_i) = fp(m_i,i,j_i,k_i) - (fp(m_i,i,j_i,k_i) - feq_p(m_i)) / tau_p enddo else if(.not.walle(i,j_i,k_i)) then rho_n = 0.d0 do m_i=0,q-1 rho_n = rho_n + fn(m_i,i,j_i,k_i) enddo Naf_o2(i,j_i,k_i) = rho_n call equilibrium_D3Q19(feq_n, rho_n, 0.d0, 0.d0, 0.d0, w) do m_i=0,q-1 ffn(m_i,i,j_i,k_i) = fn(m_i,i,j_i,k_i) - (fn(m_i,i,j_i,k_i) - feq_n(m_i)) / tau_n enddo endif enddo enddo enddo return end subroutine subroutine STREAM_D3Q19 use start use mpi integer istatus(MPI_STATUS_SIZE) double precision, dimension(0:q-1,0:ny+1,0:nz+1):: tmp1,tmp2,tmp1n,tmp2n integer left,right,tag1,tag2,tag3,tag4 integer ilen, ii, jj, kk, m_dir integer :: opp_dir integer :: i_loc, j_loc, k_loc ! MPI communication for boundary regions if(nprocsx>1) then ilen=(ny+2)*(nz+2)*q right=mxp left=mxm if(idx==0) left=MPI_PROC_NULL if(idx==nprocsx-1) right=MPI_PROC_NULL tmp1(:,:,:)=ffp(:,nx,:,:) tmp2=0.d0 tag1=10 call MPI_SENDRECV(tmp1(0,0,0),ilen,MPI_DOUBLE_PRECISION,right& &,tag1,tmp2(0,0,0),ilen,MPI_DOUBLE_PRECISION,left,tag1,nallgrp,istatus,ierr) ffp(:,0,:,:)=tmp2(:,:,:) tmp1(:,:,:)=ffp(:,1,:,:) tmp2=0.d0 tag2=20 call MPI_SENDRECV(tmp1(0,0,0),ilen,MPI_DOUBLE_PRECISION,left& &,tag2,tmp2(0,0,0),ilen,MPI_DOUBLE_PRECISION,right,tag2,nallgrp,istatus,ierr) ffp(:,nx+1,:,:)=tmp2(:,:,:) tmp1n(:,:,:)=ffn(:,nx,:,:) tmp2n=0.d0 tag3=30 call MPI_SENDRECV(tmp1n(0,0,0),ilen,MPI_DOUBLE_PRECISION,right& &,tag3,tmp2n(0,0,0),ilen,MPI_DOUBLE_PRECISION,left,tag3,nallgrp,istatus,ierr) ffn(:,0,:,:)=tmp2n(:,:,:) tmp1n(:,:,:)=ffn(:,1,:,:) tmp2n=0.d0 tag4=40 call MPI_SENDRECV(tmp1n(0,0,0),ilen,MPI_DOUBLE_PRECISION,left& &,tag4,tmp2n(0,0,0),ilen,MPI_DOUBLE_PRECISION,right,tag4,nallgrp,istatus,ierr) ffn(:,nx+1,:,:)=tmp2n(:,:,:) endif ! Periodic boundaries in y and z directions ffp(:,:,0,:)=ffp(:,:,ny,:) ffp(:,:,ny+1,:)=ffp(:,:,1,:) ffp(:,:,:,0)=ffp(:,:,:,nz) ffp(:,:,:,nz+1)=ffp(:,:,:,1) ffn(:,:,0,:)=ffn(:,:,ny,:) ffn(:,:,ny+1,:)=ffn(:,:,1,:) ffn(:,:,:,0)=ffn(:,:,:,nz) ffn(:,:,:,nz+1)=ffn(:,:,:,1) ! Perform streaming do k_loc=1,nz do j_loc=1,ny do i_loc=1,nx do m_dir=0,q-1 ii = i_loc - fcx(m_dir+1) jj = j_loc - fcy(m_dir+1) kk = k_loc - fcz(m_dir+1) ! C-style periodic boundary conditions ii = mod(ii - 1 + nx, nx) + 1 jj = mod(jj - 1 + ny, ny) + 1 kk = mod(kk - 1 + nz, nz) + 1 fp(m_dir,i_loc,j_loc,k_loc) = ffp(m_dir,ii,jj,kk) fn(m_dir,i_loc,j_loc,k_loc) = ffn(m_dir,ii,jj,kk) enddo enddo enddo enddo call Nafion_Pore_D3Q19 call Nafion_Pt_D3Q19 ! Bounce-back for solid boundaries do k_loc=1,nz do j_loc=1,ny do i_loc=1,nx if(wallg(i_loc,j_loc,k_loc)) then do m_dir=0,q-1 opp_dir = find_opposite(m_dir) ffp(m_dir,i_loc,j_loc,k_loc) = fp(opp_dir,i_loc,j_loc,k_loc) enddo endif if(walle(i_loc,j_loc,k_loc)) then do m_dir=0,q-1 opp_dir = find_opposite(m_dir) ffn(m_dir,i_loc,j_loc,k_loc) = fn(opp_dir,i_loc,j_loc,k_loc) enddo endif enddo enddo enddo return end subroutine integer function find_opposite(dir) integer, intent(in) :: dir integer :: opp select case(dir) case(0) opp = 0 case(1) opp = 2 case(2) opp = 1 case(3) opp = 4 case(4) opp = 3 case(5) opp = 6 case(6) opp = 5 case(7) opp = 10 case(8) opp = 9 case(9) opp = 8 case(10) opp = 7 case(11) opp = 14 case(12) opp = 13 case(13) opp = 12 case(14) opp = 11 case(15) opp = 18 case(16) opp = 17 case(17) opp = 16 case(18) opp = 15 case default opp = 0 end select find_opposite = opp end function subroutine LBM_BOUNDARY_D3Q19 use start use mpi integer istatus(MPI_STATUS_SIZE) integer :: m_dir, opp_dir integer :: j_loc, k_loc if(idx==0) then i=1 ! Inlet boundary - constant concentration do m_dir=0,q-1 if(fcx(m_dir+1) > 0) then opp_dir = find_opposite(m_dir) do k_loc=1,nz do j_loc=1,ny fp(m_dir,i,j_loc,k_loc) = EtaD*chgin - fp(opp_dir,i,j_loc,k_loc) enddo enddo endif enddo ! Nafion phase at inlet do m_dir=0,q-1 do k_loc=1,nz do j_loc=1,ny fn(m_dir,i,j_loc,k_loc) = ffn(m_dir,i,j_loc,k_loc) enddo enddo enddo endif if(idx==nprocsx-1) then i=nx ! Outlet boundary - zero gradient do m_dir=0,q-1 if(fcx(m_dir+1) < 0) then opp_dir = find_opposite(m_dir) do k_loc=1,nz do j_loc=1,ny fp(m_dir,i,j_loc,k_loc) = fp(opp_dir,i,j_loc,k_loc) enddo enddo endif enddo do m_dir=0,q-1 do k_loc=1,nz do j_loc=1,ny fn(m_dir,i,j_loc,k_loc) = ffn(m_dir,i,j_loc,k_loc) enddo enddo enddo endif return end subroutine Subroutine output use start use mpi double precision::act0,lamda0 integer istatus(MPI_STATUS_SIZE) integer :: IU IU=myid+101 write(name,'(I3.3)') IU fname='Rate_1_'//TRIM(ADJUSTL(name))//'.DAT' open (301, file=fname) do i=1,nx do j=1,ny do k=1,nz write (301,*) 4.d0*farady*chmcal(i,j,k)/dt enddo enddo enddo close(301) fname='Cone_1_'//TRIM(ADJUSTL(name))//'.DAT' open (301, file=fname) do i=1,nx do j=1,ny do k=1,nz write (301,*) Naf_o2(i,j,k) enddo enddo enddo close(301) fname='Conc_1_'//TRIM(ADJUSTL(name))//'.DAT' open (301, file=fname) do i=1,nx do j=1,ny do k=1,nz write (301,*) Por_o2(i,j,k) enddo enddo enddo close(301) fname='ChPt_1_'//TRIM(ADJUSTL(name))//'.DAT' open (301, file=fname) do i=1,nx do j=1,ny do k=1,nz write (301,*) chPt(i,j,k) enddo enddo enddo close(301) return End Subroutine subroutine MACRO_BOUNDARY use start use mpi integer istatus(MPI_STATUS_SIZE) integer :: m_dir do k=1,nz do j=1,ny do i=1,nx if(.not.wallg(i,j,k)) then Por_o2(i,j,k)=0.d0 do m_dir=0,q-1 Por_o2(i,j,k)=Por_o2(i,j,k)+fp(m_dir,i,j,k) enddo Por_o2(i,j,k)=max(Por_o2(i,j,k),0.d0) Naf_o2(i,j,k)=0.d0 else if(.not.walle(i,j,k)) then Naf_o2(i,j,k)=0.d0 do m_dir=0,q-1 Naf_o2(i,j,k)=Naf_o2(i,j,k)+fn(m_dir,i,j,k) enddo Naf_o2(i,j,k)=max(Naf_o2(i,j,k),0.d0) Por_o2(i,j,k)=0.d0 else Por_o2(i,j,k)=0.d0 Naf_o2(i,j,k)=0.d0 endif enddo enddo enddo return end subroutine subroutine Nafion_Pore_D3Q19 use start use mpi integer istatus(MPI_STATUS_SIZE) integer :: m_dir, ii, jj, kk, opp do k=1,NZ do j=1,NY do i=1,NX if(.not.walle(i,j,k)) then do m_dir=0,q-1 ii = i + fcx(m_dir+1) jj = j + fcy(m_dir+1) kk = k + fcz(m_dir+1) if(ii >= 1 .and. ii <= nx .and. jj >= 1 .and. jj <= ny .and. kk >= 1 .and. kk <= nz) then if(.not.wallg(ii,jj,kk)) then ! Interface between Nafion and pore opp = find_opposite(m_dir) ffn(m_dir,ii,jj,kk) = Nup1_pn/GNa_pn*ffn(opp,i,j,k) + Nup2_pn/GNa_pn*ffp(m_dir,ii,jj,kk) endif endif enddo endif if(.not.wallg(i,j,k)) then do m_dir=0,q-1 ii = i + fcx(m_dir+1) jj = j + fcy(m_dir+1) kk = k + fcz(m_dir+1) if(ii >= 1 .and. ii <= nx .and. jj >= 1 .and. jj <= ny .and. kk >= 1 .and. kk <= nz) then if(.not.walle(ii,jj,kk)) then ! Interface between pore and Nafion opp = find_opposite(m_dir) ffp(m_dir,ii,jj,kk) = Gup1_pn/GNa_pn*ffp(opp,i,j,k) + Gup2_pn/GNa_pn*ffn(m_dir,ii,jj,kk) endif endif enddo endif enddo enddo enddo return end subroutine subroutine Nafion_Pt_D3Q19 use start use mpi implicit none integer istatus(MPI_STATUS_SIZE) double precision::che,cmol integer::nchPt, opp integer::ii, jj, kk, m_dir SumRate=0.d0; chmcal=0.d0; chPt=0.d0 do k=1,nz do j=1,ny do i=1,nx nchPt=0 if(.not.walle(i,j,k)) then do m_dir=0,q-1 ii = i + fcx(m_dir+1) jj = j + fcy(m_dir+1) kk = k + fcz(m_dir+1) if(ii >= 1 .and. ii <= nx .and. jj >= 1 .and. jj <= ny .and. kk >= 1 .and. kk <= nz) then if(.not.wallPt(ii,jj,kk)) then ! Interface between Nafion and Pt opp = find_opposite(m_dir) che=2.d0*ffn(opp,i,j,k)/two cmol=(2.d0*ffn(opp,i,j,k)*Rpt0+che)/(EtaD*Rpt0+1.d0) ffn(m_dir,ii,jj,kk)=EtaD*cmol-ffn(opp,i,j,k) SumRate=SumRate+(ffn(opp,i,j,k)-ffn(m_dir,ii,jj,kk)) chmcal(i,j,k)=chmcal(i,j,k)+(ffn(opp,i,j,k)-ffn(m_dir,ii,jj,kk)) chPt(i,j,k)=chPt(i,j,k)+che nchPt=nchPt+1 endif endif enddo endif chPt(i,j,k)=chPt(i,j,k)/(float(nchPt)+1.E-20) enddo enddo enddo call MPI_REDUCE(SumRate,TotRate,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,nallgrp,ierr) call MPI_BCAST (TotRate,1,MPI_DOUBLE_PRECISION,0,nallgrp,ierr) TotRate=TotRate*dx/dt return end subroutine end PROGRAM MAIN
最新发布
11-01
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值