1508 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1509 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1510 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1511 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1512 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1513 double PSTAB = args.
scalar<
double>(
"PSTAB");
1514 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1515 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
1516 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1517 int nDOF_per_element_pressure = args.
scalar<
int>(
"nDOF_per_element_pressure");
1518 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
1519 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
1520 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
1521 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
1522 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
1523 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
1524 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
1525 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
1526 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
1527 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
1528 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
1529 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
1530 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
1531 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1532 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1533 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1534 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
1535 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
1536 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
1537 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
1538 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
1539 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
1540 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
1541 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
1542 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1543 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1544 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
1545 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1546 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1547 double hFactor = args.
scalar<
double>(
"hFactor");
1548 int nElements_global = args.
scalar<
int>(
"nElements_global");
1549 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
1550 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
1551 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
1552 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
1553 double useRBLES = args.
scalar<
double>(
"useRBLES");
1554 double useMetrics = args.
scalar<
double>(
"useMetrics");
1555 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1556 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
1557 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
1558 double sigma = args.
scalar<
double>(
"sigma");
1563 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
1564 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
1565 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
1566 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
1567 double C_dc = args.
scalar<
double>(
"C_dc");
1568 double C_b = args.
scalar<
double>(
"C_b");
1569 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
1570 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
1571 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
1572 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
1573 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
1574 xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
1575 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
1576 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
1577 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
1578 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
1579 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
1580 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
1581 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
1582 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
1583 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
1584 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
1585 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
1586 xt::pyarray<double>& q_eddy_viscosity = args.
array<
double>(
"q_eddy_viscosity");
1587 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
1588 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
1589 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
1590 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1591 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
1592 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
1593 xt::pyarray<double>& u_dof_old = args.
array<
double>(
"u_dof_old");
1594 xt::pyarray<double>& v_dof_old = args.
array<
double>(
"v_dof_old");
1595 xt::pyarray<double>& w_dof_old = args.
array<
double>(
"w_dof_old");
1596 xt::pyarray<double>& u_dof_old_old = args.
array<
double>(
"u_dof_old_old");
1597 xt::pyarray<double>& v_dof_old_old = args.
array<
double>(
"v_dof_old_old");
1598 xt::pyarray<double>& w_dof_old_old = args.
array<
double>(
"w_dof_old_old");
1599 xt::pyarray<double>& uStar_dof = args.
array<
double>(
"uStar_dof");
1600 xt::pyarray<double>& vStar_dof = args.
array<
double>(
"vStar_dof");
1601 xt::pyarray<double>& wStar_dof = args.
array<
double>(
"wStar_dof");
1602 xt::pyarray<double>& g = args.
array<
double>(
"g");
1603 const double useVF = args.
scalar<
double>(
"useVF");
1604 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
1605 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
1606 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
1607 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
1608 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
1609 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
1610 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
1611 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
1612 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
1613 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
1614 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
1615 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
1616 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
1617 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
1618 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
1619 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
1620 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
1621 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1622 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
1623 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
1624 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1625 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
1626 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
1627 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
1628 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
1629 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
1630 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
1631 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
1632 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
1633 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
1634 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
1635 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
1636 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
1637 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
1638 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
1639 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
1640 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
1641 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
1642 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
1643 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
1644 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
1645 int offset_p = args.
scalar<
int>(
"offset_p");
1646 int offset_u = args.
scalar<
int>(
"offset_u");
1647 int offset_v = args.
scalar<
int>(
"offset_v");
1648 int offset_w = args.
scalar<
int>(
"offset_w");
1649 int stride_p = args.
scalar<
int>(
"stride_p");
1650 int stride_u = args.
scalar<
int>(
"stride_u");
1651 int stride_v = args.
scalar<
int>(
"stride_v");
1652 int stride_w = args.
scalar<
int>(
"stride_w");
1653 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1654 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1655 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1656 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
1657 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1658 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1659 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
1660 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
1661 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
1662 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
1663 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
1664 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
1665 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
1666 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
1667 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
1668 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
1669 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1670 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
1671 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
1672 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
1673 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
1674 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
1675 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
1676 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
1677 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
1678 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
1679 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
1680 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
1681 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
1682 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
1683 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
1684 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1685 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
1686 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1687 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
1688 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
1689 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
1690 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
1691 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
1692 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
1693 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
1694 xt::pyarray<double>& q_grad_u = args.
array<
double>(
"q_grad_u");
1695 xt::pyarray<double>& q_grad_v = args.
array<
double>(
"q_grad_v");
1696 xt::pyarray<double>& q_grad_w = args.
array<
double>(
"q_grad_w");
1697 xt::pyarray<double>& q_divU = args.
array<
double>(
"q_divU");
1698 xt::pyarray<double>& ebqe_grad_u = args.
array<
double>(
"ebqe_grad_u");
1699 xt::pyarray<double>& ebqe_grad_v = args.
array<
double>(
"ebqe_grad_v");
1700 xt::pyarray<double>& ebqe_grad_w = args.
array<
double>(
"ebqe_grad_w");
1701 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
1702 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
1703 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
1704 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
1705 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
1706 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
1707 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
1708 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
1709 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
1710 xt::pyarray<double>& q_rho = args.
array<
double>(
"q_rho");
1711 xt::pyarray<double>& ebqe_rho = args.
array<
double>(
"ebqe_rho");
1712 xt::pyarray<double>& q_nu = args.
array<
double>(
"q_nu");
1713 xt::pyarray<double>& ebqe_nu = args.
array<
double>(
"ebqe_nu");
1714 int nParticles = args.
scalar<
int>(
"nParticles");
1715 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
1716 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
1717 double particle_beta = args.
scalar<
double>(
"particle_beta");
1718 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
1719 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
1720 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
1721 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
1722 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
1723 xt::pyarray<double>& particle_netForces = args.
array<
double>(
"particle_netForces");
1724 xt::pyarray<double>& particle_netMoments = args.
array<
double>(
"particle_netMoments");
1725 xt::pyarray<double>& particle_surfaceArea = args.
array<
double>(
"particle_surfaceArea");
1726 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
1727 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
1728 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
1729 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
1730 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
1731 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
1732 xt::pyarray<double>& phisError = args.
array<
double>(
"phisError");
1733 xt::pyarray<double>& phisErrorNodal = args.
array<
double>(
"phisErrorNodal");
1734 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
1735 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
1737 double cE = args.
scalar<
double>(
"cE");
1738 int MULTIPLY_EXTERNAL_FORCE_BY_DENSITY = args.
scalar<
int>(
"MULTIPLY_EXTERNAL_FORCE_BY_DENSITY");
1739 xt::pyarray<double>& forcex = args.
array<
double>(
"forcex");
1740 xt::pyarray<double>& forcey = args.
array<
double>(
"forcey");
1741 xt::pyarray<double>& forcez = args.
array<
double>(
"forcez");
1742 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
1743 double dt = args.
scalar<
double>(
"dt");
1744 xt::pyarray<double>& quantDOFs = args.
array<
double>(
"quantDOFs");
1745 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
1746 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
1747 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
1748 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
1749 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
1750 double order_polynomial = args.
scalar<
double>(
"order_polynomial");
1751 xt::pyarray<double>& isActiveDOF = args.
array<
double>(
"isActiveDOF");
1752 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
1753 xt::pyarray<double>& ncDrag = args.
array<
double>(
"ncDrag");
1754 xt::pyarray<double>& betaDrag = args.
array<
double>(
"betaDrag");
1755 xt::pyarray<double>& vos_vel_nodes = args.
array<
double>(
"vos_vel_nodes");
1756 xt::pyarray<double>& entropyResidualPerNode = args.
array<
double>(
"entropyResidualPerNode");
1757 xt::pyarray<double>& laggedEntropyResidualPerNode = args.
array<
double>(
"laggedEntropyResidualPerNode");
1758 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
1759 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
1760 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
1761 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
1762 int NNZ_1D = args.
scalar<
int>(
"NNZ_1D");
1763 xt::pyarray<int>& csrRowIndeces_1D = args.
array<
int>(
"csrRowIndeces_1D");
1764 xt::pyarray<int>& csrColumnOffsets_1D = args.
array<
int>(
"csrColumnOffsets_1D");
1765 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
1766 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
1767 xt::pyarray<double>& isBoundary_1D = args.
array<
double>(
"isBoundary_1D");
1768 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
1771 double element_uStar_He[nElements_global], element_vStar_He[nElements_global], element_wStar_He[nElements_global];
1775 den_hi.resize(numDOFs_1D,0.0);
1788 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1790 for (
int i=0; i<NNZ_1D; i++)
1792 uStar_dMatrix[i]=0.;
1793 vStar_dMatrix[i]=0.;
1794 wStar_dMatrix[i]=0.;
1798 for (
int i=0; i<numDOFs_1D; i++)
1803 entropyResidualPerNode[i]=0.;
1814 double mesh_volume_conservation=0.0,
1815 mesh_volume_conservation_weak=0.0,
1816 mesh_volume_conservation_err_max=0.0,
1817 mesh_volume_conservation_err_max_weak=0.0;
1818 double globalConservationError=0.0;
1819 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
1820 for(
int eN=0;eN<nElements_global;eN++)
1822 double elementTransport[nDOF_test_element][nDOF_trial_element];
1823 double elementTransposeTransport[nDOF_test_element][nDOF_trial_element];
1825 double elementResidual_p[nDOF_test_element],elementResidual_mesh[nDOF_test_element],
1826 elementResidual_u[nDOF_test_element],
1827 elementResidual_v[nDOF_test_element],
1828 mom_u_source_i[nDOF_test_element],
1829 mom_v_source_i[nDOF_test_element],
1830 mom_w_source_i[nDOF_test_element],
1831 betaDrag_i[nDOF_test_element],
1832 vos_i[nDOF_test_element],
1833 phisErrorElement[nDOF_test_element],
1834 elementResidual_w[nDOF_test_element],
1835 elementEntropyResidual[nDOF_test_element],
1837 double element_active=1.0;
1838 double mesh_volume_conservation_element=0.0,
1839 mesh_volume_conservation_element_weak=0.0;
1841 double linVisc_eN = 0, nlinVisc_eN_num = 0, nlinVisc_eN_den = 0;
1843 double det_hess_uStar_Ke=0.0, det_hess_vStar_Ke=0.0, det_hess_wStar_Ke=0.0, area_Ke=0.0;
1844 for (
int i=0;i<nDOF_test_element;i++)
1846 int eN_i = eN*nDOF_test_element+i;
1847 elementResidual_p_save[eN_i]=0.0;
1848 elementResidual_mesh[i]=0.0;
1849 elementResidual_p[i]=0.0;
1850 elementResidual_u[i]=0.0;
1851 elementResidual_v[i]=0.0;
1852 mom_u_source_i[i]=0.0;
1853 mom_v_source_i[i]=0.0;
1854 mom_w_source_i[i]=0.0;
1857 phisErrorElement[i]=0.0;
1858 elementResidual_w[i]=0.0;
1859 elementEntropyResidual[i]=0.0;
1860 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
1862 for (
int j=0;j<nDOF_trial_element;j++)
1864 elementTransport[i][j]=0.0;
1865 elementTransposeTransport[i][j]=0.0;
1870 if(use_ball_as_particle==1)
1872 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1874 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1875 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1876 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1877 phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]]);
1886 double _distance[nDOF_mesh_trial_element]={0.0};
1888 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1890 if(use_ball_as_particle==1)
1893 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
1894 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
1895 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
1900 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
1902 if ( _distance[I] >= 0)
1905 if (pos_counter == 3)
1909 for (
int I=0;I<nDOF_mesh_trial_element;I++)
1911 if (_distance[I] < 0)
1914 assert(opp_node >=0);
1915 assert(opp_node <nDOF_mesh_trial_element);
1920 const int ebN = elementBoundariesArray[eN*nDOF_mesh_trial_element+opp_node];
1921 const int eN_oppo = (eN == elementBoundaryElementsArray[ebN*2+0])?elementBoundaryElementsArray[ebN*2+1]:elementBoundaryElementsArray[ebN*2+0];
1922 if((mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]<nNodes_owned
1923 || mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]<nNodes_owned
1924 || mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]<nNodes_owned
1931 if (eN == elementBoundaryElementsArray[ebN*2+0])
1943 double distance=1e10, distance_to_ith_particle;
1944 if(use_ball_as_particle==1)
1946 double middle_point_coord[3]={0.0};
1947 double middle_point_distance;
1948 middle_point_coord[0] = (mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]+0]
1949 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]+0]
1950 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]+0])/3.0;
1951 middle_point_coord[1] = (mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]+1]
1952 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]+1]
1953 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]+1])/3.0;
1954 middle_point_coord[2] = (mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+1)%4]+2]
1955 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+2)%4]+2]
1956 +mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+(opp_node+3)%4]+2])/3.0;
1958 middle_point_coord[0],middle_point_coord[1],middle_point_coord[2],
1959 middle_point_distance);
1963 for (
int i=0;i<nParticles;++i)
1965 distance_to_ith_particle=particle_signed_distances[i*nElements_global*nQuadraturePoints_element
1966 +eN*nQuadraturePoints_element
1968 if (distance_to_ith_particle<distance)
1970 distance = distance_to_ith_particle;
1980 if(ebN<nElementBoundaries_owned)
1982 assert(eN_oppo==-1);
1986 else if (pos_counter == 4)
1989 for (
int i=0;i<nDOF_test_element;i++)
1991 isActiveDOF[offset_u+stride_u*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
1992 isActiveDOF[offset_v+stride_v*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
1993 isActiveDOF[offset_w+stride_w*vel_l2g[eN*nDOF_trial_element + i]]=1.0;
2001 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
2002 for (
int j=0;j<nDOF_mesh_trial_element;j++)
2004 int eN_j = eN*nDOF_mesh_trial_element+j;
2005 element_phi[j] = phi_dof[p_l2g[eN_j]];
2006 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
2008 double element_nodes[nDOF_mesh_trial_element*3];
2009 for (
int i=0;i<nDOF_mesh_trial_element;i++)
2011 int eN_i=eN*nDOF_mesh_trial_element+i;
2012 for(
int I=0;I<3;I++)
2013 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
2015 gf_s.
calculate(element_phi_s, element_nodes, x_ref.data(),
false);
2016 gf.
calculate(element_phi, element_nodes, x_ref.data(),
false);
2020 for(
int k=0;k<nQuadraturePoints_element;k++)
2025 int eN_k = eN*nQuadraturePoints_element+k,
2026 eN_k_nSpace = eN_k*nSpace,
2028 eN_nDOF_trial_element = eN*nDOF_trial_element;
2029 double p=0.0,
u=0.0,
v=0.0,
w=0.0,un=0.0,vn=0.0,wn=0.0,
2030 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
2039 dmass_adv_u[nSpace],
2040 dmass_adv_v[nSpace],
2041 dmass_adv_w[nSpace],
2043 dmom_u_adv_u[nSpace],
2044 dmom_u_adv_v[nSpace],
2045 dmom_u_adv_w[nSpace],
2047 dmom_v_adv_u[nSpace],
2048 dmom_v_adv_v[nSpace],
2049 dmom_v_adv_w[nSpace],
2051 dmom_w_adv_u[nSpace],
2052 dmom_w_adv_v[nSpace],
2053 dmom_w_adv_w[nSpace],
2054 mom_uu_diff_ten[nSpace],
2055 mom_vv_diff_ten[nSpace],
2056 mom_ww_diff_ten[nSpace],
2067 dmom_u_ham_grad_p[nSpace],
2068 dmom_u_ham_grad_u[nSpace],
2070 dmom_v_ham_grad_p[nSpace],
2071 dmom_v_ham_grad_v[nSpace],
2073 dmom_w_ham_grad_p[nSpace],
2074 dmom_w_ham_grad_w[nSpace],
2085 Lstar_u_p[nDOF_test_element],
2086 Lstar_v_p[nDOF_test_element],
2087 Lstar_w_p[nDOF_test_element],
2088 Lstar_u_u[nDOF_test_element],
2089 Lstar_v_v[nDOF_test_element],
2090 Lstar_w_w[nDOF_test_element],
2091 Lstar_p_u[nDOF_test_element],
2092 Lstar_p_v[nDOF_test_element],
2093 Lstar_p_w[nDOF_test_element],
2098 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
2099 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
2102 jacInv[nSpace*nSpace],
2103 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
2104 vel_hess_trial[nDOF_trial_element*
nSpace2],
2105 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
2106 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
2107 u_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2108 v_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2109 w_times_vel_grad_test_dV[nDOF_test_element*nSpace],
2115 dmom_u_source[nSpace],
2116 dmom_v_source[nSpace],
2117 dmom_w_source[nSpace],
2121 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
2123 ck.calculateMapping_element(eN,
2127 mesh_trial_ref.data(),
2128 mesh_grad_trial_ref.data(),
2133 ck.calculateH_element(eN,
2135 nodeDiametersArray.data(),
2137 mesh_trial_ref.data(),
2139 ck.calculateMappingVelocity_element(eN,
2141 mesh_velocity_dof.data(),
2143 mesh_trial_ref.data(),
2148 dV = fabs(jacDet)*dV_ref[k];
2149 ck.calculateG(jacInv,G,G_dd_G,tr_G);
2152 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2153 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2154 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
2158 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
2159 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
2164 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
2165 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
2166 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
w);
2168 ck.valFromDOF(u_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],un);
2169 ck.valFromDOF(v_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],vn);
2170 ck.valFromDOF(w_dof_old.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],wn);
2173 for (
int I=0;I<nSpace;I++)
2174 grad_p[I] = q_grad_p[eN_k_nSpace + I];
2175 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
2176 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
2177 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_w);
2178 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
2179 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
2180 ck.hessFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_w);
2181 ck.hessFromDOF(uStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_uStar);
2182 ck.hessFromDOF(vStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_vStar);
2183 ck.hessFromDOF(wStar_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_wStar);
2185 for (
int j=0;j<nDOF_trial_element;j++)
2188 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
2189 for (
int I=0;I<nSpace;I++)
2192 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
2193 if (ARTIFICIAL_VISCOSITY==4)
2196 u_times_vel_grad_test_dV[j*nSpace+I] =
2197 u*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_u[I];
2198 v_times_vel_grad_test_dV[j*nSpace+I] =
2199 v*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_v[I];
2200 w_times_vel_grad_test_dV[j*nSpace+I] =
2201 w*vel_grad_trial[j*nSpace+I]*dV + vel_test_dV[j]*grad_w[I];
2206 if (ARTIFICIAL_VISCOSITY==3)
2208 det_hess_uStar_Ke +=
2209 (hess_uStar[0]*(hess_uStar[4]*hess_uStar[8] - hess_uStar[7]*hess_uStar[5])
2210 -hess_uStar[1]*(hess_uStar[3]*hess_uStar[8] - hess_uStar[6]*hess_uStar[5])
2211 +hess_uStar[2]*(hess_uStar[3]*hess_uStar[7] - hess_uStar[6]*hess_uStar[4]))*dV;
2212 det_hess_vStar_Ke +=
2213 (hess_vStar[0]*(hess_vStar[4]*hess_vStar[8] - hess_vStar[7]*hess_vStar[5])
2214 -hess_vStar[1]*(hess_vStar[3]*hess_vStar[8] - hess_vStar[6]*hess_vStar[5])
2215 +hess_vStar[2]*(hess_vStar[3]*hess_vStar[7] - hess_vStar[6]*hess_vStar[4]))*dV;
2216 det_hess_wStar_Ke +=
2217 (hess_wStar[0]*(hess_wStar[4]*hess_wStar[8] - hess_wStar[7]*hess_wStar[5])
2218 -hess_wStar[1]*(hess_wStar[3]*hess_wStar[8] - hess_wStar[6]*hess_wStar[5])
2219 +hess_wStar[2]*(hess_wStar[3]*hess_wStar[7] - hess_wStar[6]*hess_wStar[4]))*dV;
2223 double div_mesh_velocity=0.0;
2224 int NDOF_MESH_TRIAL_ELEMENT=4;
2225 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
2227 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
2228 div_mesh_velocity +=
2229 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*nSpace+0] +
2230 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*nSpace+1] +
2231 mesh_velocity_dof[mesh_l2g[eN_j]*3+2]*vel_grad_trial[j*nSpace+2];
2233 mesh_volume_conservation_element += (alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity)*dV;
2234 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
2236 porosity = 1.0 - q_vos[eN_k];
2245 double distance_to_omega_solid = 1e10;
2246 if(use_ball_as_particle==1)
2250 distance_to_omega_solid);
2254 for (
int i = 0; i < nParticles; i++)
2256 double distance_to_i_th_solid = particle_signed_distances[i * nElements_global * nQuadraturePoints_element + eN_k];
2257 distance_to_omega_solid = (distance_to_i_th_solid < distance_to_omega_solid)?distance_to_i_th_solid:distance_to_omega_solid;
2260 phi_solid[eN_k] = distance_to_omega_solid;
2272 elementDiameter[eN],
2273 smagorinskyConstant,
2274 turbulenceClosureModel,
2279 &normal_phi[eN_k_nSpace],
2280 distance_to_omega_solid,
2293 q_velocity_sge[eN_k_nSpace+0],
2294 q_velocity_sge[eN_k_nSpace+1],
2295 q_velocity_sge[eN_k_nSpace+2],
2296 q_eddy_viscosity[eN_k],
2343 MULTIPLY_EXTERNAL_FORCE_BY_DENSITY,
2347 MATERIAL_PARAMETERS_AS_FUNCTION,
2348 density_as_function[eN_k],
2349 dynamic_viscosity_as_function[eN_k],
2352 use_ball_as_particle,
2355 ball_velocity.data(),
2356 ball_angular_velocity.data(),
2357 INT_BY_PARTS_PRESSURE);
2360 mass_source = q_mass_source[eN_k];
2361 for (
int I=0;I<nSpace;I++)
2363 dmom_u_source[I] = 0.0;
2364 dmom_v_source[I] = 0.0;
2365 dmom_w_source[I] = 0.0;
2376 q_eddy_viscosity[eN_k],
2383 q_velocity_sge[eN_k_nSpace+0],
2384 q_velocity_sge[eN_k_nSpace+1],
2385 q_velocity_sge[eN_k_nSpace+2],
2386 eps_solid[elementFlags[eN]],
2388 q_velocity_solid[eN_k_nSpace+0],
2389 q_velocity_solid[eN_k_nSpace+1],
2390 q_velocity_solid[eN_k_nSpace+2],
2391 q_velocityStar_solid[eN_k_nSpace+0],
2392 q_velocityStar_solid[eN_k_nSpace+1],
2393 q_velocityStar_solid[eN_k_nSpace+2],
2400 q_grad_vos[eN_k_nSpace+0],
2401 q_grad_vos[eN_k_nSpace+1],
2402 q_grad_vos[eN_k_nSpace+2]);
2403 double C_particles = 0.0;
2404 if (nParticles > 0 && USE_SBM==0)
2409 nQuadraturePoints_global,
2410 &particle_signed_distances[eN_k],
2411 &particle_signed_distance_normals[eN_k_nSpace],
2412 &particle_velocities[eN_k_nSpace],
2413 particle_centroids.data(),
2414 use_ball_as_particle,
2417 ball_velocity.data(),
2418 ball_angular_velocity.data(),
2420 particle_penalty_constant/h_phi,
2421 particle_alpha/h_phi,
2422 particle_beta/h_phi,
2439 q_velocity_sge[eN_k_nSpace + 0],
2440 q_velocity_sge[eN_k_nSpace + 1],
2441 q_velocity_sge[eN_k_nSpace + 2],
2464 particle_netForces.data(),
2465 particle_netMoments.data(),
2466 particle_surfaceArea.data());
2468 if (turbulenceClosureModel >= 3)
2470 const double c_mu = 0.09;
2485 &q_turb_var_grad_0[eN_k_nSpace],
2486 q_eddy_viscosity[eN_k],
2504 q_mom_u_acc[eN_k] = mom_u_acc;
2505 q_mom_v_acc[eN_k] = mom_v_acc;
2506 q_mom_w_acc[eN_k] = mom_w_acc;
2508 q_mass_adv[eN_k_nSpace+0] =
u;
2509 q_mass_adv[eN_k_nSpace+1] =
v;
2510 q_mass_adv[eN_k_nSpace+2] =
w;
2514 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
2515 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
2516 mom_u_adv[2] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*zt;
2517 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
2518 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
2519 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
2521 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
2522 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
2523 mom_v_adv[2] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*zt;
2524 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
2525 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
2526 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
2528 mom_w_adv[0] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*
xt;
2529 mom_w_adv[1] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*yt;
2530 mom_w_adv[2] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*zt;
2531 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
2532 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
2533 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
2538 if (q_dV_last[eN_k] <= -100)
2539 q_dV_last[eN_k] = dV;
2542 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2548 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2554 q_mom_w_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
2560 mom_u_acc_t *= dmom_u_acc_u;
2561 mom_v_acc_t *= dmom_v_acc_v;
2562 mom_w_acc_t *= dmom_w_acc_w;
2570 ck.Mass_strong(-q_dvos_dt[eN_k]) +
2571 ck.Advection_strong(dmass_adv_u,grad_u) +
2572 ck.Advection_strong(dmass_adv_v,grad_v) +
2573 ck.Advection_strong(dmass_adv_w,grad_w) +
2574 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
2576 ck.Reaction_strong(mass_source);
2579 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
2580 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
2581 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
2584 ck.Mass_strong(mom_u_acc_t) +
2585 ck.Advection_strong(dmom_adv_sge,grad_u) +
2586 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
2587 ck.Reaction_strong(mom_u_source) -
2588 ck.Reaction_strong(
u*div_mesh_velocity);
2591 ck.Mass_strong(mom_v_acc_t) +
2592 ck.Advection_strong(dmom_adv_sge,grad_v) +
2593 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
2594 ck.Reaction_strong(mom_v_source) -
2595 ck.Reaction_strong(
v*div_mesh_velocity);
2598 ck.Mass_strong(mom_w_acc_t) +
2599 ck.Advection_strong(dmom_adv_sge,grad_w) +
2600 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
2601 ck.Reaction_strong(mom_w_source) -
2602 ck.Reaction_strong(
w*div_mesh_velocity);
2606 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
2608 elementDiameter[eN],
2613 dmom_u_ham_grad_p[0],
2623 dmom_u_ham_grad_p[0],
2628 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
2629 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
2642 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
2643 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
2644 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
2646 mom_u_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge[eN_k_nSpace+0]);
2647 mom_u_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge[eN_k_nSpace+1]);
2648 mom_u_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge[eN_k_nSpace+2]);
2651 for (
int i=0;i<nDOF_test_element;i++)
2653 int i_nSpace = i*nSpace;
2658 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2659 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2660 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2661 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2662 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2663 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2666 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
2667 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
2668 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
2672 if (ARTIFICIAL_VISCOSITY==0 || ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2674 q_numDiff_u[eN_k] = 0;
2675 q_numDiff_v[eN_k] = 0;
2676 q_numDiff_w[eN_k] = 0;
2678 else if (ARTIFICIAL_VISCOSITY==1)
2680 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v + pdeResidual_w*pdeResidual_w);
2681 q_numDiff_u[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
2682 (1.0-useMetrics)*hFactor*hFactor*elementDiameter[eN]*elementDiameter[eN]);
2683 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
2684 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
2688 double rho = q_rho[eN_k];
2689 double mu = q_rho[eN_k]*q_nu[eN_k];
2691 double vel2 =
u*
u +
v*
v +
w*
w;
2695 porosity*rho*((
u-un)/dt + (
u*grad_u[0]+
v*grad_u[1]+
w*grad_u[2]) - g[0])
2696 + (KILL_PRESSURE_TERM == 1 ? 0 : 1.)*grad_p[0]
2697 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcex[eN_k]
2698 - mu*(hess_u[0] + hess_u[4] + hess_u[8])
2699 - mu*(hess_u[0] + hess_v[1] + hess_w[2]);
2701 porosity*rho*((
v-vn)/dt + (
u*grad_v[0]+
v*grad_v[1]+
w*grad_v[2]) - g[1])
2702 + (KILL_PRESSURE_TERM == 1 ? 0 : 1.)*grad_p[1]
2703 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcey[eN_k]
2704 - mu*(hess_v[0] + hess_v[4] + hess_v[8])
2705 - mu*(hess_u[1] + hess_v[4] + hess_w[5]);
2707 porosity*rho*((
w-wn)/dt + (
u*grad_w[0]+
v*grad_w[1]+
w*grad_w[2]) - g[2])
2708 + (KILL_PRESSURE_TERM == 1 ? 0 : 1.)*grad_p[2]
2709 - (MULTIPLY_EXTERNAL_FORCE_BY_DENSITY == 1 ? porosity*rho : 1.0)*forcez[eN_k]
2710 - mu*(hess_w[0] + hess_w[4] + hess_w[8])
2711 - mu*(hess_u[2] + hess_v[5] + hess_w[8]);
2714 double entRes_times_u = Res_in_x*
u + Res_in_y*
v + Res_in_z*
w;
2716 double hK = elementDiameter[eN]/order_polynomial;
2717 q_numDiff_u[eN_k] = fmin(
cMax*porosity*rho*hK*std::sqrt(vel2),
2718 cE*hK*hK*fabs(entRes_times_u)/(vel2+1E-10));
2719 q_numDiff_v[eN_k] = q_numDiff_u[eN_k];
2720 q_numDiff_w[eN_k] = q_numDiff_u[eN_k];
2724 linVisc_eN = fmax(porosity*rho*std::sqrt(vel2),linVisc_eN);
2725 nlinVisc_eN_num = fmax(fabs(entRes_times_u),nlinVisc_eN_num);
2726 nlinVisc_eN_den = fmax(vel2,nlinVisc_eN_den);
2737 q_velocity[eN_k_nSpace+0]=
u;
2738 q_velocity[eN_k_nSpace+1]=
v;
2739 q_velocity[eN_k_nSpace+2]=
w;
2740 for (
int I=0;I<nSpace;I++)
2742 q_grad_u[eN_k_nSpace+I] = grad_u[I];
2743 q_grad_v[eN_k_nSpace+I] = grad_v[I];
2744 q_grad_w[eN_k_nSpace+I] = grad_w[I];
2747 q_divU[eN_k] = q_grad_u[eN_k_nSpace+0] + q_grad_v[eN_k_nSpace+1] + q_grad_w[eN_k_nSpace+2];
2750 double unit_normal[nSpace];
2751 double norm_grad_phi = 0.;
2752 for (
int I=0;I<nSpace;I++)
2753 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
2754 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
2755 for (
int I=0;I<nSpace;I++)
2756 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
2760 v1[0]=1.-unit_normal[0]*unit_normal[0];
2761 v1[1]=-unit_normal[0]*unit_normal[1];
2762 v1[2]=-unit_normal[0]*unit_normal[2];
2765 v2[0]=-unit_normal[0]*unit_normal[1];
2766 v2[1]=1.-unit_normal[1]*unit_normal[1];
2767 v2[2]=-unit_normal[1]*unit_normal[2];
2770 v3[0]=-unit_normal[0]*unit_normal[2];
2771 v3[1]=-unit_normal[1]*unit_normal[2];
2772 v3[2]=1.-unit_normal[2]*unit_normal[2];
2773 double delta =
gf.
D(eps_mu,
phi[eN_k]);
2774 double vel_tgrad_test_i[nSpace],
2775 tgrad_u[nSpace], tgrad_v[nSpace], tgrad_w[nSpace];
2787 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2789 velStar[0] = q_velocity_sge[eN_k_nSpace+0];
2790 velStar[1] = q_velocity_sge[eN_k_nSpace+1];
2791 velStar[2] = q_velocity_sge[eN_k_nSpace+2];
2793 for(
int i=0;i<nDOF_test_element;i++)
2795 int i_nSpace=i*nSpace;
2797 &vel_grad_trial[i_nSpace],
2799 phisErrorElement[i]+=std::abs(phisError[eN_k_nSpace+0])*p_test_dV[i];
2817 elementResidual_u[i] +=
2818 ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
2819 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
2820 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
2821 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
2822 ck.Diffusion_weak(sdInfo_u_w_rowptr.data(),sdInfo_u_w_colind.data(),mom_uw_diff_ten,grad_w,&vel_grad_test_dV[i_nSpace]) +
2823 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
2824 ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
2825 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+0] : 0.) +
2827 USE_SUPG*
ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
2828 ck.NumericalDiffusion(q_numDiff_u_last[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]) +
2830 ck.NumericalDiffusion(delta*sigma*dV,v1,vel_tgrad_test_i) +
2831 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_u,vel_tgrad_test_i);
2832 mom_u_source_i[i] +=
ck.Reaction_weak(mom_u_source,vel_test_dV[i]);
2833 betaDrag_i[i] +=
ck.Reaction_weak(dmom_u_source[0],
2835 vos_i[i] +=
ck.Reaction_weak(1.0-porosity,
2838 elementResidual_v[i] +=
2839 ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
2840 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
2841 ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
2842 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
2843 ck.Diffusion_weak(sdInfo_v_w_rowptr.data(),sdInfo_v_w_colind.data(),mom_vw_diff_ten,grad_w,&vel_grad_test_dV[i_nSpace]) +
2844 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
2845 ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
2846 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+1] : 0.) +
2848 USE_SUPG*
ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
2849 ck.NumericalDiffusion(q_numDiff_v_last[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]) +
2851 ck.NumericalDiffusion(delta*sigma*dV,v2,vel_tgrad_test_i) +
2852 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_v,vel_tgrad_test_i);
2853 mom_v_source_i[i] +=
ck.Reaction_weak(mom_v_source,vel_test_dV[i]);
2855 elementResidual_w[i] +=
2856 ck.Mass_weak(mom_w_acc_t,vel_test_dV[i]) +
2857 ck.Advection_weak(mom_w_adv,&vel_grad_test_dV[i_nSpace]) +
2858 ck.Diffusion_weak(sdInfo_w_u_rowptr.data(),sdInfo_w_u_colind.data(),mom_wu_diff_ten,grad_u,&vel_grad_test_dV[i_nSpace]) +
2859 ck.Diffusion_weak(sdInfo_w_v_rowptr.data(),sdInfo_w_v_colind.data(),mom_wv_diff_ten,grad_v,&vel_grad_test_dV[i_nSpace]) +
2860 ck.Diffusion_weak(sdInfo_w_w_rowptr.data(),sdInfo_w_w_colind.data(),mom_ww_diff_ten,grad_w,&vel_grad_test_dV[i_nSpace]) +
2861 ck.Reaction_weak(mom_w_source,vel_test_dV[i]) +
2862 ck.Hamiltonian_weak(mom_w_ham,vel_test_dV[i]) +
2863 (INT_BY_PARTS_PRESSURE==1 ? -1.0*p*vel_grad_test_dV[i_nSpace+2] : 0.) +
2865 USE_SUPG*
ck.SubgridError(subgridError_w,Lstar_w_w[i]) +
2866 ck.NumericalDiffusion(q_numDiff_w_last[eN_k],grad_w,&vel_grad_test_dV[i_nSpace]) +
2868 ck.NumericalDiffusion(delta*sigma*dV,v3,vel_tgrad_test_i) +
2869 ck.NumericalDiffusion(dt*delta*sigma*dV,tgrad_w,vel_tgrad_test_i);
2870 mom_w_source_i[i] +=
ck.Reaction_weak(mom_w_source,vel_test_dV[i]);
2872 if (ARTIFICIAL_VISCOSITY==4)
2876 elementEntropyResidual[i] +=
2878 ck.Mass_weak(mom_u_acc_t,
u*vel_test_dV[i]) +
2879 ck.Advection_weak(mom_u_adv,&u_times_vel_grad_test_dV[i_nSpace])+
2880 ck.Diffusion_weak(sdInfo_u_u_rowptr.data(),
2881 sdInfo_u_u_colind.data(),
2884 &u_times_vel_grad_test_dV[i_nSpace]) +
2885 ck.Diffusion_weak(sdInfo_u_v_rowptr.data(),
2886 sdInfo_u_v_colind.data(),
2889 &u_times_vel_grad_test_dV[i_nSpace]) +
2890 ck.Diffusion_weak(sdInfo_u_w_rowptr.data(),
2891 sdInfo_u_w_colind.data(),
2894 &u_times_vel_grad_test_dV[i_nSpace]) +
2895 ck.Reaction_weak(mom_u_source,
u*vel_test_dV[i]) +
2896 ck.Hamiltonian_weak(mom_u_ham,
u*vel_test_dV[i])
2898 ck.Mass_weak(mom_v_acc_t,
v*vel_test_dV[i]) +
2899 ck.Advection_weak(mom_v_adv,&v_times_vel_grad_test_dV[i_nSpace])+
2900 ck.Diffusion_weak(sdInfo_v_u_rowptr.data(),
2901 sdInfo_v_u_colind.data(),
2904 &v_times_vel_grad_test_dV[i_nSpace])+
2905 ck.Diffusion_weak(sdInfo_v_v_rowptr.data(),
2906 sdInfo_v_v_colind.data(),
2909 &v_times_vel_grad_test_dV[i_nSpace])+
2910 ck.Diffusion_weak(sdInfo_v_w_rowptr.data(),
2911 sdInfo_v_w_colind.data(),
2914 &v_times_vel_grad_test_dV[i_nSpace]) +
2915 ck.Reaction_weak(mom_v_source,
v*vel_test_dV[i]) +
2916 ck.Hamiltonian_weak(mom_v_ham,
v*vel_test_dV[i])
2918 ck.Mass_weak(mom_w_acc_t,
w*vel_test_dV[i]) +
2919 ck.Advection_weak(mom_w_adv,&w_times_vel_grad_test_dV[i_nSpace])+
2920 ck.Diffusion_weak(sdInfo_w_u_rowptr.data(),
2921 sdInfo_w_u_colind.data(),
2924 &w_times_vel_grad_test_dV[i_nSpace]) +
2925 ck.Diffusion_weak(sdInfo_w_v_rowptr.data(),
2926 sdInfo_w_v_colind.data(),
2929 &w_times_vel_grad_test_dV[i_nSpace]) +
2930 ck.Diffusion_weak(sdInfo_w_w_rowptr.data(),
2931 sdInfo_w_w_colind.data(),
2934 &w_times_vel_grad_test_dV[i_nSpace]) +
2935 ck.Reaction_weak(mom_w_source,
w*vel_test_dV[i]) +
2936 ck.Hamiltonian_weak(mom_w_ham,
w*vel_test_dV[i]);
2938 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
2940 for(
int j=0;j<nDOF_trial_element;j++)
2942 int j_nSpace = j*nSpace;
2943 int i_nSpace = i*nSpace;
2944 elementTransport[i][j] +=
2945 q_rho[eN_k]*porosity*
2946 ck.AdvectionJacobian_strong(velStar,
2947 &vel_grad_test_dV[j_nSpace])
2948 *vel_trial_ref[k*nDOF_trial_element+i];
2949 elementTransposeTransport[i][j] +=
2950 q_rho[eN_k]*porosity*
2951 ck.AdvectionJacobian_strong(velStar,
2952 &vel_grad_test_dV[i_nSpace])
2953 *vel_trial_ref[k*nDOF_trial_element+j];
2958 element_uStar_He[eN] = det_hess_uStar_Ke/area_Ke;
2959 element_vStar_He[eN] = det_hess_vStar_Ke/area_Ke;
2960 element_wStar_He[eN] = det_hess_wStar_Ke/area_Ke;
2965 double hK = elementDiameter[eN];
2966 double artVisc = fmin(
cMax*hK*linVisc_eN,
2967 cE*hK*hK*nlinVisc_eN_num/(nlinVisc_eN_den+1E-10));
2968 for(
int k=0;k<nQuadraturePoints_element;k++)
2970 int eN_k = eN*nQuadraturePoints_element+k;
2971 q_numDiff_u[eN_k] = artVisc;
2972 q_numDiff_v[eN_k] = artVisc;
2973 q_numDiff_w[eN_k] = artVisc;
2979 for(
int i=0;i<nDOF_test_element;i++)
2981 int eN_i=eN*nDOF_test_element+i;
2982 phisErrorNodal[vel_l2g[eN_i]]+= element_active*phisErrorElement[i];
2986 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=element_active*elementResidual_u[i];
2987 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=element_active*elementResidual_v[i];
2988 globalResidual[offset_w+stride_w*vel_l2g[eN_i]]+=element_active*elementResidual_w[i];
2989 ncDrag[offset_u+stride_u*vel_l2g[eN_i]]+=mom_u_source_i[i];
2990 ncDrag[offset_v+stride_v*vel_l2g[eN_i]]+=mom_v_source_i[i];
2991 ncDrag[offset_w+stride_w*vel_l2g[eN_i]]+=mom_w_source_i[i];
2992 betaDrag[vel_l2g[eN_i]] += betaDrag_i[i];
2993 vos_vel_nodes[vel_l2g[eN_i]] += vos_i[i];
2996 if (ARTIFICIAL_VISCOSITY==3)
2998 uStar_hi[vel_l2g[eN_i]] += element_uStar_He[eN];
2999 vStar_hi[vel_l2g[eN_i]] += element_vStar_He[eN];
3000 wStar_hi[vel_l2g[eN_i]] += element_wStar_He[eN];
3001 den_hi[vel_l2g[eN_i]] += 1;
3003 if (ARTIFICIAL_VISCOSITY==4)
3006 entropyResidualPerNode[vel_l2g[eN_i]] += elementEntropyResidual[i];
3008 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3010 for (
int j=0;j<nDOF_trial_element;j++)
3012 int eN_i_j = eN_i*nDOF_trial_element+j;
3014 + csrColumnOffsets_1D[eN_i_j]]
3015 += elementTransport[i][j];
3018 + csrColumnOffsets_1D[eN_i_j]]
3019 += elementTransposeTransport[i][j];
3030 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
3033 for (
int i=0; i<numDOFs_1D; i++)
3035 if (ARTIFICIAL_VISCOSITY==4)
3038 double max_u2i = (std::pow(u_dof[i],2.) +
3039 std::pow(v_dof[i],2.) +
3040 std::pow(w_dof[i],2.));
3041 double min_u2i = max_u2i;
3042 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3044 int j = colind_1D[offset];
3045 double u2j = (std::pow(u_dof[j],2.) +
3046 std::pow(v_dof[j],2.) +
3047 std::pow(w_dof[j],2.));
3048 max_u2i = fmax(max_u2i,u2j);
3049 min_u2i = fmin(min_u2i,u2j);
3051 double normi = 0.5*(max_u2i + min_u2i) + 1E-10;
3052 entropyResidualPerNode[i] = fabs(entropyResidualPerNode[i])/normi;
3057 double uStari = uStar_dof[i];
3058 double vStari = vStar_dof[i];
3059 double wStari = wStar_dof[i];
3061 double u_beta_numerator = 0., u_beta_denominator = 0.;
3062 double v_beta_numerator = 0., v_beta_denominator = 0.;
3063 double w_beta_numerator = 0., w_beta_denominator = 0.;
3066 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3068 int j = colind_1D[offset];
3069 double uStarj = uStar_dof[j];
3070 double vStarj = vStar_dof[j];
3071 double wStarj = wStar_dof[j];
3074 u_beta_numerator += (uStarj - uStari);
3075 u_beta_denominator += fabs(uStarj - uStari);
3077 v_beta_numerator += (vStarj - vStari);
3078 v_beta_denominator += fabs(vStarj - vStari);
3080 w_beta_numerator += (wStarj - wStari);
3081 w_beta_denominator += fabs(wStarj - wStari);
3083 double u_beta = fabs(u_beta_numerator)/(u_beta_denominator+1E-10);
3084 double v_beta = fabs(v_beta_numerator)/(v_beta_denominator+1E-10);
3085 double w_beta = fabs(w_beta_numerator)/(w_beta_denominator+1E-10);
3107 if (ARTIFICIAL_VISCOSITY==3)
3109 for(
int eN=0;eN<nElements_global;eN++)
3111 double uStar_He = element_uStar_He[eN];
3112 double vStar_He = element_vStar_He[eN];
3113 double wStar_He = element_wStar_He[eN];
3114 for(
int i=0;i<nDOF_test_element;i++)
3116 int eN_i=eN*nDOF_test_element+i;
3117 int gi = vel_l2g[eN_i];
3126 if (ARTIFICIAL_VISCOSITY==3)
3128 for (
int i=0; i<numDOFs_1D; i++)
3134 if (isBoundary_1D[i] == 1)
3162 for (
int i=0; i<numDOFs_1D; i++)
3165 double uStar_dii = 0;
3166 double vStar_dii = 0;
3167 double wStar_dii = 0;
3168 double ui = u_dof[i];
3169 double vi = v_dof[i];
3170 double wi = w_dof[i];
3172 double ith_u_dissipative_term = 0;
3173 double ith_v_dissipative_term = 0;
3174 double ith_w_dissipative_term = 0;
3180 for (
int offset=rowptr_1D[i]; offset<rowptr_1D[i+1]; offset++)
3182 int j = colind_1D[offset];
3185 double uj = u_dof[j];
3186 double vj = v_dof[j];
3187 double wj = w_dof[j];
3193 if (ARTIFICIAL_VISCOSITY==4)
3195 double dEVij = fmax(laggedEntropyResidualPerNode[i],
3196 laggedEntropyResidualPerNode[j]);
3199 uStar_dMatrix[ij] = fmin(dLij,
cE*dEVij);
3200 vStar_dMatrix[i] = uStar_dMatrix[ij];
3201 wStar_dMatrix[i] = uStar_dMatrix[ij];
3212 uStar_dii -= uStar_dMatrix[ij];
3213 vStar_dii -= vStar_dMatrix[ij];
3214 wStar_dii -= wStar_dMatrix[ij];
3216 ith_u_dissipative_term += uStar_dMatrix[ij]*(uj-ui);
3217 ith_v_dissipative_term += vStar_dMatrix[ij]*(vj-vi);
3218 ith_w_dissipative_term += wStar_dMatrix[ij]*(wj-wi);
3227 uStar_dMatrix[ii] = uStar_dii;
3228 vStar_dMatrix[ii] = vStar_dii;
3229 wStar_dMatrix[ii] = wStar_dii;
3230 globalResidual[offset_u+stride_u*i] += -ith_u_dissipative_term;
3231 globalResidual[offset_v+stride_v*i] += -ith_v_dissipative_term;
3232 globalResidual[offset_w+stride_w*i] += -ith_w_dissipative_term;
3241 std::memset(particle_netForces.data(),0,nParticles*3*
sizeof(
double));
3242 std::memset(particle_netMoments.data(),0,nParticles*3*
sizeof(
double));
3248 eN_nDOF_trial_element = eN*nDOF_trial_element;
3249 double elementResidual_mesh[nDOF_test_element],
3250 elementResidual_p[nDOF_test_element],
3251 elementResidual_u[nDOF_test_element],
3252 elementResidual_v[nDOF_test_element],
3253 elementResidual_w[nDOF_test_element],
3258 for (
int i=0;i<nDOF_test_element;i++)
3260 elementResidual_mesh[i]=0.0;
3261 elementResidual_p[i]=0.0;
3262 elementResidual_u[i]=0.0;
3263 elementResidual_v[i]=0.0;
3264 elementResidual_w[i]=0.0;
3266 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3268 int ebN_kb = ebN*nQuadraturePoints_elementBoundary+kb,
3270 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3271 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3272 double u_ext=0.0, v_ext=0.0, w_ext=0.0,
3273 bc_u_ext=0.0, bc_v_ext=0.0, bc_w_ext=0.0,
3274 grad_u_ext[nSpace], grad_v_ext[nSpace], grad_w_ext[nSpace],
3275 jac_ext[nSpace*nSpace],
3277 jacInv_ext[nSpace*nSpace],
3278 boundaryJac[nSpace*(nSpace-1)],
3279 metricTensor[(nSpace-1)*(nSpace-1)],
3280 metricTensorDetSqrt,
3282 p_test_dS[nDOF_test_element],p_grad_trial_trace[nDOF_trial_element*nSpace],
3283 vel_test_dS[nDOF_test_element],
3284 vel_grad_trial_trace[nDOF_trial_element*nSpace],
3285 vel_grad_test_dS[nDOF_trial_element*nSpace],
3287 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3288 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3289 force_x,force_y,force_z,
3290 force_p_x,force_p_y,force_p_z,
3291 force_v_x,force_v_y,force_v_z,
3294 ck.calculateMapping_elementBoundary(eN,
3300 mesh_trial_trace_ref.data(),
3301 mesh_grad_trial_trace_ref.data(),
3302 boundaryJac_ref.data(),
3308 metricTensorDetSqrt,
3312 ck.calculateMappingVelocity_elementBoundary(eN,
3316 mesh_velocity_dof.data(),
3318 mesh_trial_trace_ref.data(),
3319 xt_ext,yt_ext,zt_ext,
3324 dS = metricTensorDetSqrt*dS_ref[kb];
3326 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3329 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3331 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3332 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3333 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],w_ext);
3335 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3336 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3337 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
3338 for (
int j=0;j<nDOF_trial_element;j++)
3340 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3341 for (
int I=0;I<nSpace;I++)
3342 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3344 ck.calculateGScale(G,normal,h_penalty);
3349 double distance[3], P_normal[3], P_tangent[3];
3350 if(use_ball_as_particle==1)
3358 P_normal[0],P_normal[1],P_normal[2]);
3360 ball_velocity.data(),ball_angular_velocity.data(),
3362 x_ext-dist*P_normal[0],
3363 y_ext-dist*P_normal[1],
3364 z_ext-dist*P_normal[2],
3365 bc_u_ext,bc_v_ext,bc_w_ext);
3370 dist = ebq_global_phi_solid[ebN_kb];
3371 P_normal[0] = ebq_global_grad_phi_solid[ebN_kb*nSpace+0];
3372 P_normal[1] = ebq_global_grad_phi_solid[ebN_kb*nSpace+1];
3373 P_normal[2] = ebq_global_grad_phi_solid[ebN_kb*nSpace+2];
3374 bc_u_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+0];
3375 bc_v_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+1];
3376 bc_w_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+2];
3378 distance[0] = -P_normal[0]*dist;
3379 distance[1] = -P_normal[1]*dist;
3380 distance[2] = -P_normal[2]*dist;
3381 assert(h_penalty>0.0);
3382 if (h_penalty < std::abs(dist))
3383 h_penalty = std::abs(dist);
3387 for(
int i=0;i<3;++i)
3394 double C_adim =
C_sbm*visco/h_penalty;
3395 double beta_adim =
beta_sbm*h_penalty*visco;
3401 const double u_m_uD[3] = {u_ext - bc_u_ext,v_ext - bc_v_ext,w_ext - bc_w_ext};
3402 const double zero_vec[3]={0.,0.,0.};
3406 for (
int i=0;i<nDOF_test_element;i++)
3408 int eN_i = eN*nDOF_test_element+i;
3410 int GlobPos_u = offset_u+stride_u*vel_l2g[eN_i];
3411 int GlobPos_v = offset_v+stride_v*vel_l2g[eN_i];
3412 int GlobPos_w = offset_w+stride_w*vel_l2g[eN_i];
3413 const double phi_i = vel_test_dS[i];
3414 double *grad_phi_i = &vel_grad_test_dS[i*nSpace+0];
3418 globalResidual[GlobPos_u] += C_adim*phi_i*u_m_uD[0];
3419 globalResidual[GlobPos_v] += C_adim*phi_i*u_m_uD[1];
3420 globalResidual[GlobPos_w] += C_adim*phi_i*u_m_uD[2];
3424 globalResidual[GlobPos_u] -= visco * phi_i*res[0];
3425 globalResidual[GlobPos_v] -= visco * phi_i*res[1];
3426 globalResidual[GlobPos_w] -= visco * phi_i*res[2];
3437 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*u_m_uD[0];
3438 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*u_m_uD[1];
3439 globalResidual[GlobPos_w] += C_adim*grad_phi_i_dot_d*u_m_uD[2];
3442 globalResidual[GlobPos_u] += C_adim*grad_phi_i_dot_d*grad_u_d[0];
3443 globalResidual[GlobPos_v] += C_adim*grad_phi_i_dot_d*grad_u_d[1];
3444 globalResidual[GlobPos_w] += C_adim*grad_phi_i_dot_d*grad_u_d[2];
3447 globalResidual[GlobPos_u] += C_adim*phi_i*grad_u_d[0];
3448 globalResidual[GlobPos_v] += C_adim*phi_i*grad_u_d[1];
3449 globalResidual[GlobPos_w] += C_adim*phi_i*grad_u_d[2];
3470 for (
int i=0; i<nDOF_per_element_pressure;++i)
3472 p_ext += p_dof[p_l2g[eN*nDOF_per_element_pressure+i]]*p_trial_trace_ref[ebN_local_kb*nDOF_per_element_pressure+i];
3475 double force_quad_pt[3]={0.0,0.0,0.0},torque_quad_pt[3]={0.0,0.0,0.0},position_vector_to_mass_center[3];
3477 get_stress_in_n(grad_u_ext,grad_v_ext,grad_w_ext,P_normal,p_ext,visco,force_quad_pt);
3479 force_quad_pt[0] *= dS;
3480 force_quad_pt[1] *= dS;
3481 force_quad_pt[2] *= dS;
3482 if(use_ball_as_particle==1)
3495 if(ebN < nElementBoundaries_owned)
3515 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3517 int ebN = exteriorElementBoundariesArray[ebNE],
3518 eN = elementBoundaryElementsArray[ebN*2+0],
3519 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0],
3520 eN_nDOF_trial_element = eN*nDOF_trial_element;
3521 double elementResidual_mesh[nDOF_test_element],
3522 elementResidual_p[nDOF_test_element],
3523 elementResidual_u[nDOF_test_element],
3524 elementResidual_v[nDOF_test_element],
3525 elementResidual_w[nDOF_test_element],
3527 for (
int i=0;i<nDOF_test_element;i++)
3529 elementResidual_mesh[i]=0.0;
3530 elementResidual_p[i]=0.0;
3531 elementResidual_u[i]=0.0;
3532 elementResidual_v[i]=0.0;
3533 elementResidual_w[i]=0.0;
3535 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3537 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3538 ebNE_kb_nSpace = ebNE_kb*nSpace,
3539 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3540 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3550 dmom_u_acc_u_ext=0.0,
3552 dmom_v_acc_v_ext=0.0,
3554 dmom_w_acc_w_ext=0.0,
3555 mass_adv_ext[nSpace],
3556 dmass_adv_u_ext[nSpace],
3557 dmass_adv_v_ext[nSpace],
3558 dmass_adv_w_ext[nSpace],
3559 mom_u_adv_ext[nSpace],
3560 dmom_u_adv_u_ext[nSpace],
3561 dmom_u_adv_v_ext[nSpace],
3562 dmom_u_adv_w_ext[nSpace],
3563 mom_v_adv_ext[nSpace],
3564 dmom_v_adv_u_ext[nSpace],
3565 dmom_v_adv_v_ext[nSpace],
3566 dmom_v_adv_w_ext[nSpace],
3567 mom_w_adv_ext[nSpace],
3568 dmom_w_adv_u_ext[nSpace],
3569 dmom_w_adv_v_ext[nSpace],
3570 dmom_w_adv_w_ext[nSpace],
3571 mom_uu_diff_ten_ext[nSpace],
3572 mom_vv_diff_ten_ext[nSpace],
3573 mom_ww_diff_ten_ext[nSpace],
3574 mom_uv_diff_ten_ext[1],
3575 mom_uw_diff_ten_ext[1],
3576 mom_vu_diff_ten_ext[1],
3577 mom_vw_diff_ten_ext[1],
3578 mom_wu_diff_ten_ext[1],
3579 mom_wv_diff_ten_ext[1],
3580 mom_u_source_ext=0.0,
3581 mom_v_source_ext=0.0,
3582 mom_w_source_ext=0.0,
3584 dmom_u_ham_grad_p_ext[nSpace],
3585 dmom_u_ham_grad_u_ext[nSpace],
3587 dmom_v_ham_grad_p_ext[nSpace],
3588 dmom_v_ham_grad_v_ext[nSpace],
3590 dmom_w_ham_grad_p_ext[nSpace],
3591 dmom_w_ham_grad_w_ext[nSpace],
3592 dmom_u_adv_p_ext[nSpace],
3593 dmom_v_adv_p_ext[nSpace],
3594 dmom_w_adv_p_ext[nSpace],
3596 flux_mom_u_adv_ext=0.0,
3597 flux_mom_v_adv_ext=0.0,
3598 flux_mom_w_adv_ext=0.0,
3599 flux_mom_uu_diff_ext=0.0,
3600 flux_mom_uv_diff_ext=0.0,
3601 flux_mom_uw_diff_ext=0.0,
3602 flux_mom_vu_diff_ext=0.0,
3603 flux_mom_vv_diff_ext=0.0,
3604 flux_mom_vw_diff_ext=0.0,
3605 flux_mom_wu_diff_ext=0.0,
3606 flux_mom_wv_diff_ext=0.0,
3607 flux_mom_ww_diff_ext=0.0,
3612 bc_mom_u_acc_ext=0.0,
3613 bc_dmom_u_acc_u_ext=0.0,
3614 bc_mom_v_acc_ext=0.0,
3615 bc_dmom_v_acc_v_ext=0.0,
3616 bc_mom_w_acc_ext=0.0,
3617 bc_dmom_w_acc_w_ext=0.0,
3618 bc_mass_adv_ext[nSpace],
3619 bc_dmass_adv_u_ext[nSpace],
3620 bc_dmass_adv_v_ext[nSpace],
3621 bc_dmass_adv_w_ext[nSpace],
3622 bc_mom_u_adv_ext[nSpace],
3623 bc_dmom_u_adv_u_ext[nSpace],
3624 bc_dmom_u_adv_v_ext[nSpace],
3625 bc_dmom_u_adv_w_ext[nSpace],
3626 bc_mom_v_adv_ext[nSpace],
3627 bc_dmom_v_adv_u_ext[nSpace],
3628 bc_dmom_v_adv_v_ext[nSpace],
3629 bc_dmom_v_adv_w_ext[nSpace],
3630 bc_mom_w_adv_ext[nSpace],
3631 bc_dmom_w_adv_u_ext[nSpace],
3632 bc_dmom_w_adv_v_ext[nSpace],
3633 bc_dmom_w_adv_w_ext[nSpace],
3634 bc_mom_uu_diff_ten_ext[nSpace],
3635 bc_mom_vv_diff_ten_ext[nSpace],
3636 bc_mom_ww_diff_ten_ext[nSpace],
3637 bc_mom_uv_diff_ten_ext[1],
3638 bc_mom_uw_diff_ten_ext[1],
3639 bc_mom_vu_diff_ten_ext[1],
3640 bc_mom_vw_diff_ten_ext[1],
3641 bc_mom_wu_diff_ten_ext[1],
3642 bc_mom_wv_diff_ten_ext[1],
3643 bc_mom_u_source_ext=0.0,
3644 bc_mom_v_source_ext=0.0,
3645 bc_mom_w_source_ext=0.0,
3646 bc_mom_u_ham_ext=0.0,
3647 bc_dmom_u_ham_grad_p_ext[nSpace],
3648 bc_dmom_u_ham_grad_u_ext[nSpace],
3649 bc_mom_v_ham_ext=0.0,
3650 bc_dmom_v_ham_grad_p_ext[nSpace],
3651 bc_dmom_v_ham_grad_v_ext[nSpace],
3652 bc_mom_w_ham_ext=0.0,
3653 bc_dmom_w_ham_grad_p_ext[nSpace],
3654 bc_dmom_w_ham_grad_w_ext[nSpace],
3655 jac_ext[nSpace*nSpace],
3657 jacInv_ext[nSpace*nSpace],
3658 boundaryJac[nSpace*(nSpace-1)],
3659 metricTensor[(nSpace-1)*(nSpace-1)],
3660 metricTensorDetSqrt,
3661 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
3662 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
3663 vel_grad_test_dS[nDOF_trial_element*nSpace],
3664 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3668 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
3669 force_x,force_y,force_z,force_p_x,force_p_y,force_p_z,force_v_x,force_v_y,force_v_z,r_x,r_y,r_z;
3671 ck.calculateMapping_elementBoundary(eN,
3677 mesh_trial_trace_ref.data(),
3678 mesh_grad_trial_trace_ref.data(),
3679 boundaryJac_ref.data(),
3685 metricTensorDetSqrt,
3689 ck.calculateMappingVelocity_elementBoundary(eN,
3693 mesh_velocity_dof.data(),
3695 mesh_trial_trace_ref.data(),
3696 xt_ext,yt_ext,zt_ext,
3708 dS = metricTensorDetSqrt*dS_ref[kb];
3711 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3712 ck.calculateGScale(G,&ebqe_normal_phi_ext[ebNE_kb_nSpace],h_phi);
3714 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3715 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3716 double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
3721 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3725 p_ext = ebqe_p[ebNE_kb];
3726 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
3727 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
3728 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],w_ext);
3730 for (
int I=0;I<nSpace;I++)
3731 grad_p_ext[I] = ebqe_grad_p[ebNE_kb_nSpace + I];
3732 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3733 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3734 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
3736 for (
int j=0;j<nDOF_trial_element;j++)
3739 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
3740 for (
int I=0;I<nSpace;I++)
3741 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
3743 bc_p_ext = isDOFBoundary_p[ebNE_kb]*ebqe_bc_p_ext[ebNE_kb]+(1-isDOFBoundary_p[ebNE_kb])*p_ext;
3745 bc_u_ext = isDOFBoundary_u[ebNE_kb]*(ebqe_bc_u_ext[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u[ebNE_kb])*u_ext;
3746 bc_v_ext = isDOFBoundary_v[ebNE_kb]*(ebqe_bc_v_ext[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v[ebNE_kb])*v_ext;
3747 bc_w_ext = isDOFBoundary_w[ebNE_kb]*(ebqe_bc_w_ext[ebNE_kb] + MOVING_DOMAIN*zt_ext) + (1-isDOFBoundary_w[ebNE_kb])*w_ext;
3749 porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
3753 double distance_to_omega_solid = 1e10;
3754 if (use_ball_as_particle == 1)
3756 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(), x_ext, y_ext, z_ext, distance_to_omega_solid);
3760 for (
int i = 0; i < nParticles; i++)
3762 double distance_to_i_th_solid = ebq_global_phi_solid[i * nElementBoundaries_global * nQuadraturePoints_elementBoundary + ebNE_kb];
3763 distance_to_omega_solid = (distance_to_i_th_solid < distance_to_omega_solid)?distance_to_i_th_solid:distance_to_omega_solid;
3766 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.);
3775 elementDiameter[eN],
3776 smagorinskyConstant,
3777 turbulenceClosureModel,
3780 ebqe_vf_ext[ebNE_kb],
3781 ebqe_phi_ext[ebNE_kb],
3782 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
3783 distance_to_omega_solid,
3784 ebqe_kappa_phi_ext[ebNE_kb],
3796 ebqe_velocity_star[ebNE_kb_nSpace+0],
3797 ebqe_velocity_star[ebNE_kb_nSpace+1],
3798 ebqe_velocity_star[ebNE_kb_nSpace+2],
3822 mom_uu_diff_ten_ext,
3823 mom_vv_diff_ten_ext,
3824 mom_ww_diff_ten_ext,
3825 mom_uv_diff_ten_ext,
3826 mom_uw_diff_ten_ext,
3827 mom_vu_diff_ten_ext,
3828 mom_vw_diff_ten_ext,
3829 mom_wu_diff_ten_ext,
3830 mom_wv_diff_ten_ext,
3835 dmom_u_ham_grad_p_ext,
3836 dmom_u_ham_grad_u_ext,
3838 dmom_v_ham_grad_p_ext,
3839 dmom_v_ham_grad_v_ext,
3841 dmom_w_ham_grad_p_ext,
3842 dmom_w_ham_grad_w_ext,
3850 MATERIAL_PARAMETERS_AS_FUNCTION,
3851 ebqe_density_as_function[ebNE_kb],
3852 ebqe_dynamic_viscosity_as_function[ebNE_kb],
3855 use_ball_as_particle,
3858 ball_velocity.data(),
3859 ball_angular_velocity.data(),
3860 INT_BY_PARTS_PRESSURE);
3869 elementDiameter[eN],
3870 smagorinskyConstant,
3871 turbulenceClosureModel,
3874 bc_ebqe_vf_ext[ebNE_kb],
3875 bc_ebqe_phi_ext[ebNE_kb],
3876 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
3877 distance_to_omega_solid,
3878 ebqe_kappa_phi_ext[ebNE_kb],
3890 ebqe_velocity_star[ebNE_kb_nSpace+0],
3891 ebqe_velocity_star[ebNE_kb_nSpace+1],
3892 ebqe_velocity_star[ebNE_kb_nSpace+2],
3893 bc_eddy_viscosity_ext,
3895 bc_dmom_u_acc_u_ext,
3897 bc_dmom_v_acc_v_ext,
3899 bc_dmom_w_acc_w_ext,
3905 bc_dmom_u_adv_u_ext,
3906 bc_dmom_u_adv_v_ext,
3907 bc_dmom_u_adv_w_ext,
3909 bc_dmom_v_adv_u_ext,
3910 bc_dmom_v_adv_v_ext,
3911 bc_dmom_v_adv_w_ext,
3913 bc_dmom_w_adv_u_ext,
3914 bc_dmom_w_adv_v_ext,
3915 bc_dmom_w_adv_w_ext,
3916 bc_mom_uu_diff_ten_ext,
3917 bc_mom_vv_diff_ten_ext,
3918 bc_mom_ww_diff_ten_ext,
3919 bc_mom_uv_diff_ten_ext,
3920 bc_mom_uw_diff_ten_ext,
3921 bc_mom_vu_diff_ten_ext,
3922 bc_mom_vw_diff_ten_ext,
3923 bc_mom_wu_diff_ten_ext,
3924 bc_mom_wv_diff_ten_ext,
3925 bc_mom_u_source_ext,
3926 bc_mom_v_source_ext,
3927 bc_mom_w_source_ext,
3929 bc_dmom_u_ham_grad_p_ext,
3930 bc_dmom_u_ham_grad_u_ext,
3932 bc_dmom_v_ham_grad_p_ext,
3933 bc_dmom_v_ham_grad_v_ext,
3935 bc_dmom_w_ham_grad_p_ext,
3936 bc_dmom_w_ham_grad_w_ext,
3944 MATERIAL_PARAMETERS_AS_FUNCTION,
3945 ebqe_density_as_function[ebNE_kb],
3946 ebqe_dynamic_viscosity_as_function[ebNE_kb],
3949 use_ball_as_particle,
3952 ball_velocity.data(),
3953 ball_angular_velocity.data(),
3954 INT_BY_PARTS_PRESSURE);
3957 if (turbulenceClosureModel >= 3)
3959 const double turb_var_grad_0_dummy[3] = {0.,0.,0.};
3960 const double c_mu = 0.09;
3969 ebqe_vf_ext[ebNE_kb],
3970 ebqe_phi_ext[ebNE_kb],
3973 ebqe_turb_var_0[ebNE_kb],
3974 ebqe_turb_var_1[ebNE_kb],
3975 turb_var_grad_0_dummy,
3977 mom_uu_diff_ten_ext,
3978 mom_vv_diff_ten_ext,
3979 mom_ww_diff_ten_ext,
3980 mom_uv_diff_ten_ext,
3981 mom_uw_diff_ten_ext,
3982 mom_vu_diff_ten_ext,
3983 mom_vw_diff_ten_ext,
3984 mom_wu_diff_ten_ext,
3985 mom_wv_diff_ten_ext,
3998 bc_ebqe_vf_ext[ebNE_kb],
3999 bc_ebqe_phi_ext[ebNE_kb],
4002 ebqe_turb_var_0[ebNE_kb],
4003 ebqe_turb_var_1[ebNE_kb],
4004 turb_var_grad_0_dummy,
4005 bc_eddy_viscosity_ext,
4006 bc_mom_uu_diff_ten_ext,
4007 bc_mom_vv_diff_ten_ext,
4008 bc_mom_ww_diff_ten_ext,
4009 bc_mom_uv_diff_ten_ext,
4010 bc_mom_uw_diff_ten_ext,
4011 bc_mom_vu_diff_ten_ext,
4012 bc_mom_vw_diff_ten_ext,
4013 bc_mom_wu_diff_ten_ext,
4014 bc_mom_wv_diff_ten_ext,
4015 bc_mom_u_source_ext,
4016 bc_mom_v_source_ext,
4017 bc_mom_w_source_ext);
4024 mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
4025 mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
4026 mom_u_adv_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*zt_ext;
4027 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
4028 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
4029 dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext;
4031 mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
4032 mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
4033 mom_v_adv_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*zt_ext;
4034 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
4035 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
4036 dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext;
4038 mom_w_adv_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*xt_ext;
4039 mom_w_adv_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*yt_ext;
4040 mom_w_adv_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*zt_ext;
4041 dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext;
4042 dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext;
4043 dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext;
4047 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
4048 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
4049 bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*zt_ext;
4051 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
4052 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
4053 bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*zt_ext;
4055 bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*xt_ext;
4056 bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*yt_ext;
4057 bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*zt_ext;
4062 ck.calculateGScale(G,normal,h_penalty);
4063 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext[ebNE_kb];
4065 isDOFBoundary_u[ebNE_kb],
4066 isDOFBoundary_v[ebNE_kb],
4067 isDOFBoundary_w[ebNE_kb],
4068 isAdvectiveFluxBoundary_p[ebNE_kb],
4069 isAdvectiveFluxBoundary_u[ebNE_kb],
4070 isAdvectiveFluxBoundary_v[ebNE_kb],
4071 isAdvectiveFluxBoundary_w[ebNE_kb],
4072 dmom_u_ham_grad_p_ext[0],
4073 bc_dmom_u_ham_grad_p_ext[0],
4075 porosity_ext*ebqe_rho[ebNE_kb],
4084 ebqe_bc_flux_mass_ext[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]+zt_ext*normal[2]),
4085 ebqe_bc_flux_mom_u_adv_ext[ebNE_kb],
4086 ebqe_bc_flux_mom_v_adv_ext[ebNE_kb],
4087 ebqe_bc_flux_mom_w_adv_ext[ebNE_kb],
4115 &ebqe_velocity_star[ebNE_kb_nSpace],
4116 &ebqe_velocity[ebNE_kb_nSpace]);
4118 for (
int I=0;I<nSpace;I++)
4120 ebqe_grad_u[ebNE_kb_nSpace+I] = grad_u_ext[I];
4121 ebqe_grad_v[ebNE_kb_nSpace+I] = grad_v_ext[I];
4122 ebqe_grad_w[ebNE_kb_nSpace+I] = grad_w_ext[I];
4125 ebqe_phi_ext[ebNE_kb],
4126 sdInfo_u_u_rowptr.data(),
4127 sdInfo_u_u_colind.data(),
4128 isDOFBoundary_u[ebNE_kb],
4129 isDiffusiveFluxBoundary_u[ebNE_kb],
4131 bc_mom_uu_diff_ten_ext,
4133 ebqe_bc_flux_u_diff_ext[ebNE_kb],
4134 mom_uu_diff_ten_ext,
4138 flux_mom_uu_diff_ext);
4140 ebqe_phi_ext[ebNE_kb],
4141 sdInfo_u_v_rowptr.data(),
4142 sdInfo_u_v_colind.data(),
4143 isDOFBoundary_v[ebNE_kb],
4144 isDiffusiveFluxBoundary_v[ebNE_kb],
4146 bc_mom_uv_diff_ten_ext,
4149 mom_uv_diff_ten_ext,
4153 flux_mom_uv_diff_ext);
4155 ebqe_phi_ext[ebNE_kb],
4156 sdInfo_u_w_rowptr.data(),
4157 sdInfo_u_w_colind.data(),
4158 isDOFBoundary_w[ebNE_kb],
4159 isDiffusiveFluxBoundary_w[ebNE_kb],
4161 bc_mom_uw_diff_ten_ext,
4164 mom_uw_diff_ten_ext,
4168 flux_mom_uw_diff_ext);
4170 ebqe_phi_ext[ebNE_kb],
4171 sdInfo_v_u_rowptr.data(),
4172 sdInfo_v_u_colind.data(),
4173 isDOFBoundary_u[ebNE_kb],
4174 isDiffusiveFluxBoundary_u[ebNE_kb],
4176 bc_mom_vu_diff_ten_ext,
4179 mom_vu_diff_ten_ext,
4183 flux_mom_vu_diff_ext);
4185 ebqe_phi_ext[ebNE_kb],
4186 sdInfo_v_v_rowptr.data(),
4187 sdInfo_v_v_colind.data(),
4188 isDOFBoundary_v[ebNE_kb],
4189 isDiffusiveFluxBoundary_v[ebNE_kb],
4191 bc_mom_vv_diff_ten_ext,
4193 ebqe_bc_flux_v_diff_ext[ebNE_kb],
4194 mom_vv_diff_ten_ext,
4198 flux_mom_vv_diff_ext);
4200 ebqe_phi_ext[ebNE_kb],
4201 sdInfo_v_w_rowptr.data(),
4202 sdInfo_v_w_colind.data(),
4203 isDOFBoundary_w[ebNE_kb],
4204 isDiffusiveFluxBoundary_w[ebNE_kb],
4206 bc_mom_vw_diff_ten_ext,
4209 mom_vw_diff_ten_ext,
4213 flux_mom_vw_diff_ext);
4215 ebqe_phi_ext[ebNE_kb],
4216 sdInfo_w_u_rowptr.data(),
4217 sdInfo_w_u_colind.data(),
4218 isDOFBoundary_u[ebNE_kb],
4219 isDiffusiveFluxBoundary_u[ebNE_kb],
4221 bc_mom_wu_diff_ten_ext,
4224 mom_wu_diff_ten_ext,
4228 flux_mom_wu_diff_ext);
4230 ebqe_phi_ext[ebNE_kb],
4231 sdInfo_w_v_rowptr.data(),
4232 sdInfo_w_v_colind.data(),
4233 isDOFBoundary_v[ebNE_kb],
4234 isDiffusiveFluxBoundary_v[ebNE_kb],
4236 bc_mom_wv_diff_ten_ext,
4239 mom_wv_diff_ten_ext,
4243 flux_mom_wv_diff_ext);
4245 ebqe_phi_ext[ebNE_kb],
4246 sdInfo_w_w_rowptr.data(),
4247 sdInfo_w_w_colind.data(),
4248 isDOFBoundary_w[ebNE_kb],
4249 isDiffusiveFluxBoundary_w[ebNE_kb],
4251 bc_mom_ww_diff_ten_ext,
4253 ebqe_bc_flux_w_diff_ext[ebNE_kb],
4254 mom_ww_diff_ten_ext,
4258 flux_mom_ww_diff_ext);
4259 flux[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
4267 if (ebN < nElementBoundaries_owned)
4269 force_v_x = (flux_mom_u_adv_ext + flux_mom_uu_diff_ext + flux_mom_uv_diff_ext + flux_mom_uw_diff_ext)/dmom_u_ham_grad_p_ext[0];
4270 force_v_y = (flux_mom_v_adv_ext + flux_mom_vu_diff_ext + flux_mom_vv_diff_ext + flux_mom_vw_diff_ext)/dmom_u_ham_grad_p_ext[0];
4271 force_v_z = (flux_mom_w_adv_ext + flux_mom_wu_diff_ext + flux_mom_wv_diff_ext + flux_mom_ww_diff_ext)/dmom_u_ham_grad_p_ext[0];
4273 force_p_x = p_ext*normal[0];
4274 force_p_y = p_ext*normal[1];
4275 force_p_z = p_ext*normal[2];
4277 force_x = force_p_x + force_v_x;
4278 force_y = force_p_y + force_v_y;
4279 force_z = force_p_z + force_v_z;
4281 r_x = x_ext - barycenters[3*boundaryFlags[ebN]+0];
4282 r_y = y_ext - barycenters[3*boundaryFlags[ebN]+1];
4283 r_z = z_ext - barycenters[3*boundaryFlags[ebN]+2];
4285 wettedAreas[boundaryFlags[ebN]] += dS*(1.0-ebqe_vf_ext[ebNE_kb]);
4287 netForces_p[3*boundaryFlags[ebN]+0] += force_p_x*dS;
4288 netForces_p[3*boundaryFlags[ebN]+1] += force_p_y*dS;
4289 netForces_p[3*boundaryFlags[ebN]+2] += force_p_z*dS;
4291 netForces_v[3*boundaryFlags[ebN]+0] += force_v_x*dS;
4292 netForces_v[3*boundaryFlags[ebN]+1] += force_v_y*dS;
4293 netForces_v[3*boundaryFlags[ebN]+2] += force_v_z*dS;
4295 netMoments[3*boundaryFlags[ebN]+0] += (r_y*force_z - r_z*force_y)*dS;
4296 netMoments[3*boundaryFlags[ebN]+1] += (r_z*force_x - r_x*force_z)*dS;
4297 netMoments[3*boundaryFlags[ebN]+2] += (r_x*force_y - r_y*force_x)*dS;
4302 for (
int i=0;i<nDOF_test_element;i++)
4304 elementResidual_mesh[i] -=
ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]+zt_ext*normal[2]),p_test_dS[i]);
4305 elementResidual_p[i] +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
4306 elementResidual_p[i] -=
DM*
ck.ExteriorElementBoundaryFlux(MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]+zt_ext*normal[2]),p_test_dS[i]);
4307 globalConservationError +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
4308 elementResidual_u[i] +=
4309 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[0] : 0.) +
4310 ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
4311 ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
4312 ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
4313 ck.ExteriorElementBoundaryFlux(flux_mom_uw_diff_ext,vel_test_dS[i])+
4314 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4315 isDiffusiveFluxBoundary_u[ebNE_kb],
4320 sdInfo_u_u_rowptr.data(),
4321 sdInfo_u_u_colind.data(),
4322 mom_uu_diff_ten_ext,
4323 &vel_grad_test_dS[i*nSpace])+
4324 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4325 isDiffusiveFluxBoundary_u[ebNE_kb],
4330 sdInfo_u_v_rowptr.data(),
4331 sdInfo_u_v_colind.data(),
4332 mom_uv_diff_ten_ext,
4333 &vel_grad_test_dS[i*nSpace])+
4334 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w[ebNE_kb],
4335 isDiffusiveFluxBoundary_u[ebNE_kb],
4340 sdInfo_u_w_rowptr.data(),
4341 sdInfo_u_w_colind.data(),
4342 mom_uw_diff_ten_ext,
4343 &vel_grad_test_dS[i*nSpace]);
4344 elementResidual_v[i] +=
4345 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[1] : 0.) +
4346 ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
4347 ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
4348 ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
4349 ck.ExteriorElementBoundaryFlux(flux_mom_vw_diff_ext,vel_test_dS[i])+
4350 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4351 isDiffusiveFluxBoundary_v[ebNE_kb],
4356 sdInfo_v_u_rowptr.data(),
4357 sdInfo_v_u_colind.data(),
4358 mom_vu_diff_ten_ext,
4359 &vel_grad_test_dS[i*nSpace])+
4360 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4361 isDiffusiveFluxBoundary_v[ebNE_kb],
4366 sdInfo_v_v_rowptr.data(),
4367 sdInfo_v_v_colind.data(),
4368 mom_vv_diff_ten_ext,
4369 &vel_grad_test_dS[i*nSpace])+
4370 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w[ebNE_kb],
4371 isDiffusiveFluxBoundary_v[ebNE_kb],
4376 sdInfo_v_w_rowptr.data(),
4377 sdInfo_v_w_colind.data(),
4378 mom_vw_diff_ten_ext,
4379 &vel_grad_test_dS[i*nSpace]);
4381 elementResidual_w[i] +=
4382 (INT_BY_PARTS_PRESSURE==1 ? p_ext*vel_test_dS[i]*normal[2] : 0.) +
4383 ck.ExteriorElementBoundaryFlux(flux_mom_w_adv_ext,vel_test_dS[i]) +
4384 ck.ExteriorElementBoundaryFlux(flux_mom_wu_diff_ext,vel_test_dS[i])+
4385 ck.ExteriorElementBoundaryFlux(flux_mom_wv_diff_ext,vel_test_dS[i])+
4386 ck.ExteriorElementBoundaryFlux(flux_mom_ww_diff_ext,vel_test_dS[i])+
4387 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u[ebNE_kb],
4388 isDiffusiveFluxBoundary_w[ebNE_kb],
4393 sdInfo_w_u_rowptr.data(),
4394 sdInfo_w_u_colind.data(),
4395 mom_wu_diff_ten_ext,
4396 &vel_grad_test_dS[i*nSpace])+
4397 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v[ebNE_kb],
4398 isDiffusiveFluxBoundary_w[ebNE_kb],
4403 sdInfo_w_v_rowptr.data(),
4404 sdInfo_w_v_colind.data(),
4405 mom_wv_diff_ten_ext,
4406 &vel_grad_test_dS[i*nSpace])+
4407 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w[ebNE_kb],
4408 isDiffusiveFluxBoundary_w[ebNE_kb],
4413 sdInfo_w_w_rowptr.data(),
4414 sdInfo_w_w_colind.data(),
4415 mom_ww_diff_ten_ext,
4416 &vel_grad_test_dS[i*nSpace]);
4422 for (
int i=0;i<nDOF_test_element;i++)
4424 int eN_i = eN*nDOF_test_element+i;
4429 globalResidual[offset_u+stride_u*vel_l2g[eN_i]]+=elementResidual_u[i];
4430 globalResidual[offset_v+stride_v*vel_l2g[eN_i]]+=elementResidual_v[i];
4431 globalResidual[offset_w+stride_w*vel_l2g[eN_i]]+=elementResidual_w[i];
4446 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
4447 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
4448 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4449 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
4450 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
4451 double PSTAB = args.
scalar<
double>(
"PSTAB");
4452 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4453 xt::pyarray<double>& x_ref = args.
array<
double>(
"x_ref");
4454 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
4455 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
4456 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
4457 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
4458 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
4459 xt::pyarray<double>& q_p = args.
array<
double>(
"q_p");
4460 xt::pyarray<double>& q_grad_p = args.
array<
double>(
"q_grad_p");
4461 xt::pyarray<double>& ebqe_p = args.
array<
double>(
"ebqe_p");
4462 xt::pyarray<double>& ebqe_grad_p = args.
array<
double>(
"ebqe_grad_p");
4463 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
4464 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
4465 xt::pyarray<double>& vel_hess_trial_ref = args.
array<
double>(
"vel_hess_trial_ref");
4466 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
4467 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
4468 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4469 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4470 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
4471 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
4472 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
4473 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
4474 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
4475 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4476 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
4477 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
4478 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
4479 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4480 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4481 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
4482 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
4483 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
4484 double hFactor = args.
scalar<
double>(
"hFactor");
4485 int nElements_global = args.
scalar<
int>(
"nElements_global");
4486 int nElements_owned = args.
scalar<
int>(
"nElements_owned");
4487 int nElementBoundaries_global = args.
scalar<
int>(
"nElementBoundaries_global");
4488 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
4489 int nNodes_owned = args.
scalar<
int>(
"nNodes_owned");
4490 double useRBLES = args.
scalar<
double>(
"useRBLES");
4491 double useMetrics = args.
scalar<
double>(
"useMetrics");
4492 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
4493 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
4494 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
4495 double sigma = args.
scalar<
double>(
"sigma");
4500 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
4501 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
4502 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
4503 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
4504 double C_dg = args.
scalar<
double>(
"C_dg");
4505 double C_b = args.
scalar<
double>(
"C_b");
4506 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
4507 const xt::pyarray<double>& ebq_global_phi_solid = args.
array<
double>(
"ebq_global_phi_solid");
4508 const xt::pyarray<double>& ebq_global_grad_phi_solid = args.
array<
double>(
"ebq_global_grad_phi_solid");
4509 const xt::pyarray<double>& ebq_particle_velocity_solid = args.
array<
double>(
"ebq_particle_velocity_solid");
4510 xt::pyarray<double>& phi_solid_nodes = args.
array<
double>(
"phi_solid_nodes");
4511 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
4512 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
4513 const xt::pyarray<double>& q_velocityStar_solid = args.
array<
double>(
"q_velocityStar_solid");
4514 const xt::pyarray<double>& q_vos = args.
array<
double>(
"q_vos");
4515 const xt::pyarray<double>& q_dvos_dt = args.
array<
double>(
"q_dvos_dt");
4516 const xt::pyarray<double>& q_grad_vos = args.
array<
double>(
"q_grad_vos");
4517 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
4518 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
4519 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
4520 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
4521 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
4522 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
4523 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
4524 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4525 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
4526 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4527 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4528 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4529 xt::pyarray<double>& g = args.
array<
double>(
"g");
4530 const double useVF = args.
scalar<
double>(
"useVF");
4531 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
4532 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
4533 xt::pyarray<double>& phi_dof = args.
array<
double>(
"phi_dof");
4534 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
4535 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
4536 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
4537 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
4538 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
4539 xt::pyarray<double>& q_dV = args.
array<
double>(
"q_dV");
4540 xt::pyarray<double>& q_dV_last = args.
array<
double>(
"q_dV_last");
4541 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
4542 xt::pyarray<double>& ebqe_velocity_star = args.
array<
double>(
"ebqe_velocity_star");
4543 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
4544 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
4545 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
4546 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
4547 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
4548 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
4549 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
4550 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
4551 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
4552 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
4553 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
4554 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
4555 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
4556 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
4557 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
4558 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
4559 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
4560 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
4561 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
4562 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
4563 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
4564 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
4565 xt::pyarray<int>& csrRowIndeces_p_p = args.
array<
int>(
"csrRowIndeces_p_p");
4566 xt::pyarray<int>& csrColumnOffsets_p_p = args.
array<
int>(
"csrColumnOffsets_p_p");
4567 xt::pyarray<int>& csrRowIndeces_p_u = args.
array<
int>(
"csrRowIndeces_p_u");
4568 xt::pyarray<int>& csrColumnOffsets_p_u = args.
array<
int>(
"csrColumnOffsets_p_u");
4569 xt::pyarray<int>& csrRowIndeces_p_v = args.
array<
int>(
"csrRowIndeces_p_v");
4570 xt::pyarray<int>& csrColumnOffsets_p_v = args.
array<
int>(
"csrColumnOffsets_p_v");
4571 xt::pyarray<int>& csrRowIndeces_p_w = args.
array<
int>(
"csrRowIndeces_p_w");
4572 xt::pyarray<int>& csrColumnOffsets_p_w = args.
array<
int>(
"csrColumnOffsets_p_w");
4573 xt::pyarray<int>& csrRowIndeces_u_p = args.
array<
int>(
"csrRowIndeces_u_p");
4574 xt::pyarray<int>& csrColumnOffsets_u_p = args.
array<
int>(
"csrColumnOffsets_u_p");
4575 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
4576 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
4577 xt::pyarray<int>& csrRowIndeces_u_v = args.
array<
int>(
"csrRowIndeces_u_v");
4578 xt::pyarray<int>& csrColumnOffsets_u_v = args.
array<
int>(
"csrColumnOffsets_u_v");
4579 xt::pyarray<int>& csrRowIndeces_u_w = args.
array<
int>(
"csrRowIndeces_u_w");
4580 xt::pyarray<int>& csrColumnOffsets_u_w = args.
array<
int>(
"csrColumnOffsets_u_w");
4581 xt::pyarray<int>& csrRowIndeces_v_p = args.
array<
int>(
"csrRowIndeces_v_p");
4582 xt::pyarray<int>& csrColumnOffsets_v_p = args.
array<
int>(
"csrColumnOffsets_v_p");
4583 xt::pyarray<int>& csrRowIndeces_v_u = args.
array<
int>(
"csrRowIndeces_v_u");
4584 xt::pyarray<int>& csrColumnOffsets_v_u = args.
array<
int>(
"csrColumnOffsets_v_u");
4585 xt::pyarray<int>& csrRowIndeces_v_v = args.
array<
int>(
"csrRowIndeces_v_v");
4586 xt::pyarray<int>& csrColumnOffsets_v_v = args.
array<
int>(
"csrColumnOffsets_v_v");
4587 xt::pyarray<int>& csrRowIndeces_v_w = args.
array<
int>(
"csrRowIndeces_v_w");
4588 xt::pyarray<int>& csrColumnOffsets_v_w = args.
array<
int>(
"csrColumnOffsets_v_w");
4589 xt::pyarray<int>& csrRowIndeces_w_p = args.
array<
int>(
"csrRowIndeces_w_p");
4590 xt::pyarray<int>& csrColumnOffsets_w_p = args.
array<
int>(
"csrColumnOffsets_w_p");
4591 xt::pyarray<int>& csrRowIndeces_w_u = args.
array<
int>(
"csrRowIndeces_w_u");
4592 xt::pyarray<int>& csrColumnOffsets_w_u = args.
array<
int>(
"csrColumnOffsets_w_u");
4593 xt::pyarray<int>& csrRowIndeces_w_v = args.
array<
int>(
"csrRowIndeces_w_v");
4594 xt::pyarray<int>& csrColumnOffsets_w_v = args.
array<
int>(
"csrColumnOffsets_w_v");
4595 xt::pyarray<int>& csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
4596 xt::pyarray<int>& csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
4597 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
4598 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4599 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4600 xt::pyarray<int>& elementBoundariesArray = args.
array<
int>(
"elementBoundariesArray");
4601 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4602 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4603 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
4604 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
4605 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
4606 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
4607 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
4608 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
4609 const xt::pyarray<double>& ebqe_vos_ext = args.
array<
double>(
"ebqe_vos_ext");
4610 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
4611 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
4612 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
4613 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
4614 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
4615 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
4616 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
4617 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
4618 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
4619 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
4620 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
4621 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
4622 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
4623 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
4624 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
4625 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
4626 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
4627 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
4628 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
4629 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
4630 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
4631 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
4632 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
4633 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
4634 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
4635 xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.
array<
int>(
"csrColumnOffsets_eb_p_p");
4636 xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.
array<
int>(
"csrColumnOffsets_eb_p_u");
4637 xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.
array<
int>(
"csrColumnOffsets_eb_p_v");
4638 xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.
array<
int>(
"csrColumnOffsets_eb_p_w");
4639 xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.
array<
int>(
"csrColumnOffsets_eb_u_p");
4640 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
4641 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.
array<
int>(
"csrColumnOffsets_eb_u_v");
4642 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.
array<
int>(
"csrColumnOffsets_eb_u_w");
4643 xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.
array<
int>(
"csrColumnOffsets_eb_v_p");
4644 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.
array<
int>(
"csrColumnOffsets_eb_v_u");
4645 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.
array<
int>(
"csrColumnOffsets_eb_v_v");
4646 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.
array<
int>(
"csrColumnOffsets_eb_v_w");
4647 xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.
array<
int>(
"csrColumnOffsets_eb_w_p");
4648 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.
array<
int>(
"csrColumnOffsets_eb_w_u");
4649 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.
array<
int>(
"csrColumnOffsets_eb_w_v");
4650 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.
array<
int>(
"csrColumnOffsets_eb_w_w");
4651 xt::pyarray<int>& elementFlags = args.
array<
int>(
"elementFlags");
4652 int nParticles = args.
scalar<
int>(
"nParticles");
4653 double particle_epsFact = args.
scalar<
double>(
"particle_epsFact");
4654 double particle_alpha = args.
scalar<
double>(
"particle_alpha");
4655 double particle_beta = args.
scalar<
double>(
"particle_beta");
4656 double particle_penalty_constant = args.
scalar<
double>(
"particle_penalty_constant");
4657 xt::pyarray<double>& particle_signed_distances = args.
array<
double>(
"particle_signed_distances");
4658 xt::pyarray<double>& particle_signed_distance_normals = args.
array<
double>(
"particle_signed_distance_normals");
4659 xt::pyarray<double>& particle_velocities = args.
array<
double>(
"particle_velocities");
4660 xt::pyarray<double>& particle_centroids = args.
array<
double>(
"particle_centroids");
4661 double particle_nitsche = args.
scalar<
double>(
"particle_nitsche");
4662 int use_ball_as_particle = args.
scalar<
int>(
"use_ball_as_particle");
4663 xt::pyarray<double>& ball_center = args.
array<
double>(
"ball_center");
4664 xt::pyarray<double>& ball_radius = args.
array<
double>(
"ball_radius");
4665 xt::pyarray<double>& ball_velocity = args.
array<
double>(
"ball_velocity");
4666 xt::pyarray<double>& ball_angular_velocity = args.
array<
double>(
"ball_angular_velocity");
4667 int USE_SUPG = args.
scalar<
int>(
"USE_SUPG");
4668 int KILL_PRESSURE_TERM = args.
scalar<
int>(
"KILL_PRESSURE_TERM");
4669 double dt = args.
scalar<
double>(
"dt");
4670 int MATERIAL_PARAMETERS_AS_FUNCTION = args.
scalar<
int>(
"MATERIAL_PARAMETERS_AS_FUNCTION");
4671 xt::pyarray<double>& density_as_function = args.
array<
double>(
"density_as_function");
4672 xt::pyarray<double>& dynamic_viscosity_as_function = args.
array<
double>(
"dynamic_viscosity_as_function");
4673 xt::pyarray<double>& ebqe_density_as_function = args.
array<
double>(
"ebqe_density_as_function");
4674 xt::pyarray<double>& ebqe_dynamic_viscosity_as_function = args.
array<
double>(
"ebqe_dynamic_viscosity_as_function");
4675 int USE_SBM = args.
scalar<
int>(
"USE_SBM");
4676 int ARTIFICIAL_VISCOSITY = args.
scalar<
int>(
"ARTIFICIAL_VISCOSITY");
4677 xt::pyarray<double>& uStar_dMatrix = args.
array<
double>(
"uStar_dMatrix");
4678 xt::pyarray<double>& vStar_dMatrix = args.
array<
double>(
"vStar_dMatrix");
4679 xt::pyarray<double>& wStar_dMatrix = args.
array<
double>(
"wStar_dMatrix");
4680 int numDOFs_1D = args.
scalar<
int>(
"numDOFs_1D");
4681 int offset_u = args.
scalar<
int>(
"offset_u");
4682 int offset_v = args.
scalar<
int>(
"offset_v");
4683 int offset_w = args.
scalar<
int>(
"offset_w");
4684 int stride_u = args.
scalar<
int>(
"stride_u");
4685 int stride_v = args.
scalar<
int>(
"stride_v");
4686 int stride_w = args.
scalar<
int>(
"stride_w");
4687 xt::pyarray<int>& rowptr_1D = args.
array<
int>(
"rowptr_1D");
4688 xt::pyarray<int>& colind_1D = args.
array<
int>(
"colind_1D");
4689 xt::pyarray<int>& rowptr = args.
array<
int>(
"rowptr");
4690 xt::pyarray<int>& colind = args.
array<
int>(
"colind");
4691 int INT_BY_PARTS_PRESSURE = args.
scalar<
int>(
"INT_BY_PARTS_PRESSURE");
4695 std::valarray<double> particle_surfaceArea(nParticles), particle_netForces(nParticles * 3), particle_netMoments(nParticles * 3);
4696 const int nQuadraturePoints_global(nElements_global*nQuadraturePoints_element);
4697 for(
int eN=0;eN<nElements_global;eN++)
4699 double eps_rho,eps_mu;
4700 double element_active=1.0;
4702 double elementJacobian_p_p[nDOF_test_element][nDOF_trial_element],
4703 elementJacobian_p_u[nDOF_test_element][nDOF_trial_element],
4704 elementJacobian_p_v[nDOF_test_element][nDOF_trial_element],
4705 elementJacobian_p_w[nDOF_test_element][nDOF_trial_element],
4706 elementJacobian_u_p[nDOF_test_element][nDOF_trial_element],
4707 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
4708 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
4709 elementJacobian_u_w[nDOF_test_element][nDOF_trial_element],
4710 elementJacobian_v_p[nDOF_test_element][nDOF_trial_element],
4711 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
4712 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element],
4713 elementJacobian_v_w[nDOF_test_element][nDOF_trial_element],
4714 elementJacobian_w_p[nDOF_test_element][nDOF_trial_element],
4715 elementJacobian_w_u[nDOF_test_element][nDOF_trial_element],
4716 elementJacobian_w_v[nDOF_test_element][nDOF_trial_element],
4717 elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
4718 for (
int i=0;i<nDOF_test_element;i++)
4719 for (
int j=0;j<nDOF_trial_element;j++)
4721 elementJacobian_p_p[i][j]=0.0;
4722 elementJacobian_p_u[i][j]=0.0;
4723 elementJacobian_p_v[i][j]=0.0;
4724 elementJacobian_p_w[i][j]=0.0;
4725 elementJacobian_u_p[i][j]=0.0;
4726 elementJacobian_u_u[i][j]=0.0;
4727 elementJacobian_u_v[i][j]=0.0;
4728 elementJacobian_u_w[i][j]=0.0;
4729 elementJacobian_v_p[i][j]=0.0;
4730 elementJacobian_v_u[i][j]=0.0;
4731 elementJacobian_v_v[i][j]=0.0;
4732 elementJacobian_v_w[i][j]=0.0;
4733 elementJacobian_w_p[i][j]=0.0;
4734 elementJacobian_w_u[i][j]=0.0;
4735 elementJacobian_w_v[i][j]=0.0;
4736 elementJacobian_w_w[i][j]=0.0;
4743 double _distance[nDOF_mesh_trial_element]={0.0};
4745 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4747 if(use_ball_as_particle==1)
4750 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+0],
4751 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+1],
4752 mesh_dof[3*mesh_l2g[eN*nDOF_mesh_trial_element+I]+2],
4757 _distance[I] = phi_solid_nodes[mesh_l2g[eN*nDOF_mesh_trial_element+I]];
4759 if ( _distance[I] >= 0)
4762 if (pos_counter == 3)
4766 for (
int I=0;I<nDOF_mesh_trial_element;I++)
4768 if (_distance[I] < 0)
4771 assert(opp_node >=0);
4772 assert(opp_node <nDOF_mesh_trial_element);
4774 else if (pos_counter == 4)
4783 double element_phi[nDOF_mesh_trial_element], element_phi_s[nDOF_mesh_trial_element];
4784 for (
int j=0;j<nDOF_mesh_trial_element;j++)
4786 int eN_j = eN*nDOF_mesh_trial_element+j;
4787 element_phi[j] = phi_dof[p_l2g[eN_j]];
4788 element_phi_s[j] = phi_solid_nodes[p_l2g[eN_j]];
4790 double element_nodes[nDOF_mesh_trial_element*3];
4791 for (
int i=0;i<nDOF_mesh_trial_element;i++)
4793 int eN_i=eN*nDOF_mesh_trial_element+i;
4794 for(
int I=0;I<3;I++)
4795 element_nodes[i*3 + I] = mesh_dof[mesh_l2g[eN_i]*3 + I];
4797 gf_s.
calculate(element_phi_s, element_nodes, x_ref.data(),
false);
4798 gf.
calculate(element_phi, element_nodes, x_ref.data(),
false);
4799 for (
int k=0;k<nQuadraturePoints_element;k++)
4803 int eN_k = eN*nQuadraturePoints_element+k,
4804 eN_k_nSpace = eN_k*nSpace,
4805 eN_nDOF_trial_element = eN*nDOF_trial_element;
4808 double p=0.0,
u=0.0,
v=0.0,
w=0.0,
4809 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
4818 dmass_adv_u[nSpace],
4819 dmass_adv_v[nSpace],
4820 dmass_adv_w[nSpace],
4822 dmom_u_adv_u[nSpace],
4823 dmom_u_adv_v[nSpace],
4824 dmom_u_adv_w[nSpace],
4826 dmom_v_adv_u[nSpace],
4827 dmom_v_adv_v[nSpace],
4828 dmom_v_adv_w[nSpace],
4830 dmom_w_adv_u[nSpace],
4831 dmom_w_adv_v[nSpace],
4832 dmom_w_adv_w[nSpace],
4833 mom_uu_diff_ten[nSpace],
4834 mom_vv_diff_ten[nSpace],
4835 mom_ww_diff_ten[nSpace],
4846 dmom_u_ham_grad_p[nSpace],
4847 dmom_u_ham_grad_u[nSpace],
4849 dmom_v_ham_grad_p[nSpace],
4850 dmom_v_ham_grad_v[nSpace],
4852 dmom_w_ham_grad_p[nSpace],
4853 dmom_w_ham_grad_w[nSpace],
4864 dpdeResidual_p_u[nDOF_trial_element],dpdeResidual_p_v[nDOF_trial_element],dpdeResidual_p_w[nDOF_trial_element],
4865 dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_trial_element],
4866 dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_trial_element],
4867 dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_trial_element],
4868 Lstar_u_p[nDOF_test_element],
4869 Lstar_v_p[nDOF_test_element],
4870 Lstar_w_p[nDOF_test_element],
4871 Lstar_u_u[nDOF_test_element],
4872 Lstar_v_v[nDOF_test_element],
4873 Lstar_w_w[nDOF_test_element],
4874 Lstar_p_u[nDOF_test_element],
4875 Lstar_p_v[nDOF_test_element],
4876 Lstar_p_w[nDOF_test_element],
4881 dsubgridError_p_u[nDOF_trial_element],
4882 dsubgridError_p_v[nDOF_trial_element],
4883 dsubgridError_p_w[nDOF_trial_element],
4884 dsubgridError_u_p[nDOF_trial_element],
4885 dsubgridError_u_u[nDOF_trial_element],
4886 dsubgridError_v_p[nDOF_trial_element],
4887 dsubgridError_v_v[nDOF_trial_element],
4888 dsubgridError_w_p[nDOF_trial_element],
4889 dsubgridError_w_w[nDOF_trial_element],
4890 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
4891 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
4894 jacInv[nSpace*nSpace],
4895 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
4896 vel_hess_trial[nDOF_trial_element*
nSpace2],
4898 p_test_dV[nDOF_test_element],vel_test_dV[nDOF_test_element],
4899 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
4904 dmom_u_source[nSpace],
4905 dmom_v_source[nSpace],
4906 dmom_w_source[nSpace],
4909 G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
4911 ck.calculateMapping_element(eN,
4915 mesh_trial_ref.data(),
4916 mesh_grad_trial_ref.data(),
4921 ck.calculateH_element(eN,
4923 nodeDiametersArray.data(),
4925 mesh_trial_ref.data(),
4927 ck.calculateMappingVelocity_element(eN,
4929 mesh_velocity_dof.data(),
4931 mesh_trial_ref.data(),
4936 dV = fabs(jacDet)*dV_ref[k];
4937 ck.calculateG(jacInv,G,G_dd_G,tr_G);
4940 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
4941 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
4942 const double particle_eps = particle_epsFact*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
4946 ck.gradTrialFromRef(&vel_grad_trial_ref[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
4947 ck.hessTrialFromRef(&vel_hess_trial_ref[k*nDOF_trial_element*
nSpace2],jacInv,vel_hess_trial);
4951 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
u);
4952 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
v);
4953 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_ref[k*nDOF_trial_element],
w);
4956 for (
int I=0;I<nSpace;I++)
4957 grad_p[I] = q_grad_p[eN_k_nSpace+I];
4958 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_u);
4959 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_v);
4960 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial,grad_w);
4961 ck.hessFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_u);
4962 ck.hessFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_v);
4963 ck.hessFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_hess_trial,hess_w);
4965 for (
int j=0;j<nDOF_trial_element;j++)
4968 vel_test_dV[j] = vel_test_ref[k*nDOF_trial_element+j]*dV;
4969 for (
int I=0;I<nSpace;I++)
4972 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
4976 double div_mesh_velocity=0.0;
4977 int NDOF_MESH_TRIAL_ELEMENT=4;
4978 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
4980 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
4981 div_mesh_velocity +=
4982 mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*vel_grad_trial[j*3+0] +
4983 mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*vel_grad_trial[j*3+1] +
4984 mesh_velocity_dof[mesh_l2g[eN_j]*3+2]*vel_grad_trial[j*3+2];
4986 div_mesh_velocity =
DM3*div_mesh_velocity + (1.0-
DM3)*alphaBDF*(dV-q_dV_last[eN_k])/dV;
4989 porosity = 1.0 - q_vos[eN_k];
4994 double distance_to_omega_solid = phi_solid[eN_k];
4995 double eddy_viscosity(0.),rhoSave,nuSave;
5004 elementDiameter[eN],
5005 smagorinskyConstant,
5006 turbulenceClosureModel,
5011 &normal_phi[eN_k_nSpace],
5012 distance_to_omega_solid,
5025 q_velocity_sge[eN_k_nSpace+0],
5026 q_velocity_sge[eN_k_nSpace+1],
5027 q_velocity_sge[eN_k_nSpace+2],
5079 MATERIAL_PARAMETERS_AS_FUNCTION,
5080 density_as_function[eN_k],
5081 dynamic_viscosity_as_function[eN_k],
5084 use_ball_as_particle,
5087 ball_velocity.data(),
5088 ball_angular_velocity.data(),
5089 INT_BY_PARTS_PRESSURE);
5091 mass_source = q_mass_source[eN_k];
5092 for (
int I=0;I<nSpace;I++)
5094 dmom_u_source[I] = 0.0;
5095 dmom_v_source[I] = 0.0;
5096 dmom_w_source[I] = 0.0;
5117 q_velocity_sge[eN_k_nSpace+0],
5118 q_velocity_sge[eN_k_nSpace+1],
5119 q_velocity_sge[eN_k_nSpace+2],
5120 eps_solid[elementFlags[eN]],
5122 q_velocity_solid[eN_k_nSpace+0],
5123 q_velocity_solid[eN_k_nSpace+1],
5124 q_velocity_solid[eN_k_nSpace+2],
5125 q_velocityStar_solid[eN_k_nSpace+0],
5126 q_velocityStar_solid[eN_k_nSpace+1],
5127 q_velocityStar_solid[eN_k_nSpace+2],
5134 q_grad_vos[eN_k_nSpace+0],
5135 q_grad_vos[eN_k_nSpace+1],
5136 q_grad_vos[eN_k_nSpace+2]);
5137 double C_particles=0.0;
5138 if(nParticles > 0 && USE_SBM==0)
5143 nQuadraturePoints_global,
5144 &particle_signed_distances[eN_k],
5145 &particle_signed_distance_normals[eN_k_nSpace],
5146 &particle_velocities[eN_k_nSpace],
5147 particle_centroids.data(),
5148 use_ball_as_particle,
5151 ball_velocity.data(),
5152 ball_angular_velocity.data(),
5154 particle_penalty_constant/h_phi,
5155 particle_alpha/h_phi,
5156 particle_beta/h_phi,
5173 q_velocity_sge[eN_k_nSpace+0],
5174 q_velocity_sge[eN_k_nSpace+1],
5175 q_velocity_sge[eN_k_nSpace+2],
5198 &particle_netForces[0],
5199 &particle_netMoments[0],
5200 &particle_surfaceArea[0]);
5202 if (turbulenceClosureModel >= 3)
5204 const double c_mu = 0.09;
5219 &q_turb_var_grad_0[eN_k_nSpace],
5239 mom_u_adv[0] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*
xt;
5240 mom_u_adv[1] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*yt;
5241 mom_u_adv[2] -= MOVING_DOMAIN*dmom_u_acc_u*mom_u_acc*zt;
5242 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
5243 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
5244 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
5246 mom_v_adv[0] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*
xt;
5247 mom_v_adv[1] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*yt;
5248 mom_v_adv[2] -= MOVING_DOMAIN*dmom_v_acc_v*mom_v_acc*zt;
5249 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
5250 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
5251 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
5253 mom_w_adv[0] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*
xt;
5254 mom_w_adv[1] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*yt;
5255 mom_w_adv[2] -= MOVING_DOMAIN*dmom_w_acc_w*mom_w_acc*zt;
5256 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
5257 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
5258 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
5263 q_mom_u_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5269 q_mom_v_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5275 q_mom_w_acc_beta_bdf[eN_k]*q_dV_last[eN_k]/dV,
5281 mom_u_acc_t *= dmom_u_acc_u;
5282 mom_v_acc_t *= dmom_v_acc_v;
5283 mom_w_acc_t *= dmom_w_acc_w;
5288 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
5289 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
5290 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
5295 ck.Mass_strong(-q_dvos_dt[eN_k]) +
5296 ck.Advection_strong(dmass_adv_u,grad_u) +
5297 ck.Advection_strong(dmass_adv_v,grad_v) +
5298 ck.Advection_strong(dmass_adv_w,grad_w) +
5299 DM2*MOVING_DOMAIN*
ck.Reaction_strong(alphaBDF*(dV-q_dV_last[eN_k])/dV - div_mesh_velocity) +
5301 ck.Reaction_strong(mass_source);
5305 ck.Mass_strong(mom_u_acc_t) +
5306 ck.Advection_strong(dmom_adv_sge,grad_u) +
5307 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
5308 ck.Reaction_strong(mom_u_source) -
5309 ck.Reaction_strong(
u*div_mesh_velocity);
5312 ck.Mass_strong(mom_v_acc_t) +
5313 ck.Advection_strong(dmom_adv_sge,grad_v) +
5314 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
5315 ck.Reaction_strong(mom_v_source) -
5316 ck.Reaction_strong(
v*div_mesh_velocity);
5318 pdeResidual_w =
ck.Mass_strong(mom_w_acc_t) +
5319 ck.Advection_strong(dmom_adv_sge,grad_w) +
5320 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
5321 ck.Reaction_strong(mom_w_source) -
5322 ck.Reaction_strong(
w*div_mesh_velocity);
5325 for (
int j=0;j<nDOF_trial_element;j++)
5327 int j_nSpace = j*nSpace;
5328 dpdeResidual_p_u[j]=
ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial[j_nSpace]);
5329 dpdeResidual_p_v[j]=
ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial[j_nSpace]);
5330 dpdeResidual_p_w[j]=
ck.AdvectionJacobian_strong(dmass_adv_w,&vel_grad_trial[j_nSpace]);
5332 dpdeResidual_u_p[j]=
ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace]);
5333 dpdeResidual_u_u[j]=
ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5334 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5335 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5337 dpdeResidual_v_p[j]=
ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace]);
5338 dpdeResidual_v_v[j]=
ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5339 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5340 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5342 dpdeResidual_w_p[j]=
ck.HamiltonianJacobian_strong(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace]);
5343 dpdeResidual_w_w[j]=
ck.MassJacobian_strong(dmom_w_acc_w_t,vel_trial_ref[k*nDOF_trial_element+j]) +
5344 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]) -
5345 ck.ReactionJacobian_strong(div_mesh_velocity,vel_trial_ref[k*nDOF_trial_element+j]);
5348 dpdeResidual_u_u[j]+=
ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial_ref[k*nDOF_trial_element+j]);
5349 dpdeResidual_v_v[j]+=
ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial_ref[k*nDOF_trial_element+j]);
5350 dpdeResidual_w_w[j]+=
ck.ReactionJacobian_strong(dmom_w_source[2],vel_trial_ref[k*nDOF_trial_element+j]);
5355 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
5357 elementDiameter[eN],
5362 dmom_u_ham_grad_p[0],
5372 dmom_u_ham_grad_p[0],
5378 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
5379 tau_p = KILL_PRESSURE_TERM == 1 ? 0. : PSTAB*(useMetrics*tau_p1+(1.0-useMetrics)*tau_p0);
5412 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
5413 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
5414 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
5417 for (
int i=0;i<nDOF_test_element;i++)
5419 int i_nSpace = i*nSpace;
5420 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
5421 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
5422 Lstar_w_p[i]=
ck.Advection_adjoint(dmass_adv_w,&p_grad_test_dV[i_nSpace]);
5423 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5424 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5425 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
5426 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5427 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5428 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
5430 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
5431 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
5432 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
5436 dmom_u_adv_u[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5437 dmom_u_adv_u[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5438 dmom_u_adv_u[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
5440 dmom_v_adv_v[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5441 dmom_v_adv_v[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5442 dmom_v_adv_v[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
5444 dmom_w_adv_w[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
5445 dmom_w_adv_w[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
5446 dmom_w_adv_w[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
5449 double unit_normal[nSpace];
5450 double norm_grad_phi = 0.;
5451 for (
int I=0;I<nSpace;I++)
5452 norm_grad_phi += normal_phi[eN_k_nSpace+I]*normal_phi[eN_k_nSpace+I];
5453 norm_grad_phi = std::sqrt(norm_grad_phi) + 1E-10;
5454 for (
int I=0;I<nSpace;I++)
5455 unit_normal[I] = normal_phi[eN_k_nSpace+I]/norm_grad_phi;
5456 double delta =
gf.
D(eps_mu,
phi[eN_k]);
5457 double vel_tgrad_test_i[nSpace], vel_tgrad_test_j[nSpace];
5461 for(
int i=0;i<nDOF_test_element;i++)
5463 int i_nSpace = i*nSpace;
5465 &vel_grad_trial[i_nSpace],
5467 for(
int j=0;j<nDOF_trial_element;j++)
5469 int j_nSpace = j*nSpace;
5471 &vel_grad_trial[j_nSpace],
5487 elementJacobian_u_u[i][j] +=
5488 ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5489 ck.HamiltonianJacobian_weak(dmom_u_ham_grad_u,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
5490 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5491 ck.SimpleDiffusionJacobian_weak(sdInfo_u_u_rowptr.data(),sdInfo_u_u_colind.data(),mom_uu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5493 ck.ReactionJacobian_weak(dmom_u_source[0],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5496 USE_SUPG*
ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
5497 ck.NumericalDiffusionJacobian(q_numDiff_u_last[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5499 ck.NumericalDiffusion(dt*delta*sigma*dV,
5503 elementJacobian_u_v[i][j] +=
5504 ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5505 ck.SimpleDiffusionJacobian_weak(sdInfo_u_v_rowptr.data(),sdInfo_u_v_colind.data(),mom_uv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5507 ck.ReactionJacobian_weak(dmom_u_source[1],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5510 elementJacobian_u_w[i][j] +=
5511 ck.AdvectionJacobian_weak(dmom_u_adv_w,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5512 ck.SimpleDiffusionJacobian_weak(sdInfo_u_w_rowptr.data(),sdInfo_u_w_colind.data(),mom_uw_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5514 ck.ReactionJacobian_weak(dmom_u_source[2],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5521 elementJacobian_v_u[i][j] +=
5522 ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5523 ck.SimpleDiffusionJacobian_weak(sdInfo_v_u_rowptr.data(),sdInfo_v_u_colind.data(),mom_vu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5525 ck.ReactionJacobian_weak(dmom_v_source[0],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5528 elementJacobian_v_v[i][j] +=
5529 ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5530 ck.HamiltonianJacobian_weak(dmom_v_ham_grad_v,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
5531 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5532 ck.SimpleDiffusionJacobian_weak(sdInfo_v_v_rowptr.data(),sdInfo_v_v_colind.data(),mom_vv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5534 ck.ReactionJacobian_weak(dmom_v_source[1],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5537 USE_SUPG*
ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
5538 ck.NumericalDiffusionJacobian(q_numDiff_v_last[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5540 ck.NumericalDiffusion(dt*delta*sigma*dV,
5544 elementJacobian_v_w[i][j] +=
5545 ck.AdvectionJacobian_weak(dmom_v_adv_w,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5546 ck.SimpleDiffusionJacobian_weak(sdInfo_v_w_rowptr.data(),sdInfo_v_w_colind.data(),mom_vw_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5548 ck.ReactionJacobian_weak(dmom_v_source[2],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5554 elementJacobian_w_u[i][j] +=
5555 ck.AdvectionJacobian_weak(dmom_w_adv_u,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5556 ck.SimpleDiffusionJacobian_weak(sdInfo_w_u_rowptr.data(),sdInfo_w_u_colind.data(),mom_wu_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5558 ck.ReactionJacobian_weak(dmom_w_source[0],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5561 elementJacobian_w_v[i][j] +=
5562 ck.AdvectionJacobian_weak(dmom_w_adv_v,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5563 ck.SimpleDiffusionJacobian_weak(sdInfo_w_v_rowptr.data(),sdInfo_w_v_colind.data(),mom_wv_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5565 ck.ReactionJacobian_weak(dmom_w_source[1],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i])
5568 elementJacobian_w_w[i][j] +=
5569 ck.MassJacobian_weak(dmom_w_acc_w_t,vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5570 ck.HamiltonianJacobian_weak(dmom_w_ham_grad_w,&vel_grad_trial[j_nSpace],vel_test_dV[i]) +
5571 ck.AdvectionJacobian_weak(dmom_w_adv_w,vel_trial_ref[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
5572 ck.SimpleDiffusionJacobian_weak(sdInfo_w_w_rowptr.data(),sdInfo_w_w_colind.data(),mom_ww_diff_ten,&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5574 ck.ReactionJacobian_weak(dmom_w_source[2],vel_trial_ref[k*nDOF_trial_element+j],vel_test_dV[i]) +
5577 USE_SUPG*
ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_w[i]) +
5578 ck.NumericalDiffusionJacobian(q_numDiff_w_last[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]) +
5580 ck.NumericalDiffusion(dt*delta*sigma*dV,
5589 for (
int i=0;i<nDOF_test_element;i++)
5591 int eN_i = eN*nDOF_test_element+i;
5592 for (
int j=0;j<nDOF_trial_element;j++)
5594 int eN_i_j = eN_i*nDOF_trial_element+j;
5601 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += element_active*elementJacobian_u_u[i][j];
5602 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_u_v[eN_i_j]] += element_active*elementJacobian_u_v[i][j];
5603 globalJacobian[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_u_w[eN_i_j]] += element_active*elementJacobian_u_w[i][j];
5606 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_v_u[eN_i_j]] += element_active*elementJacobian_v_u[i][j];
5607 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_v_v[eN_i_j]] += element_active*elementJacobian_v_v[i][j];
5608 globalJacobian[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_v_w[eN_i_j]] += element_active*elementJacobian_v_w[i][j];
5611 globalJacobian[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_w_u[eN_i_j]] += element_active*elementJacobian_w_u[i][j];
5612 globalJacobian[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_w_v[eN_i_j]] += element_active*elementJacobian_w_v[i][j];
5613 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_w_w[eN_i_j]] += element_active*elementJacobian_w_w[i][j];
5619 if (ARTIFICIAL_VISCOSITY==3 || ARTIFICIAL_VISCOSITY==4)
5622 for (
int i=0; i<numDOFs_1D; i++)
5625 int u_gi = offset_u+stride_u*i;
5626 int v_gi = offset_v+stride_v*i;
5627 int w_gi = offset_w+stride_w*i;
5630 int u_ith_row_ptr = rowptr[u_gi];
5631 int v_ith_row_ptr = rowptr[v_gi];
5632 int w_ith_row_ptr = rowptr[w_gi];
5635 int numDOFs_ith_row = rowptr_1D[i+1]-rowptr_1D[i];
5636 for (
int counter = 0; counter < numDOFs_ith_row; counter++)
5639 int uu_ij = u_ith_row_ptr + (offset_u + counter*stride_u);
5640 int vv_ij = v_ith_row_ptr + (offset_v + counter*stride_v);
5641 int ww_ij = w_ith_row_ptr + (offset_w + counter*stride_w);
5644 double uStar_dij = uStar_dMatrix[ij];
5645 double vStar_dij = vStar_dMatrix[ij];
5646 double wStar_dij = wStar_dMatrix[ij];
5649 globalJacobian[uu_ij] -= uStar_dij;
5650 globalJacobian[vv_ij] -= vStar_dij;
5651 globalJacobian[ww_ij] -= wStar_dij;
5669 eN_nDOF_trial_element = eN*nDOF_trial_element;
5670 double eps_rho,eps_mu;
5671 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5673 int ebN_kb = ebN*nQuadraturePoints_elementBoundary+kb,
5674 ebN_kb_nSpace = ebN_kb*nSpace,
5675 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5676 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5687 jac_ext[nSpace*nSpace],
5689 jacInv_ext[nSpace*nSpace],
5690 boundaryJac[nSpace*(nSpace-1)],
5691 metricTensor[(nSpace-1)*(nSpace-1)],
5692 metricTensorDetSqrt,
5693 vel_grad_trial_trace[nDOF_trial_element*nSpace],
5695 vel_test_dS[nDOF_test_element],
5697 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
5698 vel_grad_test_dS[nDOF_trial_element*nSpace],
5699 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
5700 ck.calculateMapping_elementBoundary(eN,
5706 mesh_trial_trace_ref.data(),
5707 mesh_grad_trial_trace_ref.data(),
5708 boundaryJac_ref.data(),
5714 metricTensorDetSqrt,
5718 ck.calculateMappingVelocity_elementBoundary(eN,
5722 mesh_velocity_dof.data(),
5724 mesh_trial_trace_ref.data(),
5725 xt_ext,yt_ext,zt_ext,
5730 dS = metricTensorDetSqrt*dS_ref[kb];
5731 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
5734 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
5736 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
5737 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
5738 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],w_ext);
5740 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
5741 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
5742 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
5744 for (
int j=0;j<nDOF_trial_element;j++)
5746 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
5747 for (
int I=0;I<nSpace;I++)
5748 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
5753 ck.calculateGScale(G,normal,h_penalty);
5754 penalty = h_penalty;
5759 double distance[3], P_normal[3], P_tangent[3]={0.0};
5760 if(use_ball_as_particle==1)
5768 P_normal[0],P_normal[1],P_normal[2]);
5770 ball_velocity.data(),ball_angular_velocity.data(),
5772 x_ext-dist*P_normal[0],
5773 y_ext-dist*P_normal[1],
5774 z_ext-dist*P_normal[2],
5775 bc_u_ext,bc_v_ext,bc_w_ext);
5780 dist = ebq_global_phi_solid[ebN_kb];
5781 P_normal[0] = ebq_global_grad_phi_solid[ebN_kb*nSpace+0];
5782 P_normal[1] = ebq_global_grad_phi_solid[ebN_kb*nSpace+1];
5783 P_normal[2] = ebq_global_grad_phi_solid[ebN_kb*nSpace+2];
5784 bc_u_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+0];
5785 bc_v_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+1];
5786 bc_w_ext = ebq_particle_velocity_solid [ebN_kb*nSpace+2];
5788 distance[0] = -P_normal[0]*dist;
5789 distance[1] = -P_normal[1]*dist;
5790 distance[2] = -P_normal[2]*dist;
5797 for(
int i=0;i<3;++i)
5802 assert(h_penalty>0.0);
5803 if (h_penalty < std::abs(dist))
5804 h_penalty = std::abs(dist);
5807 double C_adim =
C_sbm*visco/h_penalty;
5808 double beta_adim =
beta_sbm*h_penalty*visco;
5811 const double zero_vec[3]={0.,0.,0.};
5812 for (
int i=0;i<nDOF_test_element;i++)
5814 int eN_i = eN*nDOF_test_element+i;
5815 double phi_i = vel_test_dS[i];
5816 double* grad_phi_i = &vel_grad_test_dS[i*nSpace+0];
5818 for (
int j=0;j<nDOF_trial_element;j++)
5823 + i*nDOF_trial_element
5826 double phi_j = vel_test_dS[j]/dS;
5827 const double grad_phi_j[3]={vel_grad_test_dS[j*nSpace+0]/dS,
5828 vel_grad_test_dS[j*nSpace+1]/dS,
5829 vel_grad_test_dS[j*nSpace+2]/dS};
5834 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
5836 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
5838 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] +=
5842 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
5843 visco * phi_i * res[0];
5844 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
5845 visco * phi_i * res[1];
5846 globalJacobian[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_eb_u_w[ebN_i_j]] -=
5847 visco * phi_i * res[2];
5850 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
5851 visco * phi_i * res[0] ;
5852 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
5853 visco * phi_i * res[1];
5854 globalJacobian[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_eb_v_w[ebN_i_j]] -=
5855 visco * phi_i * res[2];
5858 globalJacobian[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_eb_w_u[ebN_i_j]] -=
5859 visco * phi_i * res[0] ;
5860 globalJacobian[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_eb_w_v[ebN_i_j]] -=
5861 visco * phi_i * res[1];
5862 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] -=
5863 visco * phi_i * res[2];
5867 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
5868 visco * phi_j * res[0];
5869 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
5870 visco * phi_j * res[1];
5871 globalJacobian[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_eb_u_w[ebN_i_j]] -=
5872 visco * phi_j * res[2];
5875 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
5876 visco * phi_j * res[0] ;
5877 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
5878 visco * phi_j * res[1];
5879 globalJacobian[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_eb_v_w[ebN_i_j]] -=
5880 visco * phi_j * res[2];
5883 globalJacobian[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_eb_w_u[ebN_i_j]] -=
5884 visco * phi_j * res[0] ;
5885 globalJacobian[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_eb_w_v[ebN_i_j]] -=
5886 visco * phi_j * res[1];
5887 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] -=
5888 visco * phi_j * res[2];
5891 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
5892 C_adim*grad_phi_i_dot_d*phi_j;
5893 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
5894 C_adim*grad_phi_i_dot_d*phi_j;
5895 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] +=
5896 C_adim*grad_phi_i_dot_d*phi_j;
5898 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
5899 C_adim*grad_phi_i_dot_d*grad_phi_j_dot_d;
5900 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
5901 C_adim*grad_phi_i_dot_d*grad_phi_j_dot_d;
5902 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] +=
5903 C_adim*grad_phi_i_dot_d*grad_phi_j_dot_d;
5905 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] +=
5906 C_adim*grad_phi_j_dot_d*phi_i;
5907 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] +=
5908 C_adim*grad_phi_j_dot_d*phi_i;
5909 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] +=
5910 C_adim*grad_phi_j_dot_d*phi_i;
5913 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] -=
5914 visco * grad_phi_j_dot_d * res[0];
5915 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] -=
5916 visco * grad_phi_j_dot_d * res[1];
5917 globalJacobian[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_eb_u_w[ebN_i_j]] -=
5918 visco * grad_phi_j_dot_d * res[2];
5921 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] -=
5922 visco * grad_phi_j_dot_d * res[0] ;
5923 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] -=
5924 visco * grad_phi_j_dot_d * res[1];
5925 globalJacobian[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_eb_v_w[ebN_i_j]] -=
5926 visco * grad_phi_j_dot_d * res[2];
5929 globalJacobian[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_eb_w_u[ebN_i_j]] -=
5930 visco * grad_phi_j_dot_d * res[0] ;
5931 globalJacobian[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_eb_w_v[ebN_i_j]] -=
5932 visco * grad_phi_j_dot_d * res[1];
5933 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] -=
5934 visco * grad_phi_j_dot_d * res[2];
5961 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
5963 int ebN = exteriorElementBoundariesArray[ebNE],
5964 eN = elementBoundaryElementsArray[ebN*2+0],
5965 eN_nDOF_trial_element = eN*nDOF_trial_element,
5966 ebN_local = elementBoundaryLocalElementBoundariesArray[ebN*2+0];
5967 double eps_rho,eps_mu;
5968 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5970 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
5971 ebNE_kb_nSpace = ebNE_kb*nSpace,
5972 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5973 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5984 dmom_u_acc_u_ext=0.0,
5986 dmom_v_acc_v_ext=0.0,
5988 dmom_w_acc_w_ext=0.0,
5989 mass_adv_ext[nSpace],
5990 dmass_adv_u_ext[nSpace],
5991 dmass_adv_v_ext[nSpace],
5992 dmass_adv_w_ext[nSpace],
5993 mom_u_adv_ext[nSpace],
5994 dmom_u_adv_u_ext[nSpace],
5995 dmom_u_adv_v_ext[nSpace],
5996 dmom_u_adv_w_ext[nSpace],
5997 mom_v_adv_ext[nSpace],
5998 dmom_v_adv_u_ext[nSpace],
5999 dmom_v_adv_v_ext[nSpace],
6000 dmom_v_adv_w_ext[nSpace],
6001 mom_w_adv_ext[nSpace],
6002 dmom_w_adv_u_ext[nSpace],
6003 dmom_w_adv_v_ext[nSpace],
6004 dmom_w_adv_w_ext[nSpace],
6005 mom_uu_diff_ten_ext[nSpace],
6006 mom_vv_diff_ten_ext[nSpace],
6007 mom_ww_diff_ten_ext[nSpace],
6008 mom_uv_diff_ten_ext[1],
6009 mom_uw_diff_ten_ext[1],
6010 mom_vu_diff_ten_ext[1],
6011 mom_vw_diff_ten_ext[1],
6012 mom_wu_diff_ten_ext[1],
6013 mom_wv_diff_ten_ext[1],
6014 mom_u_source_ext=0.0,
6015 mom_v_source_ext=0.0,
6016 mom_w_source_ext=0.0,
6018 dmom_u_ham_grad_p_ext[nSpace],
6019 dmom_u_ham_grad_u_ext[nSpace],
6021 dmom_v_ham_grad_p_ext[nSpace],
6022 dmom_v_ham_grad_v_ext[nSpace],
6024 dmom_w_ham_grad_p_ext[nSpace],
6025 dmom_w_ham_grad_w_ext[nSpace],
6026 dmom_u_adv_p_ext[nSpace],
6027 dmom_v_adv_p_ext[nSpace],
6028 dmom_w_adv_p_ext[nSpace],
6029 dflux_mass_u_ext=0.0,
6030 dflux_mass_v_ext=0.0,
6031 dflux_mass_w_ext=0.0,
6032 dflux_mom_u_adv_p_ext=0.0,
6033 dflux_mom_u_adv_u_ext=0.0,
6034 dflux_mom_u_adv_v_ext=0.0,
6035 dflux_mom_u_adv_w_ext=0.0,
6036 dflux_mom_v_adv_p_ext=0.0,
6037 dflux_mom_v_adv_u_ext=0.0,
6038 dflux_mom_v_adv_v_ext=0.0,
6039 dflux_mom_v_adv_w_ext=0.0,
6040 dflux_mom_w_adv_p_ext=0.0,
6041 dflux_mom_w_adv_u_ext=0.0,
6042 dflux_mom_w_adv_v_ext=0.0,
6043 dflux_mom_w_adv_w_ext=0.0,
6048 bc_mom_u_acc_ext=0.0,
6049 bc_dmom_u_acc_u_ext=0.0,
6050 bc_mom_v_acc_ext=0.0,
6051 bc_dmom_v_acc_v_ext=0.0,
6052 bc_mom_w_acc_ext=0.0,
6053 bc_dmom_w_acc_w_ext=0.0,
6054 bc_mass_adv_ext[nSpace],
6055 bc_dmass_adv_u_ext[nSpace],
6056 bc_dmass_adv_v_ext[nSpace],
6057 bc_dmass_adv_w_ext[nSpace],
6058 bc_mom_u_adv_ext[nSpace],
6059 bc_dmom_u_adv_u_ext[nSpace],
6060 bc_dmom_u_adv_v_ext[nSpace],
6061 bc_dmom_u_adv_w_ext[nSpace],
6062 bc_mom_v_adv_ext[nSpace],
6063 bc_dmom_v_adv_u_ext[nSpace],
6064 bc_dmom_v_adv_v_ext[nSpace],
6065 bc_dmom_v_adv_w_ext[nSpace],
6066 bc_mom_w_adv_ext[nSpace],
6067 bc_dmom_w_adv_u_ext[nSpace],
6068 bc_dmom_w_adv_v_ext[nSpace],
6069 bc_dmom_w_adv_w_ext[nSpace],
6070 bc_mom_uu_diff_ten_ext[nSpace],
6071 bc_mom_vv_diff_ten_ext[nSpace],
6072 bc_mom_ww_diff_ten_ext[nSpace],
6073 bc_mom_uv_diff_ten_ext[1],
6074 bc_mom_uw_diff_ten_ext[1],
6075 bc_mom_vu_diff_ten_ext[1],
6076 bc_mom_vw_diff_ten_ext[1],
6077 bc_mom_wu_diff_ten_ext[1],
6078 bc_mom_wv_diff_ten_ext[1],
6079 bc_mom_u_source_ext=0.0,
6080 bc_mom_v_source_ext=0.0,
6081 bc_mom_w_source_ext=0.0,
6082 bc_mom_u_ham_ext=0.0,
6083 bc_dmom_u_ham_grad_p_ext[nSpace],
6084 bc_dmom_u_ham_grad_u_ext[nSpace],
6085 bc_mom_v_ham_ext=0.0,
6086 bc_dmom_v_ham_grad_p_ext[nSpace],
6087 bc_dmom_v_ham_grad_v_ext[nSpace],
6088 bc_mom_w_ham_ext=0.0,
6089 bc_dmom_w_ham_grad_p_ext[nSpace],
6090 bc_dmom_w_ham_grad_w_ext[nSpace],
6091 fluxJacobian_p_p[nDOF_trial_element],
6092 fluxJacobian_p_u[nDOF_trial_element],
6093 fluxJacobian_p_v[nDOF_trial_element],
6094 fluxJacobian_p_w[nDOF_trial_element],
6095 fluxJacobian_u_p[nDOF_trial_element],
6096 fluxJacobian_u_u[nDOF_trial_element],
6097 fluxJacobian_u_v[nDOF_trial_element],
6098 fluxJacobian_u_w[nDOF_trial_element],
6099 fluxJacobian_v_p[nDOF_trial_element],
6100 fluxJacobian_v_u[nDOF_trial_element],
6101 fluxJacobian_v_v[nDOF_trial_element],
6102 fluxJacobian_v_w[nDOF_trial_element],
6103 fluxJacobian_w_p[nDOF_trial_element],
6104 fluxJacobian_w_u[nDOF_trial_element],
6105 fluxJacobian_w_v[nDOF_trial_element],
6106 fluxJacobian_w_w[nDOF_trial_element],
6107 jac_ext[nSpace*nSpace],
6109 jacInv_ext[nSpace*nSpace],
6110 boundaryJac[nSpace*(nSpace-1)],
6111 metricTensor[(nSpace-1)*(nSpace-1)],
6112 metricTensorDetSqrt,
6113 p_grad_trial_trace[nDOF_trial_element*nSpace],
6114 vel_grad_trial_trace[nDOF_trial_element*nSpace],
6116 p_test_dS[nDOF_test_element],
6117 vel_test_dS[nDOF_test_element],
6119 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
6120 vel_grad_test_dS[nDOF_trial_element*nSpace],
6124 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
6125 ck.calculateMapping_elementBoundary(eN,
6131 mesh_trial_trace_ref.data(),
6132 mesh_grad_trial_trace_ref.data(),
6133 boundaryJac_ref.data(),
6139 metricTensorDetSqrt,
6143 ck.calculateMappingVelocity_elementBoundary(eN,
6147 mesh_velocity_dof.data(),
6149 mesh_trial_trace_ref.data(),
6150 xt_ext,yt_ext,zt_ext,
6156 dS = metricTensorDetSqrt*dS_ref[kb];
6157 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
6158 ck.calculateGScale(G,&ebqe_normal_phi_ext[ebNE_kb_nSpace],h_phi);
6160 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
6161 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter[eN]);
6162 const double particle_eps = particle_epsFact * (useMetrics * h_phi + (1.0 - useMetrics) * elementDiameter[eN]);
6167 ck.gradTrialFromRef(&vel_grad_trial_trace_ref[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
6170 p_ext = ebqe_p[ebNE_kb];
6171 ck.valFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],u_ext);
6172 ck.valFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],v_ext);
6173 ck.valFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],&vel_trial_trace_ref[ebN_local_kb*nDOF_test_element],w_ext);
6175 for (
int I=0;I<nSpace;I++)
6176 grad_p_ext[I] = ebqe_grad_p[ebNE_kb_nSpace+I];
6177 ck.gradFromDOF(u_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
6178 ck.gradFromDOF(v_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
6179 ck.gradFromDOF(w_dof.data(),&vel_l2g[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
6181 for (
int j=0;j<nDOF_trial_element;j++)
6184 vel_test_dS[j] = vel_test_trace_ref[ebN_local_kb*nDOF_test_element+j]*dS;
6185 for (
int I=0;I<nSpace;I++)
6186 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
6191 bc_p_ext = isDOFBoundary_p[ebNE_kb]*ebqe_bc_p_ext[ebNE_kb]+(1-isDOFBoundary_p[ebNE_kb])*p_ext;
6193 bc_u_ext = isDOFBoundary_u[ebNE_kb]*(ebqe_bc_u_ext[ebNE_kb] + MOVING_DOMAIN*xt_ext) + (1-isDOFBoundary_u[ebNE_kb])*u_ext;
6194 bc_v_ext = isDOFBoundary_v[ebNE_kb]*(ebqe_bc_v_ext[ebNE_kb] + MOVING_DOMAIN*yt_ext) + (1-isDOFBoundary_v[ebNE_kb])*v_ext;
6195 bc_w_ext = isDOFBoundary_w[ebNE_kb]*(ebqe_bc_w_ext[ebNE_kb] + MOVING_DOMAIN*zt_ext) + (1-isDOFBoundary_w[ebNE_kb])*w_ext;
6197 porosity_ext = 1.0 - ebqe_vos_ext[ebNE_kb];
6201 double distance_to_omega_solid = 1e10;
6202 if (use_ball_as_particle == 1)
6204 get_distance_to_ball(nParticles, ball_center.data(), ball_radius.data(), x_ext, y_ext, z_ext, distance_to_omega_solid);
6208 for (
int i = 0; i < nParticles; i++)
6210 double distance_to_i_th_solid = ebq_global_phi_solid[i * nElementBoundaries_global * nQuadraturePoints_elementBoundary + ebNE_kb];
6211 distance_to_omega_solid = (distance_to_i_th_solid < distance_to_omega_solid)?distance_to_i_th_solid:distance_to_omega_solid;
6214 double eddy_viscosity_ext(0.),bc_eddy_viscosity_ext(0.),rhoSave, nuSave;
6223 elementDiameter[eN],
6224 smagorinskyConstant,
6225 turbulenceClosureModel,
6228 ebqe_vf_ext[ebNE_kb],
6229 ebqe_phi_ext[ebNE_kb],
6230 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
6231 distance_to_omega_solid,
6232 ebqe_kappa_phi_ext[ebNE_kb],
6244 ebqe_velocity_star[ebNE_kb_nSpace+0],
6245 ebqe_velocity_star[ebNE_kb_nSpace+1],
6246 ebqe_velocity_star[ebNE_kb_nSpace+2],
6270 mom_uu_diff_ten_ext,
6271 mom_vv_diff_ten_ext,
6272 mom_ww_diff_ten_ext,
6273 mom_uv_diff_ten_ext,
6274 mom_uw_diff_ten_ext,
6275 mom_vu_diff_ten_ext,
6276 mom_vw_diff_ten_ext,
6277 mom_wu_diff_ten_ext,
6278 mom_wv_diff_ten_ext,
6283 dmom_u_ham_grad_p_ext,
6284 dmom_u_ham_grad_u_ext,
6286 dmom_v_ham_grad_p_ext,
6287 dmom_v_ham_grad_v_ext,
6289 dmom_w_ham_grad_p_ext,
6290 dmom_w_ham_grad_w_ext,
6298 MATERIAL_PARAMETERS_AS_FUNCTION,
6299 ebqe_density_as_function[ebNE_kb],
6300 ebqe_dynamic_viscosity_as_function[ebNE_kb],
6303 use_ball_as_particle,
6306 ball_velocity.data(),
6307 ball_angular_velocity.data(),
6308 INT_BY_PARTS_PRESSURE);
6317 elementDiameter[eN],
6318 smagorinskyConstant,
6319 turbulenceClosureModel,
6322 bc_ebqe_vf_ext[ebNE_kb],
6323 bc_ebqe_phi_ext[ebNE_kb],
6324 &ebqe_normal_phi_ext[ebNE_kb_nSpace],
6325 distance_to_omega_solid,
6326 ebqe_kappa_phi_ext[ebNE_kb],
6338 ebqe_velocity_star[ebNE_kb_nSpace+0],
6339 ebqe_velocity_star[ebNE_kb_nSpace+1],
6340 ebqe_velocity_star[ebNE_kb_nSpace+2],
6341 bc_eddy_viscosity_ext,
6343 bc_dmom_u_acc_u_ext,
6345 bc_dmom_v_acc_v_ext,
6347 bc_dmom_w_acc_w_ext,
6353 bc_dmom_u_adv_u_ext,
6354 bc_dmom_u_adv_v_ext,
6355 bc_dmom_u_adv_w_ext,
6357 bc_dmom_v_adv_u_ext,
6358 bc_dmom_v_adv_v_ext,
6359 bc_dmom_v_adv_w_ext,
6361 bc_dmom_w_adv_u_ext,
6362 bc_dmom_w_adv_v_ext,
6363 bc_dmom_w_adv_w_ext,
6364 bc_mom_uu_diff_ten_ext,
6365 bc_mom_vv_diff_ten_ext,
6366 bc_mom_ww_diff_ten_ext,
6367 bc_mom_uv_diff_ten_ext,
6368 bc_mom_uw_diff_ten_ext,
6369 bc_mom_vu_diff_ten_ext,
6370 bc_mom_vw_diff_ten_ext,
6371 bc_mom_wu_diff_ten_ext,
6372 bc_mom_wv_diff_ten_ext,
6373 bc_mom_u_source_ext,
6374 bc_mom_v_source_ext,
6375 bc_mom_w_source_ext,
6377 bc_dmom_u_ham_grad_p_ext,
6378 bc_dmom_u_ham_grad_u_ext,
6380 bc_dmom_v_ham_grad_p_ext,
6381 bc_dmom_v_ham_grad_v_ext,
6383 bc_dmom_w_ham_grad_p_ext,
6384 bc_dmom_w_ham_grad_w_ext,
6392 MATERIAL_PARAMETERS_AS_FUNCTION,
6393 ebqe_density_as_function[ebNE_kb],
6394 ebqe_dynamic_viscosity_as_function[ebNE_kb],
6397 use_ball_as_particle,
6400 ball_velocity.data(),
6401 ball_angular_velocity.data(),
6402 INT_BY_PARTS_PRESSURE);
6404 if (turbulenceClosureModel >= 3)
6406 const double turb_var_grad_0_dummy[3] = {0.,0.,0.};
6407 const double c_mu = 0.09;
6416 ebqe_vf_ext[ebNE_kb],
6417 ebqe_phi_ext[ebNE_kb],
6420 ebqe_turb_var_0[ebNE_kb],
6421 ebqe_turb_var_1[ebNE_kb],
6422 turb_var_grad_0_dummy,
6424 mom_uu_diff_ten_ext,
6425 mom_vv_diff_ten_ext,
6426 mom_ww_diff_ten_ext,
6427 mom_uv_diff_ten_ext,
6428 mom_uw_diff_ten_ext,
6429 mom_vu_diff_ten_ext,
6430 mom_vw_diff_ten_ext,
6431 mom_wu_diff_ten_ext,
6432 mom_wv_diff_ten_ext,
6445 ebqe_vf_ext[ebNE_kb],
6446 ebqe_phi_ext[ebNE_kb],
6449 ebqe_turb_var_0[ebNE_kb],
6450 ebqe_turb_var_1[ebNE_kb],
6451 turb_var_grad_0_dummy,
6452 bc_eddy_viscosity_ext,
6453 bc_mom_uu_diff_ten_ext,
6454 bc_mom_vv_diff_ten_ext,
6455 bc_mom_ww_diff_ten_ext,
6456 bc_mom_uv_diff_ten_ext,
6457 bc_mom_uw_diff_ten_ext,
6458 bc_mom_vu_diff_ten_ext,
6459 bc_mom_vw_diff_ten_ext,
6460 bc_mom_wu_diff_ten_ext,
6461 bc_mom_wv_diff_ten_ext,
6462 bc_mom_u_source_ext,
6463 bc_mom_v_source_ext,
6464 bc_mom_w_source_ext);
6469 mom_u_adv_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*xt_ext;
6470 mom_u_adv_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*yt_ext;
6471 mom_u_adv_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*mom_u_acc_ext*zt_ext;
6472 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
6473 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
6474 dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext;
6476 mom_v_adv_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*xt_ext;
6477 mom_v_adv_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*yt_ext;
6478 mom_v_adv_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*mom_v_acc_ext*zt_ext;
6479 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
6480 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
6481 dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext;
6483 mom_w_adv_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*xt_ext;
6484 mom_w_adv_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*yt_ext;
6485 mom_w_adv_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*mom_w_acc_ext*zt_ext;
6486 dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext;
6487 dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext;
6488 dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext;
6492 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*xt_ext;
6493 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*yt_ext;
6494 bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_u_acc_u_ext*bc_mom_u_acc_ext*zt_ext;
6496 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*xt_ext;
6497 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*yt_ext;
6498 bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_v_acc_v_ext*bc_mom_v_acc_ext*zt_ext;
6500 bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*xt_ext;
6501 bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*yt_ext;
6502 bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*bc_dmom_w_acc_w_ext*bc_mom_w_acc_ext*zt_ext;
6507 isDOFBoundary_u[ebNE_kb],
6508 isDOFBoundary_v[ebNE_kb],
6509 isDOFBoundary_w[ebNE_kb],
6510 isAdvectiveFluxBoundary_p[ebNE_kb],
6511 isAdvectiveFluxBoundary_u[ebNE_kb],
6512 isAdvectiveFluxBoundary_v[ebNE_kb],
6513 isAdvectiveFluxBoundary_w[ebNE_kb],
6514 dmom_u_ham_grad_p_ext[0],
6516 porosity_ext*rhoSave,
6525 ebqe_bc_flux_mass_ext[ebNE_kb]+MOVING_DOMAIN*(xt_ext*normal[0]+yt_ext*normal[1]+zt_ext*normal[2]),
6526 ebqe_bc_flux_mom_u_adv_ext[ebNE_kb],
6527 ebqe_bc_flux_mom_v_adv_ext[ebNE_kb],
6528 ebqe_bc_flux_mom_w_adv_ext[ebNE_kb],
6555 dflux_mom_u_adv_p_ext,
6556 dflux_mom_u_adv_u_ext,
6557 dflux_mom_u_adv_v_ext,
6558 dflux_mom_u_adv_w_ext,
6559 dflux_mom_v_adv_p_ext,
6560 dflux_mom_v_adv_u_ext,
6561 dflux_mom_v_adv_v_ext,
6562 dflux_mom_v_adv_w_ext,
6563 dflux_mom_w_adv_p_ext,
6564 dflux_mom_w_adv_u_ext,
6565 dflux_mom_w_adv_v_ext,
6566 dflux_mom_w_adv_w_ext,
6567 &ebqe_velocity_star[ebNE_kb_nSpace]);
6571 ck.calculateGScale(G,normal,h_penalty);
6572 penalty = useMetrics*C_b/h_penalty + (1.0-useMetrics)*ebqe_penalty_ext[ebNE_kb];
6573 for (
int j=0;j<nDOF_trial_element;j++)
6575 int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
6582 fluxJacobian_u_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6584 ebqe_phi_ext[ebNE_kb],
6585 sdInfo_u_u_rowptr.data(),
6586 sdInfo_u_u_colind.data(),
6587 isDOFBoundary_u[ebNE_kb],
6588 isDiffusiveFluxBoundary_u[ebNE_kb],
6590 mom_uu_diff_ten_ext,
6591 vel_trial_trace_ref[ebN_local_kb_j],
6592 &vel_grad_trial_trace[j_nSpace],
6594 fluxJacobian_u_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6596 ebqe_phi_ext[ebNE_kb],
6597 sdInfo_u_v_rowptr.data(),
6598 sdInfo_u_v_colind.data(),
6599 isDOFBoundary_v[ebNE_kb],
6600 isDiffusiveFluxBoundary_v[ebNE_kb],
6602 mom_uv_diff_ten_ext,
6603 vel_trial_trace_ref[ebN_local_kb_j],
6604 &vel_grad_trial_trace[j_nSpace],
6606 fluxJacobian_u_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_w_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6608 ebqe_phi_ext[ebNE_kb],
6609 sdInfo_u_w_rowptr.data(),
6610 sdInfo_u_w_colind.data(),
6611 isDOFBoundary_w[ebNE_kb],
6612 isDiffusiveFluxBoundary_u[ebNE_kb],
6614 mom_uw_diff_ten_ext,
6615 vel_trial_trace_ref[ebN_local_kb_j],
6616 &vel_grad_trial_trace[j_nSpace],
6620 fluxJacobian_v_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6622 ebqe_phi_ext[ebNE_kb],
6623 sdInfo_v_u_rowptr.data(),
6624 sdInfo_v_u_colind.data(),
6625 isDOFBoundary_u[ebNE_kb],
6626 isDiffusiveFluxBoundary_u[ebNE_kb],
6628 mom_vu_diff_ten_ext,
6629 vel_trial_trace_ref[ebN_local_kb_j],
6630 &vel_grad_trial_trace[j_nSpace],
6632 fluxJacobian_v_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6634 ebqe_phi_ext[ebNE_kb],
6635 sdInfo_v_v_rowptr.data(),
6636 sdInfo_v_v_colind.data(),
6637 isDOFBoundary_v[ebNE_kb],
6638 isDiffusiveFluxBoundary_v[ebNE_kb],
6640 mom_vv_diff_ten_ext,
6641 vel_trial_trace_ref[ebN_local_kb_j],
6642 &vel_grad_trial_trace[j_nSpace],
6644 fluxJacobian_v_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_w_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6646 ebqe_phi_ext[ebNE_kb],
6647 sdInfo_v_w_rowptr.data(),
6648 sdInfo_v_w_colind.data(),
6649 isDOFBoundary_w[ebNE_kb],
6650 isDiffusiveFluxBoundary_v[ebNE_kb],
6652 mom_vw_diff_ten_ext,
6653 vel_trial_trace_ref[ebN_local_kb_j],
6654 &vel_grad_trial_trace[j_nSpace],
6658 fluxJacobian_w_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_u_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6660 ebqe_phi_ext[ebNE_kb],
6661 sdInfo_w_u_rowptr.data(),
6662 sdInfo_w_u_colind.data(),
6663 isDOFBoundary_u[ebNE_kb],
6664 isDiffusiveFluxBoundary_w[ebNE_kb],
6666 mom_wu_diff_ten_ext,
6667 vel_trial_trace_ref[ebN_local_kb_j],
6668 &vel_grad_trial_trace[j_nSpace],
6670 fluxJacobian_w_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_v_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6672 ebqe_phi_ext[ebNE_kb],
6673 sdInfo_w_v_rowptr.data(),
6674 sdInfo_w_v_colind.data(),
6675 isDOFBoundary_v[ebNE_kb],
6676 isDiffusiveFluxBoundary_w[ebNE_kb],
6678 mom_wv_diff_ten_ext,
6679 vel_trial_trace_ref[ebN_local_kb_j],
6680 &vel_grad_trial_trace[j_nSpace],
6682 fluxJacobian_w_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_w_ext,vel_trial_trace_ref[ebN_local_kb_j]) +
6684 ebqe_phi_ext[ebNE_kb],
6685 sdInfo_w_w_rowptr.data(),
6686 sdInfo_w_w_colind.data(),
6687 isDOFBoundary_w[ebNE_kb],
6688 isDiffusiveFluxBoundary_w[ebNE_kb],
6690 mom_ww_diff_ten_ext,
6691 vel_trial_trace_ref[ebN_local_kb_j],
6692 &vel_grad_trial_trace[j_nSpace],
6698 for (
int i=0;i<nDOF_test_element;i++)
6700 int eN_i = eN*nDOF_test_element+i;
6701 for (
int j=0;j<nDOF_trial_element;j++)
6703 int ebN_i_j = ebN*4*
nDOF_test_X_trial_element + i*nDOF_trial_element + j,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
6711 globalJacobian[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*vel_test_dS[i]+
6712 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u[ebNE_kb],
6713 isDiffusiveFluxBoundary_u[ebNE_kb],
6715 vel_trial_trace_ref[ebN_local_kb_j],
6717 sdInfo_u_u_rowptr.data(),
6718 sdInfo_u_u_colind.data(),
6719 mom_uu_diff_ten_ext,
6720 &vel_grad_test_dS[i*nSpace]);
6721 globalJacobian[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] += fluxJacobian_u_v[j]*vel_test_dS[i]+
6722 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v[ebNE_kb],
6723 isDiffusiveFluxBoundary_u[ebNE_kb],
6725 vel_trial_trace_ref[ebN_local_kb_j],
6727 sdInfo_u_v_rowptr.data(),
6728 sdInfo_u_v_colind.data(),
6729 mom_uv_diff_ten_ext,
6730 &vel_grad_test_dS[i*nSpace]);
6731 globalJacobian[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_eb_u_w[ebN_i_j]] += fluxJacobian_u_w[j]*vel_test_dS[i]+
6732 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w[ebNE_kb],
6733 isDiffusiveFluxBoundary_u[ebNE_kb],
6735 vel_trial_trace_ref[ebN_local_kb_j],
6737 sdInfo_u_w_rowptr.data(),
6738 sdInfo_u_w_colind.data(),
6739 mom_uw_diff_ten_ext,
6740 &vel_grad_test_dS[i*nSpace]);
6743 globalJacobian[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] += fluxJacobian_v_u[j]*vel_test_dS[i]+
6744 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u[ebNE_kb],
6745 isDiffusiveFluxBoundary_v[ebNE_kb],
6747 vel_trial_trace_ref[ebN_local_kb_j],
6749 sdInfo_v_u_rowptr.data(),
6750 sdInfo_v_u_colind.data(),
6751 mom_vu_diff_ten_ext,
6752 &vel_grad_test_dS[i*nSpace]);
6753 globalJacobian[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] += fluxJacobian_v_v[j]*vel_test_dS[i]+
6754 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v[ebNE_kb],
6755 isDiffusiveFluxBoundary_v[ebNE_kb],
6757 vel_trial_trace_ref[ebN_local_kb_j],
6759 sdInfo_v_v_rowptr.data(),
6760 sdInfo_v_v_colind.data(),
6761 mom_vv_diff_ten_ext,
6762 &vel_grad_test_dS[i*nSpace]);
6763 globalJacobian[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_eb_v_w[ebN_i_j]] += fluxJacobian_v_w[j]*vel_test_dS[i]+
6764 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w[ebNE_kb],
6765 isDiffusiveFluxBoundary_v[ebNE_kb],
6767 vel_trial_trace_ref[ebN_local_kb_j],
6769 sdInfo_v_w_rowptr.data(),
6770 sdInfo_v_w_colind.data(),
6771 mom_vw_diff_ten_ext,
6772 &vel_grad_test_dS[i*nSpace]);
6775 globalJacobian[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_eb_w_u[ebN_i_j]] += fluxJacobian_w_u[j]*vel_test_dS[i]+
6776 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u[ebNE_kb],
6777 isDiffusiveFluxBoundary_w[ebNE_kb],
6779 vel_trial_trace_ref[ebN_local_kb_j],
6781 sdInfo_w_u_rowptr.data(),
6782 sdInfo_w_u_colind.data(),
6783 mom_wu_diff_ten_ext,
6784 &vel_grad_test_dS[i*nSpace]);
6785 globalJacobian[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_eb_w_v[ebN_i_j]] += fluxJacobian_w_v[j]*vel_test_dS[i]+
6786 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v[ebNE_kb],
6787 isDiffusiveFluxBoundary_w[ebNE_kb],
6789 vel_trial_trace_ref[ebN_local_kb_j],
6791 sdInfo_w_v_rowptr.data(),
6792 sdInfo_w_v_colind.data(),
6793 mom_wv_diff_ten_ext,
6794 &vel_grad_test_dS[i*nSpace]);
6795 globalJacobian[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] += fluxJacobian_w_w[j]*vel_test_dS[i]+
6796 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w[ebNE_kb],
6797 isDiffusiveFluxBoundary_w[ebNE_kb],
6799 vel_trial_trace_ref[ebN_local_kb_j],
6801 sdInfo_w_w_rowptr.data(),
6802 sdInfo_w_w_colind.data(),
6803 mom_ww_diff_ten_ext,
6804 &vel_grad_test_dS[i*nSpace]);