obsidian_backup/多体+耦合求解器/多体+水动HD+系泊MD.md
2026-02-12 23:22:38 +08:00

399 lines
19 KiB
Markdown
Raw Permalink Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

# 配置文件
- 配置文件新增MD和MD输入文件路径加在气动后 done
- 新增p_fast.comp_mooring done
# init HD MD初始化
对照FAST_subs.f90 774行附近
```fortran
IF ( p_FAST%CompHydro == Module_HD ) THEN
```
对照FAST_subs.f90 919行附近
```fortran
! initialize CompMooring modules
```
**补充代码** coupled_solver_sub.rs
```rust
pub fn fast_init_with_yaml(){}
    if p_fast.comp_hydro == 1 {
    to be done
    }
   
    if p_fast.comp_mooring == 1{
    to be done
    }
   
}
```
模型文件路径data\HEROWIND_models\NREL_5MW_OC4Semi\
配置文件Herowind_NREL5MW_ppl_hdmd_calc_config.yaml中气动控制HD MD绝对路径需按需修改
路径使用示例代码:
配置文件参数
```rust
calc_config.hydrofile.hydro_input_file
calc_config.mooringfile.mooring_input_file
```
代码
```rust
// 1 获取运行程序exe所在的路径
    let exe_path: PathBuf = env::current_exe()
        .expect("Unable to get current executable path");
    let exe_dir = exe_path
        .parent()
        .expect("Executable must have a parent directory");
       
    // 2 路径拼接
let dll_path = &calc_config.controller_input.dll_path;
let dll_full_path = exe_dir.join(dll_path);
let dll_full_path_str: &str = dll_full_path.to_str().expect("proxy_path is not valid UTF8"); 
     
```
## InitModuleMappings
mapping code
```fortran
if (p_FAST%CompElast == Module_ED) then
NumBl = SIZE(ED%y%BladeRootMotion,1)
PlatformMotion => ED%y%PlatformPtMesh
PlatformLoads => ED%Input(1)%PlatformPtMesh
SubstructureMotion2HD => ED%y%PlatformPtMesh
SubstructureMotion => ED%y%PlatformPtMesh
SubstructureLoads => ED%Input(1)%PlatformPtMesh
IF ( SrvD%Input(1)%PtfmMotionMesh%Committed ) THEN
CALL MeshMapCreate( PlatformMotion, SrvD%Input(1)%PtfmMotionMesh, MeshMapData%ED_P_2_SrvD_P_P, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':ED_P_2_SrvD_P_P' )
ENDIF
IF ( p_FAST%CompHydro == Module_HD ) THEN ! HydroDyn-{ElastoDyn or SubDyn}
! Regardless of the offshore configuration, ED platform motions will be mapped to the PRPMesh of HD
! we're just going to assume PlatformLoads and PlatformMotion are committed
CALL MeshMapCreate( PlatformMotion, HD%Input(1)%PRPMesh, MeshMapData%ED_P_2_HD_PRP_P, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':ED_P_2_HD_PRP_P' )
IF ( HD%y%WAMITMesh%Committed ) THEN ! meshes for floating
! HydroDyn WAMIT point mesh to/from ElastoDyn or SD point mesh
CALL MeshMapCreate( HD%y%WAMITMesh, SubstructureLoads, MeshMapData%HD_W_P_2_SubStructure, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':HD_W_P_2_SubStructure' )
CALL MeshMapCreate( SubstructureMotion2HD, HD%Input(1)%WAMITMesh, MeshMapData%SubStructure_2_HD_W_P, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':SubStructure_2_HD_W_P' )
END IF
IF ( HD%Input(1)%Morison%Mesh%Committed ) THEN
CALL MeshMapCreate( HD%y%Morison%Mesh, SubstructureLoads, MeshMapData%HD_M_P_2_SubStructure, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':HD_M_P_2_SubStructure' )
CALL MeshMapCreate( SubstructureMotion2HD, HD%Input(1)%Morison%Mesh, MeshMapData%SubStructure_2_HD_M_P, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':SubStructure_2_HD_M_P' )
END IF
IF ( p_FAST%CompMooring == Module_MD ) THEN
!-------------------------
! SubDyn/ElastoDyn <-> MoorDyn
!-------------------------
! MoorDyn point mesh to/from SubDyn or ElastoDyn point mesh
CALL MeshMapCreate( MD%y%CoupledLoads(1), SubstructureLoads, MeshMapData%Mooring_2_Structure, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':Mooring_2_Structure' )
CALL MeshMapCreate( SubstructureMotion, MD%Input(1)%CoupledKinematics(1), MeshMapData%Structure_2_Mooring, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':Structure_2_Mooring' )
END IF
IF ( p_FAST%SolveOption == Solve_SimplifiedOpt1 ) THEN
CALL AllocAry( MeshMapData%Jacobian_Opt1, SizeJac_ED_HD, SizeJac_ED_HD, 'Jacobian for Ptfm-HD coupling', ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
IF ( ALLOCATED( MeshMapData%Jacobian_Opt1 ) ) THEN
CALL AllocAry( MeshMapData%Jacobian_pivot, SIZE(MeshMapData%Jacobian_Opt1), 'Pivot array for Jacobian LU decomposition', ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
END IF
CALL ResetRemapFlags(p_FAST, ED, SED, BD, AD, ExtLd, HD, SD, ExtPtfm, SrvD, MAPp, FEAM, MD, Orca, IceF, IceD )
CALL MeshCopy ( ED%Input(1)%HubPtLoad, MeshMapData%u_ED_HubPtLoad, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':u_ED_HubPtLoad' )
CALL MeshCopy ( SubStructureLoads, MeshMapData%SubstructureLoads_Tmp, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':SubstructureLoads_Tmp' )
CALL MeshCopy ( SubStructureLoads, MeshMapData%SubstructureLoads_Tmp2, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':SubstructureLoads_Tmp2' )
CALL MeshCopy ( PlatformLoads, MeshMapData%PlatformLoads_Tmp, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':PlatformLoads_Tmp' )
CALL MeshCopy ( PlatformLoads, MeshMapData%PlatformLoads_Tmp2, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':PlatformLoads_Tmp2' )
IF ( p_FAST%CompHydro == Module_HD ) THEN
!TODO: GJH Is this needed, I created it as a place holder, 5/11/2020
!CALL MeshCopy ( HD%Input(1)%PRPMesh, MeshMapData%u_HD_PRP_Mesh, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
! CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':u_HD_PRP_Mesh' )
CALL MeshCopy ( HD%Input(1)%WAMITMesh, MeshMapData%u_HD_W_Mesh, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':u_HD_W_Mesh' )
CALL MeshCopy ( HD%Input(1)%Morison%Mesh, MeshMapData%u_HD_M_Mesh, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':u_HD_M_Mesh' )
END IF
END SUBROUTINE InitModuleMappings
```
### SUBROUTINE ResetRemapFlags
```fortran
IF ( p_FAST%CompHydro == Module_HD ) THEN
HD%Input(1)%PRPMesh%RemapFlag = .FALSE.
IF (HD%Input(1)%WAMITMesh%Committed) THEN
HD%Input(1)%WAMITMesh%RemapFlag = .FALSE.
HD%y%WAMITMesh%RemapFlag = .FALSE.
END IF
IF (HD%Input(1)%Morison%Mesh%Committed) THEN
HD%Input(1)%Morison%Mesh%RemapFlag = .FALSE.
HD%y%Morison%Mesh%RemapFlag = .FALSE.
END IF
END IF
IF ( p_FAST%CompMooring == Module_MD ) THEN
MD%Input(1)%CoupledKinematics(1)%RemapFlag = .FALSE.
MD%y%CoupledLoads(1)%RemapFlag = .FALSE.
```
# solution t0
```fortran
if ( P_FAST%CompSeaSt == Module_SeaSt .and. y_FAST%WriteThisStep) then
! note: SeaState has no inputs and only calculates WriteOutputs, so we don't need to call CalcOutput unless we are writing to the file
call SeaSt_CalcOutput( t_initial, SeaSt%u, SeaSt%p, SeaSt%x(1), SeaSt%xd(1), SeaSt%z(1), SeaSt%OtherSt(1), SeaSt%y, SeaSt%m, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
end if
```
## CALL CalcOutputs_And_SolveForInputs()
### solver_option_2基本没变化
新增 transfer_ServD_to_SD_MD
传递u_MD%DeltaL    =  y_SrvD%CableDeltaL  和  u_MD%DeltaLdot =  y_SrvD%CableDeltaLdot
新增CALL Transfer_Structure_to_Opt1Inputs
### Transfer_Structure_to_Opt1Inputs
```fortran
PlatformMotion => y_ED%PlatformPtMesh
IF (p_FAST%CompSub == Module_SD) THEN
SubstructureMotion => SD%y%y3Mesh
SubstructureMotion2HD => SD%y%y2Mesh
ELSE
SubstructureMotion => PlatformMotion
SubstructureMotion2HD => PlatformMotion
ENDIF
IF ( p_FAST%CompHydro == Module_HD ) THEN
CALL Transfer_Point_to_Point( PlatformMotion, u_HD%PRPMesh, MeshMapData%ED_P_2_HD_PRP_P, ErrStat2, ErrMsg2 )
CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat, ErrMsg, RoutineName//' (u_HD%PRPMesh)' )
! if we don't have a call to SD_CalcOutput, we need to check that p_FAST%CompSub /= Module_SD before this:
! IF (p_FAST%CompSub /= Module_SD) THEN
CALL Transfer_SubStructureMotion_to_HD( SubstructureMotion2HD, u_HD%WAMITMesh, u_HD%Morison%Mesh, MeshMapData, ErrStat2, ErrMsg2 )
CALL SetErrStat(ErrStat2,ErrMsg2, ErrStat, ErrMsg, RoutineName)
!END IF ! don't transfer for SubDyn unless we have called SD_CalcOutput
END IF
ELSEIF ( p_FAST%CompMooring == Module_MD ) THEN
! motions:
CALL Transfer_Point_to_Point( SubstructureMotion, u_MD%CoupledKinematics(1), MeshMapData%Structure_2_Mooring, ErrStat2, ErrMsg2 )
CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat, ErrMsg,RoutineName//'u_MD%CoupledKinematics' )
```
### solveoption1
```fortran
SubstructureMotion => ED%y%PlatformPtMesh
ELSEIF ( p_FAST%CompMooring == Module_MD ) THEN
CALL MD_CalcOutput( this_time, MD%Input(1), MD%p, MD%x(this_state), MD%xd(this_state), MD%z(this_state), &
MD%OtherSt(this_state), MD%y, MD%m, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
ELSEIF ( p_FAST%SolveOption == Solve_SimplifiedOpt1 ) THEN ! No substructure model
CALL ED_HD_InputOutputSolve( this_time, p_FAST, calcJacobian &
, ED%Input(1), ED%p, ED%x(this_state), ED%xd(this_state), ED%z(this_state), ED%OtherSt(this_state), ED%y, ED%m &
, HD%Input(1), HD%p, HD%x(this_state), HD%xd(this_state), HD%z(this_state), HD%OtherSt(this_state), HD%y, HD%m &
, MAPp%Input(1), MAPp%y, FEAM%Input(1), FEAM%y, MD%Input(1), MD%y, SrvD%Input(1), SrvD%y &
, MeshMapData , ErrStat2, ErrMsg2, WriteThisStep )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
```
#### ED_HD_InputOutputSolve
```fortran
CALL ED_CopyOutput( y_ED, y_ED_input, MESH_NEWCOPY, ErrStat2, ErrMsg2)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
IF ( calcJacobian ) THEN
CALL ED_CopyInput( u_ED, u_ED_perturb, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL ED_CopyOutput( y_ED, y_ED_perturb, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL HydroDyn_CopyInput( u_HD, u_HD_perturb, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL HydroDyn_CopyOutput( y_HD, y_HD_perturb, MESH_NEWCOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL Transfer_PlatformMotion_to_HD(y_ED_input%PlatformPtMesh, u_HD, MeshMapData, ErrStat2, ErrMsg2 ) ! get u_HD from y_ED_input
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
u( 1: 3) = u_ED%PlatformPtMesh%Force(:,1) / p_FAST%UJacSclFact
u( 4: 6) = u_ED%PlatformPtMesh%Moment(:,1) / p_FAST%UJacSclFact
u( 7: 9) = y_ED_input%PlatformPtMesh%TranslationAcc(:,1)
u(10:12) = y_ED_input%PlatformPtMesh%RotationAcc(:,1)
K = 0
DO
!-------------------------------------------------------------------------------------------------
! Calculate outputs at this_time, based on inputs at this_time
!-------------------------------------------------------------------------------------------------
CALL ED_CalcOutput( this_time, u_ED, p_ED, x_ED, xd_ED, z_ED, OtherSt_ED, y_ED, m_ED, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL HydroDyn_CalcOutput( this_time, u_HD, p_HD, x_HD, xd_HD, z_HD, OtherSt_HD, y_HD, m_HD, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
IF ( K >= p_FAST%KMax ) EXIT
!-------------------------------------------------------------------------------------------------
! Calculate Jacobian: partial U/partial u:
! (note that we don't want to change u_ED or u_HD here)
!-------------------------------------------------------------------------------------------------
CALL U_ED_HD_Residual(y_ED, y_HD, u, Fn_U_Resid) ! U_ED_HD_Residual checks for error
IF (ErrStat >= AbortErrLev) THEN
CALL CleanUp()
RETURN
END IF
IF ( calcJacobian ) THEN
!...............................
! Get ElastoDyn's contribution:
!...............................
DO i=1,6 !call ED_CalcOutput
CALL ED_CopyInput( u_ED, u_ED_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
u_perturb = u
CALL Perturb_u( i, u_perturb, u_ED_perturb=u_ED_perturb, perturb=ThisPerturb ) ! perturb u and u_ED by ThisPerturb [routine sets ThisPerturb]
! calculate outputs with perturbed inputs:
CALL ED_CalcOutput( this_time, u_ED_perturb, p_ED, x_ED, xd_ED, z_ED, OtherSt_ED, y_ED_perturb, m_ED, ErrStat2, ErrMsg2 ) !calculate y_ED_perturb
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL U_ED_HD_Residual(y_ED_perturb, y_HD, u_perturb, Fn_U_perturb) ! get this perturbation, U_perturb
IF ( ErrStat >= AbortErrLev ) RETURN ! U_ED_HD_Residual checks for error
IF (ErrStat >= AbortErrLev) THEN
CALL CleanUp()
RETURN
END IF
MeshMapData%Jacobian_Opt1(:,i) = (Fn_U_perturb - Fn_U_Resid) / ThisPerturb
END DO ! ElastoDyn contribution ( columns 1-6 )
!...............................
! Get HydroDyn's contribution:
!...............................
DO i=7,12 !call HD_CalcOutput
! we want to perturb u_HD, but we're going to perturb the input y_ED and transfer that to HD to get u_HD
CALL ED_CopyOutput( y_ED_input, y_ED_perturb, MESH_UPDATECOPY, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
u_perturb = u
CALL Perturb_u( i, u_perturb, y_ED_perturb=y_ED_perturb, perturb=ThisPerturb ) ! perturb u and y_ED by ThisPerturb [routine sets ThisPerturb]
CALL Transfer_PlatformMotion_to_HD( y_ED_perturb%PlatformPtMesh, u_HD_perturb, MeshMapData, ErrStat2, ErrMsg2 ) ! get u_HD_perturb
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
! calculate outputs with perturbed inputs:
CALL HydroDyn_CalcOutput( this_time, u_HD_perturb, p_HD, x_HD, xd_HD, z_HD, OtherSt_HD, y_HD_perturb, m_HD, ErrStat2, ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL U_ED_HD_Residual(y_ED, y_HD_perturb, u_perturb, Fn_U_perturb) ! get this perturbation ! U_ED_HD_Residual checks for error
IF ( ErrStat >= AbortErrLev ) THEN
CALL CleanUp()
RETURN
END IF
MeshMapData%Jacobian_Opt1(:,i) = (Fn_U_perturb - Fn_U_Resid) / ThisPerturb
END DO ! HydroDyn contribution ( columns 7-12 )
CALL LAPACK_getrf( M=NumInputs, N=NumInputs, A=MeshMapData%Jacobian_Opt1, IPIV=MeshMapData%Jacobian_pivot, ErrStat=ErrStat2, ErrMsg=ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
u_delta = -Fn_U_Resid
CALL LAPACK_getrs( TRANS='N', N=NumInputs, A=MeshMapData%Jacobian_Opt1, IPIV=MeshMapData%Jacobian_pivot, B=u_delta, &
ErrStat=ErrStat2, ErrMsg=ErrMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
u = u + u_delta
u_ED%PlatformPtMesh%Force( :,1) = u_ED%PlatformPtMesh%Force( :,1) + u_delta( 1: 3) * p_FAST%UJacSclFact
u_ED%PlatformPtMesh%Moment(:,1) = u_ED%PlatformPtMesh%Moment(:,1) + u_delta( 4: 6) * p_FAST%UJacSclFact
y_ED_input%PlatformPtMesh%TranslationAcc(:,1) = y_ED_input%PlatformPtMesh%TranslationAcc(:,1) + u_delta( 7: 9)
y_ED_input%PlatformPtMesh%RotationAcc( :,1) = y_ED_input%PlatformPtMesh%RotationAcc( :,1) + u_delta(10:12)
CALL Transfer_PlatformMotion_to_HD( y_ED_input%PlatformPtMesh, u_HD, MeshMapData, ErrStat2, ErrMsg2 ) ! get u_HD with u_delta changes
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
K = K + 1
END DO ! K
```
CALL MD_CalcOutput