Main loop
subroutine main_loop() !! Main loop real(dp) :: error_val call timer%set(9,description='Main loop start') if(sim%rollback) then ! Rollback in progress: state already set, skip state advance sim%rollback = .false. h5%steps_to_revert = 3 else if(sim%rejected) then ! Retrieve previous state channels%StVar = channels%StVarOld strands%StVar = strands%StVarOld solids%StVar = solids%StVarOld mesh2Ds%StVar = mesh2Ds%StVarOld else ! Advance state variables channels%StVarOld6 = channels%StVarOld5 channels%StVarOld5 = channels%StVarOld4 channels%StVarOld4 = channels%StVarOld3 channels%StVarOld3 = channels%StVarOld2 channels%StVarOld2 = channels%StVarOld channels%StVarOld = channels%StVar strands%StVarOld6 = strands%StVarOld5 strands%StVarOld5 = strands%StVarOld4 strands%StVarOld4 = strands%StVarOld3 strands%StVarOld3 = strands%StVarOld2 strands%StVarOld2 = strands%StVarOld strands%StVarOld = strands%StVar solids%StVarOld6 = solids%StVarOld5 solids%StVarOld5 = solids%StVarOld4 solids%StVarOld4 = solids%StVarOld3 solids%StVarOld3 = solids%StVarOld2 solids%StVarOld2 = solids%StVarOld solids%StVarOld = solids%StVar mesh2Ds%StVarOld6 = mesh2Ds%StVarOld5 mesh2Ds%StVarOld5 = mesh2Ds%StVarOld4 mesh2Ds%StVarOld4 = mesh2Ds%StVarOld3 mesh2Ds%StVarOld3 = mesh2Ds%StVarOld2 mesh2Ds%StVarOld2 = mesh2Ds%StVarOld mesh2Ds%StVarOld = mesh2Ds%StVar endif ! Qext_Check_step=0.0_dp !$omp parallel do private(i) do i=1, size(channels) call friction_source_term(channels(i)) call channels_FF_flux_port_comm(channels(i)) enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1, size(junctions) call junction_resolution_from_and_to_ports(junctions(i)) enddo !$omp end parallel do if (handle_numerical_error()) return do i=1, size(circulators) call circulator_resolution_from_and_to_ports(circulators(i)) enddo if (handle_numerical_error()) return do i=1, size(boundaries) call boundary_resolution_from_and_to_ports(boundaries(i)) enddo if (handle_numerical_error()) return !$omp parallel do private(i) do i=1, size(channels) call He_Riemann_solver_channel(channels(i),sim) call flux_from_FF_flux_ports_self(channels(i)) enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1, size(strands) call scenario_update(strands(i),sim) call strands_SS_flux_port_comm(strands(i)) enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1, size(solids) call solid_scenario_update(solids(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(SSfluxLinks) call SSfluxLink_resolution_from_and_to_ports(SSfluxLinks(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(strands) call heat_diffusion_strand(strands(i)) call flux_from_SS_flux_ports_self(strands(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(solids) call heat_diffusion_solid(solids(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(mesh2Ds) call heat_diffusion_mesh2D(mesh2Ds(i),sim) call mesh2Ds_FS_port_comm(mesh2Ds(i)) call mesh2Ds_SS_src_port_comm(mesh2Ds(i)) enddo !$omp end parallel do ! Calculation of minimal time for wave propagation through the node ! and updating explicit step call sim%update_exp_dt( min( & minVal(channels%wave_time), & minVal(junctions%wave_time), & minVal(circulators%wave_time), & minVal(boundaries%wave_time))) call timer%set(10,9,'Fluxes + Riemann') call timer%set(11) if(sim%explicit) then ! explicit fluxes / implicit source terms ------------------------------ !print*,'Explicit scheme at physical time',sim%t,'(in s)' !$omp parallel do private(i) do i=1,size(channels) call explicit_scheme_for_fluxes_channel(channels(i),sim%dt) call friction_source_term(channels(i)) call source_term_definition_channel(channels(i),sim) enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1,size(strands) call explicit_scheme_for_fluxes_strand(strands(i),sim%dt) call source_term_definition_strand(strands(i),sim) !,Qext_Check_step) ! Attention in case of parallelization --> Qext_Check enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1,size(solids) call explicit_scheme_for_fluxes_solid(solids(i),sim%dt) call source_term_definition_solid(solids(i),sim) enddo !$omp end parallel do else ! full implicit ------------------------------------------------------- !$omp parallel do private(i) do i=1,size(channels) call full_physics_definition_channel(channels(i),sim) enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1,size(strands) call full_physics_definition_strand(strands(i),sim) !,Qext_Check_step) enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1,size(solids) call full_physics_definition_solid(solids(i),sim) enddo !$omp end parallel do !$omp parallel do private(i) do i=1,size(junctions) call links_junction_update(junctions(i),sim%dt) enddo !$omp end parallel do do i=1,size(circulators) call links_circulator_update(circulators(i),sim%dt) enddo !$omp parallel do private(i) do i=1,size(SSfluxLinks) call links_SSfluxLink_update(SSfluxLinks(i),sim%dt) enddo !$omp end parallel do endif ! -------------------------------------------------------------------- ! Qext_Check_case_step=0.0_dp !$omp parallel do private(i) do i=1,size(mesh2Ds) ! SurfReg=0.0_dp ! do nb_e=1,mesh2Ds(i)%M2D_Prop%nb_elements ! if(trim(mesh2Ds(i)%M2D_Prop%elem(nb_e)%PhysE)=='SS_Outer') then ! SurfReg=SurfReg+mesh2Ds(i)%M2D_Prop%elem(nb_e)%surface ! endif ! enddo call full_physics_definition_mesh2D(mesh2Ds(i),sim) !,Qext_Check_case_loc) ! Qext_Check_case_step = Qext_Check_case_step + Qext_Check_case_loc * SurfReg * 1.062_dp enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(FSlinks) call FSlink_src_reinitialization(FSlinks(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(FSlinks) call FSlink_resolution_from_and_to_ports(FSlinks(i)) call links_FSlink_update(FSlinks(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(channels) call source_from_FS_ports_self_channel(channels(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(strands) call source_from_FS_ports_self_strand(strands(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(solids) call source_from_FS_ports_self_solid(solids(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(mesh2Ds) call source_from_FS_ports_self_mesh2D(mesh2Ds(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(SSsrcLinks) call SSsrcLink_src_reinitialization(SSsrcLinks(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(SSsrcLinks) call SSsrcLink_resolution_from_and_to_ports(SSsrcLinks(i)) call links_SSsrcLink_update(SSsrcLinks(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(solids) call source_from_SS_src_ports_self_solid(solids(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(mesh2Ds) call source_from_SS_src_ports_self_mesh2D(mesh2Ds(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(FFsrcLinks) call FFsrcLink_resolution_from_and_to_ports(FFsrcLinks(i)) call links_FFsrcLink_update(FFsrcLinks(i),sim%dt) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(channels) call source_from_FF_src_ports_self_channel(channels(i),sim%dt) enddo !$omp end parallel do if (handle_numerical_error()) return ! ---------------------- L I N E A R S O L V E R ------------------------- call Pardiso_fact_solve(krn,sim%explicit,timer) if (handle_numerical_error()) return ! -------------------------------------------------------------------------- call timer%set(12) !$omp parallel do private(i) do i=1,size(channels) call from_sol_to_prim_channel(channels(i)) enddo !$omp end parallel do if (handle_numerical_error()) return !$omp parallel do private(i) do i=1,size(strands) call from_sol_to_temp_strand(strands(i),sim) enddo !$omp end parallel do !$omp parallel do private(i) do i=1,size(solids) call from_sol_to_temp_solid(solids(i),sim) enddo !$omp end parallel do !$omp parallel do private(i) do i=1,size(mesh2Ds) call from_sol_to_temp_mesh2D(mesh2Ds(i),sim) enddo !$omp end parallel do !$omp parallel do private(i) do i=1,size(channels) call channels_to_FF_src_port_comm(channels(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1, size(FFsrcLinks) call FFsrcLink_Prelax_resolution_from_and_to_ports(FFsrcLinks(i)) enddo !$omp end parallel do !$omp parallel do private(i) do i=1,size(channels) call FF_src_port_to_channels_comm(channels(i)) call last_tasks_for_channels(channels(i),sim) enddo !$omp end parallel do if (handle_numerical_error()) return if(sim%explicit) then ! pressure criteria for switching to implicit error_val = maxVal(sqrt(channels%err) / max(1.0e-6_dp, sqrt(channels%err_den))) else; error_val = & ! compute local truncation error (lte) - energy calculation as ! criteria for implicit step sqrt(sum(channels%err)) / max(sqrt(sum(channels%err_den)), sim_delta) + & sqrt(sum(strands%err)) / max(sqrt(sum(strands%err_den)), sim_delta) + & sqrt(sum(solids%err)) / max(sqrt(sum(solids%err_den)), sim_delta) + & sqrt(sum(mesh2Ds%err)) / max(sqrt(sum(mesh2Ds%err_den)), sim_delta) endif call timer%set(13,12,'Solution distribution') call sim%step_management(error_val) if(.not.sim%rejected) call h5%write(sim%t,sim%dt,sim%t_final) if(h5%write_results2D .and. .not.sim%rejected) then if(sim%t >= next2DWriteTime) then step2D = step2D + 1 do i=1,size(mesh2Ds) call writing_HDF5_2D_casing(mesh2Ds(i),sim%t,step2D) enddo next2DWriteTime = next2DWriteTime + h5%time_btw_2D_writes endif endif ! if(.not.sim%rejected) then ! mDotHOut = channels(1)%flxHe%Cons(Con_Ene,15) * channels(1)%HeProp%Area ! mDotHin = channels(1)%flxHe%Cons(Con_Ene,5)*channels(1)%HeProp%Area ! checkMDotHDt = checkMDotHDt + (mDotHOut-mDotHin)*sim%dt ! write(132,'(4(e13.6,1x))') sim%t, FullEnergy_MC, FullEnergy_2D, -checkMDotHDt ! Qext_Check = Qext_Check + Qext_Check_step ! write(133,'(2(e13.6,1x))') sim%t, Qext_Check ! Qext_Check_case = Qext_Check_case + Qext_Check_case_step ! write(133,'(2(e13.6,1x))') sim%t, Qext_Check_case ! do i=1,size(solids) ! do ii=1,solids(i)%MC_Prop%NbCells ! if(abs(solids(i)%StVarOld%MCtemp(ii))>1.0e-10_dp) then ! Q_MC = Q_MC + solids(i)%MC_Prop%ro_M * solids(i)%MC_Prop%volLoc(ii) * & ! (CpSS_Integral(solids(i)%StVar%MCtemp(ii)) - CpSS_Integral(solids(i)%StVarOld%MCtemp(ii))) ! if(i==1) then ! Q_MC1 = Q_MC1 + solids(i)%MC_Prop%ro_M * solids(i)%MC_Prop%volLoc(ii) * & ! (CpSS_Integral(solids(i)%StVar%MCtemp(ii)) - CpSS_Integral(solids(i)%StVarOld%MCtemp(ii))) ! else if(i==2) then ! Q_MC2 = Q_MC2 + solids(i)%MC_Prop%ro_M * solids(i)%MC_Prop%volLoc(ii) * & ! (CpSS_Integral(solids(i)%StVar%MCtemp(ii)) - CpSS_Integral(solids(i)%StVarOld%MCtemp(ii))) ! endif ! endif ! enddo ! enddo ! do i=1,size(mesh2Ds) ! do ii=1,mesh2Ds(i)%M2D_Prop%nb_elements ! if(abs(mesh2Ds(i)%StVarOld%temp(ii))>1.0e-10_dp) then ! Q_2D = Q_2D + 7900.0_dp * mesh2Ds(i)%M2D_Prop%elem(ii)%surface * 0.5_dp * & ! 0.5_dp for 0.5m of contact with MC ! (CpSS_Integral(mesh2Ds(i)%StVar%temp(ii)) - CpSS_Integral(mesh2Ds(i)%StVarOld%temp(ii))) ! endif ! enddo ! enddo ! write(134,'(9(e13.6,1x))') sim%t, Q_MC, Q_MC1, Q_MC2, Q_2D ! endif end subroutine main_loop