1380 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
1381 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
1382 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
1383 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
1384 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
1385 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
1386 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
1387 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
1388 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
1389 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
1390 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
1391 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
1392 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
1393 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
1394 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
1395 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
1396 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
1397 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
1398 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
1399 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
1400 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
1401 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
1402 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
1403 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
1404 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
1405 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
1406 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
1407 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
1408 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
1409 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
1410 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
1411 double hFactor = args.
scalar<
double>(
"hFactor");
1412 int nElements_global = args.
scalar<
int>(
"nElements_global");
1413 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
1414 double useRBLES = args.
scalar<
double>(
"useRBLES");
1415 double useMetrics = args.
scalar<
double>(
"useMetrics");
1416 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
1417 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
1418 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
1419 double sigma = args.
scalar<
double>(
"sigma");
1424 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
1425 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
1426 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
1427 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
1428 double C_dc = args.
scalar<
double>(
"C_dc");
1429 double C_b = args.
scalar<
double>(
"C_b");
1430 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
1431 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
1432 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
1433 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
1434 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
1435 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
1436 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
1437 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
1438 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
1439 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
1440 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
1441 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
1442 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
1443 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
1444 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
1445 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
1446 xt::pyarray<double>& g = args.
array<
double>(
"g");
1447 const double useVF = args.
scalar<
double>(
"useVF");
1448 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
1449 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
1450 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
1451 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
1452 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
1453 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
1454 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
1455 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
1456 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
1457 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
1458 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
1459 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
1460 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
1461 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
1462 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
1463 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
1464 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
1465 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
1466 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
1467 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
1468 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
1469 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
1470 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
1471 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
1472 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
1473 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
1474 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
1475 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
1476 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
1477 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
1478 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
1479 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
1480 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
1481 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
1482 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
1483 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
1484 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
1485 int offset_p = args.
scalar<
int>(
"offset_p");
1486 int offset_u = args.
scalar<
int>(
"offset_u");
1487 int offset_v = args.
scalar<
int>(
"offset_v");
1488 int offset_w = args.
scalar<
int>(
"offset_w");
1489 int stride_p = args.
scalar<
int>(
"stride_p");
1490 int stride_u = args.
scalar<
int>(
"stride_u");
1491 int stride_v = args.
scalar<
int>(
"stride_v");
1492 int stride_w = args.
scalar<
int>(
"stride_w");
1493 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
1494 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
1495 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
1496 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
1497 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
1498 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
1499 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
1500 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
1501 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
1502 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
1503 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
1504 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
1505 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
1506 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
1507 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
1508 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
1509 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
1510 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
1511 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
1512 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
1513 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
1514 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
1515 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
1516 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
1517 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
1518 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
1519 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
1520 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
1521 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
1522 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
1523 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
1524 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
1525 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
1526 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
1527 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
1528 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
1529 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
1530 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
1531 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
1532 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
1533 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
1534 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
1535 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
1536 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
1537 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
1538 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
1539 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
1540 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
1541 xt::pyarray<double>& q_dragBeam1 = args.
array<
double>(
"q_dragBeam1");
1542 xt::pyarray<double>& q_dragBeam2 = args.
array<
double>(
"q_dragBeam2");
1543 xt::pyarray<double>& q_dragBeam3 = args.
array<
double>(
"q_dragBeam3");
1544 xt::pyarray<double>& ebqe_dragBeam1 = args.
array<
double>(
"ebqe_dragBeam1");
1545 xt::pyarray<double>& ebqe_dragBeam2 = args.
array<
double>(
"ebqe_dragBeam2");
1546 xt::pyarray<double>& ebqe_dragBeam3 = args.
array<
double>(
"ebqe_dragBeam3");
1550 double globalConservationError=0.0;
1551 for(
int eN=0;eN<nElements_global;eN++)
1554 double elementResidual_p[nDOF_test_element],
1555 elementResidual_u[nDOF_test_element],
1556 elementResidual_v[nDOF_test_element],
1557 elementResidual_w[nDOF_test_element],
1559 for (
int i=0;i<nDOF_test_element;i++)
1561 int eN_i = eN*nDOF_test_element+i;
1562 elementResidual_p_save.data()[eN_i]=0.0;
1563 elementResidual_p[i]=0.0;
1564 elementResidual_u[i]=0.0;
1565 elementResidual_v[i]=0.0;
1566 elementResidual_w[i]=0.0;
1571 for(
int k=0;k<nQuadraturePoints_element;k++)
1574 int eN_k = eN*nQuadraturePoints_element+k,
1575 eN_k_nSpace = eN_k*nSpace,
1576 eN_nDOF_trial_element = eN*nDOF_trial_element;
1577 double p=0.0,
u=0.0,
v=0.0,
w=0.0,
1578 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
1586 dmass_adv_u[nSpace],
1587 dmass_adv_v[nSpace],
1588 dmass_adv_w[nSpace],
1590 dmom_u_adv_u[nSpace],
1591 dmom_u_adv_v[nSpace],
1592 dmom_u_adv_w[nSpace],
1594 dmom_v_adv_u[nSpace],
1595 dmom_v_adv_v[nSpace],
1596 dmom_v_adv_w[nSpace],
1598 dmom_w_adv_u[nSpace],
1599 dmom_w_adv_v[nSpace],
1600 dmom_w_adv_w[nSpace],
1601 mom_uu_diff_ten[nSpace],
1602 mom_vv_diff_ten[nSpace],
1603 mom_ww_diff_ten[nSpace],
1614 dmom_u_ham_grad_p[nSpace],
1616 dmom_v_ham_grad_p[nSpace],
1618 dmom_w_ham_grad_p[nSpace],
1629 Lstar_u_p[nDOF_test_element],
1630 Lstar_v_p[nDOF_test_element],
1631 Lstar_w_p[nDOF_test_element],
1632 Lstar_u_u[nDOF_test_element],
1633 Lstar_v_v[nDOF_test_element],
1634 Lstar_w_w[nDOF_test_element],
1635 Lstar_p_u[nDOF_test_element],
1636 Lstar_p_v[nDOF_test_element],
1637 Lstar_p_w[nDOF_test_element],
1642 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
1643 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
1646 jacInv[nSpace*nSpace],
1647 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
1648 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
1649 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
1655 dmom_u_source[nSpace],
1656 dmom_v_source[nSpace],
1657 dmom_w_source[nSpace],
1659 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
1661 ck.calculateMapping_element(eN,
1665 mesh_trial_ref.data(),
1666 mesh_grad_trial_ref.data(),
1671 ck.calculateH_element(eN,
1673 nodeDiametersArray.data(),
1675 mesh_trial_ref.data(),
1678 ck.calculateMappingVelocity_element(eN,
1680 mesh_velocity_dof.data(),
1682 mesh_trial_ref.data(),
1687 dV = fabs(jacDet)*dV_ref.data()[k];
1688 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1691 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1692 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
1695 ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial);
1696 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
1698 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_ref.data()[k*nDOF_trial_element],p);
1699 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
u);
1700 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
v);
1701 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
w);
1703 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial,grad_p);
1704 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
1705 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
1706 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_w);
1708 for (
int j=0;j<nDOF_trial_element;j++)
1710 p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV;
1711 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
1712 for (
int I=0;I<nSpace;I++)
1714 p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;
1715 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
1719 porosity = q_porosity.data()[eN_k];
1723 q_velocity.data()[eN_k_nSpace+0]=
u;
1724 q_velocity.data()[eN_k_nSpace+1]=
v;
1725 q_velocity.data()[eN_k_nSpace+2]=
w;
1726 q_x.data()[eN_k_nSpace+0]=x;
1727 q_x.data()[eN_k_nSpace+1]=y;
1728 q_x.data()[eN_k_nSpace+2]=
z;
1739 elementDiameter.data()[eN],
1740 smagorinskyConstant,
1741 turbulenceClosureModel,
1746 &normal_phi.data()[eN_k_nSpace],
1747 kappa_phi.data()[eN_k],
1799 q_dragBeam1.data()[eN_k],
1800 q_dragBeam2.data()[eN_k],
1801 q_dragBeam3.data()[eN_k]);
1803 mass_source = q_mass_source.data()[eN_k];
1809 q_dragAlpha.data()[eN_k],
1810 q_dragBeta.data()[eN_k],
1823 eps_solid.data()[0],
1824 phi_solid.data()[eN_k],
1825 q_velocity_solid.data()[eN_k_nSpace+0],
1826 q_velocity_solid.data()[eN_k_nSpace+1],
1827 q_velocity_solid.data()[eN_k_nSpace+2],
1836 if (turbulenceClosureModel >= 3)
1838 const double c_mu = 0.09;
1851 q_turb_var_0.data()[eN_k],
1852 q_turb_var_1.data()[eN_k],
1853 &q_turb_var_grad_0.data()[eN_k_nSpace],
1871 q_mom_u_acc.data()[eN_k] = mom_u_acc;
1872 q_mom_v_acc.data()[eN_k] = mom_v_acc;
1873 q_mom_w_acc.data()[eN_k] = mom_w_acc;
1875 q_mass_adv.data()[eN_k_nSpace+0] =
u;
1876 q_mass_adv.data()[eN_k_nSpace+1] =
v;
1877 q_mass_adv.data()[eN_k_nSpace+2] =
w;
1882 mass_adv[0] -= MOVING_DOMAIN*
xt;
1883 mass_adv[1] -= MOVING_DOMAIN*yt;
1884 mass_adv[2] -= MOVING_DOMAIN*zt;
1886 mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*
xt;
1887 mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt;
1888 mom_u_adv[2] -= MOVING_DOMAIN*mom_u_acc*zt;
1889 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
1890 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
1891 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
1893 mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*
xt;
1894 mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt;
1895 mom_v_adv[2] -= MOVING_DOMAIN*mom_v_acc*zt;
1896 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
1897 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
1898 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
1900 mom_w_adv[0] -= MOVING_DOMAIN*mom_w_acc*
xt;
1901 mom_w_adv[1] -= MOVING_DOMAIN*mom_w_acc*yt;
1902 mom_w_adv[2] -= MOVING_DOMAIN*mom_w_acc*zt;
1903 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
1904 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
1905 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
1910 q_mom_u_acc_beta_bdf.data()[eN_k],
1916 q_mom_v_acc_beta_bdf.data()[eN_k],
1922 q_mom_w_acc_beta_bdf.data()[eN_k],
1931 pdeResidual_p =
ck.Advection_strong(dmass_adv_u,grad_u) +
1933 ck.Reaction_strong(mass_source) +
1935 ck.Advection_strong(dmass_adv_v,grad_v) +
1936 ck.Advection_strong(dmass_adv_w,grad_w);
1938 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
1939 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
1940 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
1942 pdeResidual_u =
ck.Mass_strong(mom_u_acc_t) +
1943 ck.Advection_strong(dmom_adv_sge,grad_u) +
1944 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
1945 ck.Reaction_strong(mom_u_source);
1947 pdeResidual_v =
ck.Mass_strong(mom_v_acc_t) +
1948 ck.Advection_strong(dmom_adv_sge,grad_v) +
1949 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
1950 ck.Reaction_strong(mom_v_source);
1952 pdeResidual_w =
ck.Mass_strong(mom_w_acc_t) +
1953 ck.Advection_strong(dmom_adv_sge,grad_w) +
1954 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
1955 ck.Reaction_strong(mom_w_source);
1959 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
1961 elementDiameter.data()[eN],
1966 dmom_u_ham_grad_p[0],
1969 q_cfl.data()[eN_k]);
1976 dmom_u_ham_grad_p[0],
1979 q_cfl.data()[eN_k]);
1981 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
1982 tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0;
1995 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
1996 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
1997 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
1999 mom_u_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+0]);
2000 mom_u_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+0]);
2001 mom_u_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge.data()[eN_k_nSpace+0]);
2003 mom_v_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+1]);
2004 mom_v_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+1]);
2005 mom_v_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge.data()[eN_k_nSpace+1]);
2007 mom_w_adv[0] += dmom_u_acc_u*(useRBLES*subgridError_u*q_velocity_sge.data()[eN_k_nSpace+2]);
2008 mom_w_adv[1] += dmom_u_acc_u*(useRBLES*subgridError_v*q_velocity_sge.data()[eN_k_nSpace+2]);
2009 mom_w_adv[2] += dmom_u_acc_u*(useRBLES*subgridError_w*q_velocity_sge.data()[eN_k_nSpace+2]);
2012 for (
int i=0;i<nDOF_test_element;i++)
2014 int i_nSpace = i*nSpace;
2015 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
2016 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
2017 Lstar_w_p[i]=
ck.Advection_adjoint(dmass_adv_w,&p_grad_test_dV[i_nSpace]);
2019 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2020 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2021 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
2022 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2023 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2024 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
2027 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
2028 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
2029 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
2033 norm_Rv = sqrt(pdeResidual_u*pdeResidual_u + pdeResidual_v*pdeResidual_v + pdeResidual_w*pdeResidual_w);
2034 q_numDiff_u.data()[eN_k] = C_dc*norm_Rv*(useMetrics/sqrt(G_dd_G+1.0e-12) +
2035 (1.0-useMetrics)*hFactor*hFactor*elementDiameter.data()[eN]*elementDiameter.data()[eN]);
2036 q_numDiff_v.data()[eN_k] = q_numDiff_u.data()[eN_k];
2037 q_numDiff_w.data()[eN_k] = q_numDiff_u.data()[eN_k];
2041 for(
int i=0;i<nDOF_test_element;i++)
2043 int i_nSpace=i*nSpace;
2045 elementResidual_p[i] +=
ck.Advection_weak(mass_adv,&p_grad_test_dV[i_nSpace]) +
2047 ck.Reaction_weak(mass_source,p_test_dV[i]) +
2049 ck.SubgridError(subgridError_u,Lstar_u_p[i]) +
2050 ck.SubgridError(subgridError_v,Lstar_v_p[i]) +
2051 ck.SubgridError(subgridError_w,Lstar_w_p[i]);
2053 elementResidual_u[i] +=
ck.Mass_weak(mom_u_acc_t,vel_test_dV[i]) +
2054 ck.Advection_weak(mom_u_adv,&vel_grad_test_dV[i_nSpace]) +
2055 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]) +
2056 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]) +
2057 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]) +
2058 ck.Reaction_weak(mom_u_source,vel_test_dV[i]) +
2059 ck.Hamiltonian_weak(mom_u_ham,vel_test_dV[i]) +
2060 ck.SubgridError(subgridError_p,Lstar_p_u[i]) +
2061 ck.SubgridError(subgridError_u,Lstar_u_u[i]) +
2062 ck.NumericalDiffusion(q_numDiff_u_last.data()[eN_k],grad_u,&vel_grad_test_dV[i_nSpace]);
2064 elementResidual_v[i] +=
ck.Mass_weak(mom_v_acc_t,vel_test_dV[i]) +
2065 ck.Advection_weak(mom_v_adv,&vel_grad_test_dV[i_nSpace]) +
2066 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]) +
2067 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]) +
2068 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]) +
2069 ck.Reaction_weak(mom_v_source,vel_test_dV[i]) +
2070 ck.Hamiltonian_weak(mom_v_ham,vel_test_dV[i]) +
2071 ck.SubgridError(subgridError_p,Lstar_p_v[i]) +
2072 ck.SubgridError(subgridError_v,Lstar_v_v[i]) +
2073 ck.NumericalDiffusion(q_numDiff_v_last.data()[eN_k],grad_v,&vel_grad_test_dV[i_nSpace]);
2075 elementResidual_w[i] +=
ck.Mass_weak(mom_w_acc_t,vel_test_dV[i]) +
2076 ck.Advection_weak(mom_w_adv,&vel_grad_test_dV[i_nSpace]) +
2077 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]) +
2078 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]) +
2079 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]) +
2080 ck.Reaction_weak(mom_w_source,vel_test_dV[i]) +
2081 ck.Hamiltonian_weak(mom_w_ham,vel_test_dV[i]) +
2082 ck.SubgridError(subgridError_p,Lstar_p_w[i]) +
2083 ck.SubgridError(subgridError_w,Lstar_w_w[i]) +
2084 ck.NumericalDiffusion(q_numDiff_w_last.data()[eN_k],grad_w,&vel_grad_test_dV[i_nSpace]);
2090 for(
int i=0;i<nDOF_test_element;i++)
2092 int eN_i=eN*nDOF_test_element+i;
2094 elementResidual_p_save.data()[eN_i] += elementResidual_p[i];
2096 globalResidual.data()[offset_p+stride_p*p_l2g.data()[eN_i]]+=elementResidual_p[i];
2097 globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
2098 globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
2099 globalResidual.data()[offset_w+stride_w*vel_l2g.data()[eN_i]]+=elementResidual_w[i];
2108 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
2110 int ebN = exteriorElementBoundariesArray.data()[ebNE],
2111 eN = elementBoundaryElementsArray.data()[ebN*2+0],
2112 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
2113 eN_nDOF_trial_element = eN*nDOF_trial_element;
2114 double elementResidual_p[nDOF_test_element],
2115 elementResidual_u[nDOF_test_element],
2116 elementResidual_v[nDOF_test_element],
2117 elementResidual_w[nDOF_test_element],
2119 for (
int i=0;i<nDOF_test_element;i++)
2121 elementResidual_p[i]=0.0;
2122 elementResidual_u[i]=0.0;
2123 elementResidual_v[i]=0.0;
2124 elementResidual_w[i]=0.0;
2126 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
2128 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
2129 ebNE_kb_nSpace = ebNE_kb*nSpace,
2130 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
2131 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
2141 dmom_u_acc_u_ext=0.0,
2143 dmom_v_acc_v_ext=0.0,
2145 dmom_w_acc_w_ext=0.0,
2146 mass_adv_ext[nSpace],
2147 dmass_adv_u_ext[nSpace],
2148 dmass_adv_v_ext[nSpace],
2149 dmass_adv_w_ext[nSpace],
2150 mom_u_adv_ext[nSpace],
2151 dmom_u_adv_u_ext[nSpace],
2152 dmom_u_adv_v_ext[nSpace],
2153 dmom_u_adv_w_ext[nSpace],
2154 mom_v_adv_ext[nSpace],
2155 dmom_v_adv_u_ext[nSpace],
2156 dmom_v_adv_v_ext[nSpace],
2157 dmom_v_adv_w_ext[nSpace],
2158 mom_w_adv_ext[nSpace],
2159 dmom_w_adv_u_ext[nSpace],
2160 dmom_w_adv_v_ext[nSpace],
2161 dmom_w_adv_w_ext[nSpace],
2162 mom_uu_diff_ten_ext[nSpace],
2163 mom_vv_diff_ten_ext[nSpace],
2164 mom_ww_diff_ten_ext[nSpace],
2165 mom_uv_diff_ten_ext[1],
2166 mom_uw_diff_ten_ext[1],
2167 mom_vu_diff_ten_ext[1],
2168 mom_vw_diff_ten_ext[1],
2169 mom_wu_diff_ten_ext[1],
2170 mom_wv_diff_ten_ext[1],
2171 mom_u_source_ext=0.0,
2172 mom_v_source_ext=0.0,
2173 mom_w_source_ext=0.0,
2175 dmom_u_ham_grad_p_ext[nSpace],
2177 dmom_v_ham_grad_p_ext[nSpace],
2179 dmom_w_ham_grad_p_ext[nSpace],
2180 dmom_u_adv_p_ext[nSpace],
2181 dmom_v_adv_p_ext[nSpace],
2182 dmom_w_adv_p_ext[nSpace],
2184 flux_mom_u_adv_ext=0.0,
2185 flux_mom_v_adv_ext=0.0,
2186 flux_mom_w_adv_ext=0.0,
2187 flux_mom_uu_diff_ext=0.0,
2188 flux_mom_uv_diff_ext=0.0,
2189 flux_mom_uw_diff_ext=0.0,
2190 flux_mom_vu_diff_ext=0.0,
2191 flux_mom_vv_diff_ext=0.0,
2192 flux_mom_vw_diff_ext=0.0,
2193 flux_mom_wu_diff_ext=0.0,
2194 flux_mom_wv_diff_ext=0.0,
2195 flux_mom_ww_diff_ext=0.0,
2200 bc_mom_u_acc_ext=0.0,
2201 bc_dmom_u_acc_u_ext=0.0,
2202 bc_mom_v_acc_ext=0.0,
2203 bc_dmom_v_acc_v_ext=0.0,
2204 bc_mom_w_acc_ext=0.0,
2205 bc_dmom_w_acc_w_ext=0.0,
2206 bc_mass_adv_ext[nSpace],
2207 bc_dmass_adv_u_ext[nSpace],
2208 bc_dmass_adv_v_ext[nSpace],
2209 bc_dmass_adv_w_ext[nSpace],
2210 bc_mom_u_adv_ext[nSpace],
2211 bc_dmom_u_adv_u_ext[nSpace],
2212 bc_dmom_u_adv_v_ext[nSpace],
2213 bc_dmom_u_adv_w_ext[nSpace],
2214 bc_mom_v_adv_ext[nSpace],
2215 bc_dmom_v_adv_u_ext[nSpace],
2216 bc_dmom_v_adv_v_ext[nSpace],
2217 bc_dmom_v_adv_w_ext[nSpace],
2218 bc_mom_w_adv_ext[nSpace],
2219 bc_dmom_w_adv_u_ext[nSpace],
2220 bc_dmom_w_adv_v_ext[nSpace],
2221 bc_dmom_w_adv_w_ext[nSpace],
2222 bc_mom_uu_diff_ten_ext[nSpace],
2223 bc_mom_vv_diff_ten_ext[nSpace],
2224 bc_mom_ww_diff_ten_ext[nSpace],
2225 bc_mom_uv_diff_ten_ext[1],
2226 bc_mom_uw_diff_ten_ext[1],
2227 bc_mom_vu_diff_ten_ext[1],
2228 bc_mom_vw_diff_ten_ext[1],
2229 bc_mom_wu_diff_ten_ext[1],
2230 bc_mom_wv_diff_ten_ext[1],
2231 bc_mom_u_source_ext=0.0,
2232 bc_mom_v_source_ext=0.0,
2233 bc_mom_w_source_ext=0.0,
2234 bc_mom_u_ham_ext=0.0,
2235 bc_dmom_u_ham_grad_p_ext[nSpace],
2236 bc_mom_v_ham_ext=0.0,
2237 bc_dmom_v_ham_grad_p_ext[nSpace],
2238 bc_mom_w_ham_ext=0.0,
2239 bc_dmom_w_ham_grad_p_ext[nSpace],
2240 jac_ext[nSpace*nSpace],
2242 jacInv_ext[nSpace*nSpace],
2243 boundaryJac[nSpace*(nSpace-1)],
2244 metricTensor[(nSpace-1)*(nSpace-1)],
2245 metricTensorDetSqrt,
2246 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
2247 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
2248 vel_grad_test_dS[nDOF_trial_element*nSpace],
2249 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
2253 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
2254 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;
2256 ck.calculateMapping_elementBoundary(eN,
2262 mesh_trial_trace_ref.data(),
2263 mesh_grad_trial_trace_ref.data(),
2264 boundaryJac_ref.data(),
2270 metricTensorDetSqrt,
2274 ck.calculateMappingVelocity_elementBoundary(eN,
2278 mesh_velocity_dof.data(),
2280 mesh_trial_trace_ref.data(),
2281 xt_ext,yt_ext,zt_ext,
2291 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
2294 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
2295 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
2297 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2298 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
2302 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
2303 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
2306 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_ext);
2307 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
2308 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
2309 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext);
2310 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
2311 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
2312 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
2313 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
2315 for (
int j=0;j<nDOF_trial_element;j++)
2317 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2318 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
2319 for (
int I=0;I<nSpace;I++)
2320 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
2322 bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
2323 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
2324 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
2325 bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext;
2327 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
2338 elementDiameter.data()[eN],
2339 smagorinskyConstant,
2340 turbulenceClosureModel,
2343 ebqe_vf_ext.data()[ebNE_kb],
2344 ebqe_phi_ext.data()[ebNE_kb],
2345 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
2346 ebqe_kappa_phi_ext.data()[ebNE_kb],
2380 mom_uu_diff_ten_ext,
2381 mom_vv_diff_ten_ext,
2382 mom_ww_diff_ten_ext,
2383 mom_uv_diff_ten_ext,
2384 mom_uw_diff_ten_ext,
2385 mom_vu_diff_ten_ext,
2386 mom_vw_diff_ten_ext,
2387 mom_wu_diff_ten_ext,
2388 mom_wv_diff_ten_ext,
2393 dmom_u_ham_grad_p_ext,
2395 dmom_v_ham_grad_p_ext,
2397 dmom_w_ham_grad_p_ext,
2398 ebqe_dragBeam1.data()[ebNE_kb],
2399 ebqe_dragBeam2.data()[ebNE_kb],
2400 ebqe_dragBeam3.data()[ebNE_kb]);
2408 elementDiameter.data()[eN],
2409 smagorinskyConstant,
2410 turbulenceClosureModel,
2413 bc_ebqe_vf_ext.data()[ebNE_kb],
2414 bc_ebqe_phi_ext.data()[ebNE_kb],
2415 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
2416 ebqe_kappa_phi_ext.data()[ebNE_kb],
2429 bc_dmom_u_acc_u_ext,
2431 bc_dmom_v_acc_v_ext,
2433 bc_dmom_w_acc_w_ext,
2439 bc_dmom_u_adv_u_ext,
2440 bc_dmom_u_adv_v_ext,
2441 bc_dmom_u_adv_w_ext,
2443 bc_dmom_v_adv_u_ext,
2444 bc_dmom_v_adv_v_ext,
2445 bc_dmom_v_adv_w_ext,
2447 bc_dmom_w_adv_u_ext,
2448 bc_dmom_w_adv_v_ext,
2449 bc_dmom_w_adv_w_ext,
2450 bc_mom_uu_diff_ten_ext,
2451 bc_mom_vv_diff_ten_ext,
2452 bc_mom_ww_diff_ten_ext,
2453 bc_mom_uv_diff_ten_ext,
2454 bc_mom_uw_diff_ten_ext,
2455 bc_mom_vu_diff_ten_ext,
2456 bc_mom_vw_diff_ten_ext,
2457 bc_mom_wu_diff_ten_ext,
2458 bc_mom_wv_diff_ten_ext,
2459 bc_mom_u_source_ext,
2460 bc_mom_v_source_ext,
2461 bc_mom_w_source_ext,
2463 bc_dmom_u_ham_grad_p_ext,
2465 bc_dmom_v_ham_grad_p_ext,
2467 bc_dmom_w_ham_grad_p_ext,
2468 ebqe_dragBeam1.data()[ebNE_kb],
2469 ebqe_dragBeam2.data()[ebNE_kb],
2470 ebqe_dragBeam3.data()[ebNE_kb]);
2473 if (turbulenceClosureModel >= 3)
2475 const double turb_var_grad_0_dummy[3] = {0.,0.,0.};
2476 const double c_mu = 0.09;
2485 ebqe_vf_ext.data()[ebNE_kb],
2486 ebqe_phi_ext.data()[ebNE_kb],
2489 ebqe_turb_var_0.data()[ebNE_kb],
2490 ebqe_turb_var_1.data()[ebNE_kb],
2491 turb_var_grad_0_dummy,
2492 mom_uu_diff_ten_ext,
2493 mom_vv_diff_ten_ext,
2494 mom_ww_diff_ten_ext,
2495 mom_uv_diff_ten_ext,
2496 mom_uw_diff_ten_ext,
2497 mom_vu_diff_ten_ext,
2498 mom_vw_diff_ten_ext,
2499 mom_wu_diff_ten_ext,
2500 mom_wv_diff_ten_ext,
2513 bc_ebqe_vf_ext.data()[ebNE_kb],
2514 bc_ebqe_phi_ext.data()[ebNE_kb],
2517 ebqe_turb_var_0.data()[ebNE_kb],
2518 ebqe_turb_var_1.data()[ebNE_kb],
2519 turb_var_grad_0_dummy,
2520 bc_mom_uu_diff_ten_ext,
2521 bc_mom_vv_diff_ten_ext,
2522 bc_mom_ww_diff_ten_ext,
2523 bc_mom_uv_diff_ten_ext,
2524 bc_mom_uw_diff_ten_ext,
2525 bc_mom_vu_diff_ten_ext,
2526 bc_mom_vw_diff_ten_ext,
2527 bc_mom_wu_diff_ten_ext,
2528 bc_mom_wv_diff_ten_ext,
2529 bc_mom_u_source_ext,
2530 bc_mom_v_source_ext,
2531 bc_mom_w_source_ext);
2538 mass_adv_ext[0] -= MOVING_DOMAIN*xt_ext;
2539 mass_adv_ext[1] -= MOVING_DOMAIN*yt_ext;
2540 mass_adv_ext[2] -= MOVING_DOMAIN*zt_ext;
2542 mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext;
2543 mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext;
2544 mom_u_adv_ext[2] -= MOVING_DOMAIN*mom_u_acc_ext*zt_ext;
2545 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
2546 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
2547 dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext;
2549 mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext;
2550 mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext;
2551 mom_v_adv_ext[2] -= MOVING_DOMAIN*mom_v_acc_ext*zt_ext;
2552 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
2553 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
2554 dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext;
2556 mom_w_adv_ext[0] -= MOVING_DOMAIN*mom_w_acc_ext*xt_ext;
2557 mom_w_adv_ext[1] -= MOVING_DOMAIN*mom_w_acc_ext*yt_ext;
2558 mom_w_adv_ext[2] -= MOVING_DOMAIN*mom_w_acc_ext*zt_ext;
2559 dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext;
2560 dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext;
2561 dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext;
2564 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext;
2565 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext;
2566 bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*bc_mom_u_acc_ext*zt_ext;
2568 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext;
2569 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext;
2570 bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*bc_mom_v_acc_ext*zt_ext;
2572 bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*bc_mom_w_acc_ext*xt_ext;
2573 bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*bc_mom_w_acc_ext*yt_ext;
2574 bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*bc_mom_w_acc_ext*zt_ext;
2578 ck.calculateGScale(G,normal,h_penalty);
2579 penalty = useMetrics*C_b*h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
2581 isDOFBoundary_u.data()[ebNE_kb],
2582 isDOFBoundary_v.data()[ebNE_kb],
2583 isDOFBoundary_w.data()[ebNE_kb],
2584 isAdvectiveFluxBoundary_p.data()[ebNE_kb],
2585 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
2586 isAdvectiveFluxBoundary_v.data()[ebNE_kb],
2587 isAdvectiveFluxBoundary_w.data()[ebNE_kb],
2588 dmom_u_ham_grad_p_ext[0],
2589 bc_dmom_u_ham_grad_p_ext[0],
2596 ebqe_bc_flux_mass_ext.data()[ebNE_kb],
2597 ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
2598 ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
2599 ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
2624 &ebqe_velocity.data()[ebNE_kb_nSpace]);
2626 ebqe_phi_ext.data()[ebNE_kb],
2627 sdInfo_u_u_rowptr.data(),
2628 sdInfo_u_u_colind.data(),
2629 isDOFBoundary_u.data()[ebNE_kb],
2630 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2632 bc_mom_uu_diff_ten_ext,
2634 ebqe_bc_flux_u_diff_ext.data()[ebNE_kb],
2635 mom_uu_diff_ten_ext,
2639 flux_mom_uu_diff_ext);
2641 ebqe_phi_ext.data()[ebNE_kb],
2642 sdInfo_u_v_rowptr.data(),
2643 sdInfo_u_v_colind.data(),
2644 isDOFBoundary_v.data()[ebNE_kb],
2645 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2647 bc_mom_uv_diff_ten_ext,
2650 mom_uv_diff_ten_ext,
2654 flux_mom_uv_diff_ext);
2656 ebqe_phi_ext.data()[ebNE_kb],
2657 sdInfo_u_w_rowptr.data(),
2658 sdInfo_u_w_colind.data(),
2659 isDOFBoundary_w.data()[ebNE_kb],
2660 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2662 bc_mom_uw_diff_ten_ext,
2665 mom_uw_diff_ten_ext,
2669 flux_mom_uw_diff_ext);
2671 ebqe_phi_ext.data()[ebNE_kb],
2672 sdInfo_v_u_rowptr.data(),
2673 sdInfo_v_u_colind.data(),
2674 isDOFBoundary_u.data()[ebNE_kb],
2675 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2677 bc_mom_vu_diff_ten_ext,
2680 mom_vu_diff_ten_ext,
2684 flux_mom_vu_diff_ext);
2686 ebqe_phi_ext.data()[ebNE_kb],
2687 sdInfo_v_v_rowptr.data(),
2688 sdInfo_v_v_colind.data(),
2689 isDOFBoundary_v.data()[ebNE_kb],
2690 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2692 bc_mom_vv_diff_ten_ext,
2694 ebqe_bc_flux_v_diff_ext.data()[ebNE_kb],
2695 mom_vv_diff_ten_ext,
2699 flux_mom_vv_diff_ext);
2701 ebqe_phi_ext.data()[ebNE_kb],
2702 sdInfo_v_w_rowptr.data(),
2703 sdInfo_v_w_colind.data(),
2704 isDOFBoundary_w.data()[ebNE_kb],
2705 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2707 bc_mom_vw_diff_ten_ext,
2710 mom_vw_diff_ten_ext,
2714 flux_mom_vw_diff_ext);
2716 ebqe_phi_ext.data()[ebNE_kb],
2717 sdInfo_w_u_rowptr.data(),
2718 sdInfo_w_u_colind.data(),
2719 isDOFBoundary_u.data()[ebNE_kb],
2720 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2722 bc_mom_wu_diff_ten_ext,
2725 mom_wu_diff_ten_ext,
2729 flux_mom_wu_diff_ext);
2731 ebqe_phi_ext.data()[ebNE_kb],
2732 sdInfo_w_v_rowptr.data(),
2733 sdInfo_w_v_colind.data(),
2734 isDOFBoundary_v.data()[ebNE_kb],
2735 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2737 bc_mom_wv_diff_ten_ext,
2740 mom_wv_diff_ten_ext,
2744 flux_mom_wv_diff_ext);
2746 ebqe_phi_ext.data()[ebNE_kb],
2747 sdInfo_w_w_rowptr.data(),
2748 sdInfo_w_w_colind.data(),
2749 isDOFBoundary_w.data()[ebNE_kb],
2750 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2752 bc_mom_ww_diff_ten_ext,
2754 ebqe_bc_flux_w_diff_ext.data()[ebNE_kb],
2755 mom_ww_diff_ten_ext,
2759 flux_mom_ww_diff_ext);
2760 flux.data()[ebN*nQuadraturePoints_elementBoundary+kb] = flux_mass_ext;
2764 if (ebN < nElementBoundaries_owned)
2766 force_v_x = flux_mom_uu_diff_ext + flux_mom_uv_diff_ext + flux_mom_uw_diff_ext;
2767 force_v_y = flux_mom_vu_diff_ext + flux_mom_vv_diff_ext + flux_mom_vw_diff_ext;
2768 force_v_z = flux_mom_wu_diff_ext + flux_mom_wv_diff_ext + flux_mom_ww_diff_ext;
2770 force_p_x = p_ext*normal[0];
2771 force_p_y = p_ext*normal[1];
2772 force_p_z = p_ext*normal[2];
2774 force_x = force_p_x + force_v_x;
2775 force_y = force_p_y + force_v_y;
2776 force_z = force_p_z + force_v_z;
2778 r_x = x_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+0];
2779 r_y = y_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+1];
2780 r_z = z_ext - barycenters.data()[3*boundaryFlags.data()[ebN]+2];
2782 wettedAreas.data()[boundaryFlags.data()[ebN]] += dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2784 netForces_p.data()[3*boundaryFlags.data()[ebN]+0] += force_p_x*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2785 netForces_p.data()[3*boundaryFlags.data()[ebN]+1] += force_p_y*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2786 netForces_p.data()[3*boundaryFlags.data()[ebN]+2] += force_p_z*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2788 netForces_v.data()[3*boundaryFlags.data()[ebN]+0] += force_v_x*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2789 netForces_v.data()[3*boundaryFlags.data()[ebN]+1] += force_v_y*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2790 netForces_v.data()[3*boundaryFlags.data()[ebN]+2] += force_v_z*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2792 netMoments.data()[3*boundaryFlags.data()[ebN]+0] += (r_y*force_z - r_z*force_y)*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2793 netMoments.data()[3*boundaryFlags.data()[ebN]+1] += (r_z*force_x - r_x*force_z)*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2794 netMoments.data()[3*boundaryFlags.data()[ebN]+2] += (r_x*force_y - r_y*force_x)*dS*(1.0-ebqe_vf_ext.data()[ebNE_kb]);
2799 for (
int i=0;i<nDOF_test_element;i++)
2801 elementResidual_p[i] +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
2802 globalConservationError +=
ck.ExteriorElementBoundaryFlux(flux_mass_ext,p_test_dS[i]);
2804 elementResidual_u[i] +=
ck.ExteriorElementBoundaryFlux(flux_mom_u_adv_ext,vel_test_dS[i])+
2805 ck.ExteriorElementBoundaryFlux(flux_mom_uu_diff_ext,vel_test_dS[i])+
2806 ck.ExteriorElementBoundaryFlux(flux_mom_uv_diff_ext,vel_test_dS[i])+
2807 ck.ExteriorElementBoundaryFlux(flux_mom_uw_diff_ext,vel_test_dS[i])+
2808 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2809 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2814 sdInfo_u_u_rowptr.data(),
2815 sdInfo_u_u_colind.data(),
2816 mom_uu_diff_ten_ext,
2817 &vel_grad_test_dS[i*nSpace])+
2818 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2819 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2824 sdInfo_u_v_rowptr.data(),
2825 sdInfo_u_v_colind.data(),
2826 mom_uv_diff_ten_ext,
2827 &vel_grad_test_dS[i*nSpace])+
2828 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb],
2829 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
2834 sdInfo_u_w_rowptr.data(),
2835 sdInfo_u_w_colind.data(),
2836 mom_uw_diff_ten_ext,
2837 &vel_grad_test_dS[i*nSpace]);
2838 elementResidual_v[i] +=
ck.ExteriorElementBoundaryFlux(flux_mom_v_adv_ext,vel_test_dS[i]) +
2839 ck.ExteriorElementBoundaryFlux(flux_mom_vu_diff_ext,vel_test_dS[i])+
2840 ck.ExteriorElementBoundaryFlux(flux_mom_vv_diff_ext,vel_test_dS[i])+
2841 ck.ExteriorElementBoundaryFlux(flux_mom_vw_diff_ext,vel_test_dS[i])+
2842 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2843 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2848 sdInfo_v_u_rowptr.data(),
2849 sdInfo_v_u_colind.data(),
2850 mom_vu_diff_ten_ext,
2851 &vel_grad_test_dS[i*nSpace])+
2852 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2853 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2858 sdInfo_v_v_rowptr.data(),
2859 sdInfo_v_v_colind.data(),
2860 mom_vv_diff_ten_ext,
2861 &vel_grad_test_dS[i*nSpace])+
2862 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb],
2863 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
2868 sdInfo_v_w_rowptr.data(),
2869 sdInfo_v_w_colind.data(),
2870 mom_vw_diff_ten_ext,
2871 &vel_grad_test_dS[i*nSpace]);
2873 elementResidual_w[i] +=
ck.ExteriorElementBoundaryFlux(flux_mom_w_adv_ext,vel_test_dS[i]) +
2874 ck.ExteriorElementBoundaryFlux(flux_mom_wu_diff_ext,vel_test_dS[i])+
2875 ck.ExteriorElementBoundaryFlux(flux_mom_wv_diff_ext,vel_test_dS[i])+
2876 ck.ExteriorElementBoundaryFlux(flux_mom_ww_diff_ext,vel_test_dS[i])+
2877 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_u.data()[ebNE_kb],
2878 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2883 sdInfo_w_u_rowptr.data(),
2884 sdInfo_w_u_colind.data(),
2885 mom_wu_diff_ten_ext,
2886 &vel_grad_test_dS[i*nSpace])+
2887 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_v.data()[ebNE_kb],
2888 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2893 sdInfo_w_v_rowptr.data(),
2894 sdInfo_w_v_colind.data(),
2895 mom_wv_diff_ten_ext,
2896 &vel_grad_test_dS[i*nSpace])+
2897 ck.ExteriorElementBoundaryDiffusionAdjoint(isDOFBoundary_w.data()[ebNE_kb],
2898 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
2903 sdInfo_w_w_rowptr.data(),
2904 sdInfo_w_w_colind.data(),
2905 mom_ww_diff_ten_ext,
2906 &vel_grad_test_dS[i*nSpace]);
2912 for (
int i=0;i<nDOF_test_element;i++)
2914 int eN_i = eN*nDOF_test_element+i;
2916 elementResidual_p_save.data()[eN_i] += elementResidual_p[i];
2918 globalResidual.data()[offset_p+stride_p*p_l2g.data()[eN_i]]+=elementResidual_p[i];
2919 globalResidual.data()[offset_u+stride_u*vel_l2g.data()[eN_i]]+=elementResidual_u[i];
2920 globalResidual.data()[offset_v+stride_v*vel_l2g.data()[eN_i]]+=elementResidual_v[i];
2921 globalResidual.data()[offset_w+stride_w*vel_l2g.data()[eN_i]]+=elementResidual_w[i];
2928 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
2929 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
2930 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
2931 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
2932 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
2933 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
2934 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
2935 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
2936 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
2937 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
2938 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
2939 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
2940 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
2941 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
2942 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
2943 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
2944 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
2945 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
2946 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
2947 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
2948 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
2949 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
2950 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
2951 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
2952 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
2953 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
2954 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
2955 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
2956 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
2957 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
2958 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
2959 double hFactor = args.
scalar<
double>(
"hFactor");
2960 int nElements_global = args.
scalar<
int>(
"nElements_global");
2961 double useRBLES = args.
scalar<
double>(
"useRBLES");
2962 double useMetrics = args.
scalar<
double>(
"useMetrics");
2963 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
2964 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
2965 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
2966 double sigma = args.
scalar<
double>(
"sigma");
2971 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
2972 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
2973 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
2974 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
2975 double C_dg = args.
scalar<
double>(
"C_dg");
2976 double C_b = args.
scalar<
double>(
"C_b");
2977 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
2978 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
2979 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
2980 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
2981 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
2982 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
2983 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
2984 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
2985 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
2986 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
2987 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
2988 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
2989 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
2990 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
2991 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
2992 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
2993 xt::pyarray<double>& g = args.
array<
double>(
"g");
2994 const double useVF = args.
scalar<
double>(
"useVF");
2995 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
2996 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
2997 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
2998 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
2999 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
3000 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
3001 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
3002 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
3003 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
3004 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
3005 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
3006 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
3007 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
3008 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
3009 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
3010 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
3011 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
3012 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
3013 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
3014 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
3015 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
3016 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
3017 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
3018 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
3019 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
3020 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
3021 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
3022 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
3023 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
3024 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
3025 xt::pyarray<int>& csrRowIndeces_p_p = args.
array<
int>(
"csrRowIndeces_p_p");
3026 xt::pyarray<int>& csrColumnOffsets_p_p = args.
array<
int>(
"csrColumnOffsets_p_p");
3027 xt::pyarray<int>& csrRowIndeces_p_u = args.
array<
int>(
"csrRowIndeces_p_u");
3028 xt::pyarray<int>& csrColumnOffsets_p_u = args.
array<
int>(
"csrColumnOffsets_p_u");
3029 xt::pyarray<int>& csrRowIndeces_p_v = args.
array<
int>(
"csrRowIndeces_p_v");
3030 xt::pyarray<int>& csrColumnOffsets_p_v = args.
array<
int>(
"csrColumnOffsets_p_v");
3031 xt::pyarray<int>& csrRowIndeces_p_w = args.
array<
int>(
"csrRowIndeces_p_w");
3032 xt::pyarray<int>& csrColumnOffsets_p_w = args.
array<
int>(
"csrColumnOffsets_p_w");
3033 xt::pyarray<int>& csrRowIndeces_u_p = args.
array<
int>(
"csrRowIndeces_u_p");
3034 xt::pyarray<int>& csrColumnOffsets_u_p = args.
array<
int>(
"csrColumnOffsets_u_p");
3035 xt::pyarray<int>& csrRowIndeces_u_u = args.
array<
int>(
"csrRowIndeces_u_u");
3036 xt::pyarray<int>& csrColumnOffsets_u_u = args.
array<
int>(
"csrColumnOffsets_u_u");
3037 xt::pyarray<int>& csrRowIndeces_u_v = args.
array<
int>(
"csrRowIndeces_u_v");
3038 xt::pyarray<int>& csrColumnOffsets_u_v = args.
array<
int>(
"csrColumnOffsets_u_v");
3039 xt::pyarray<int>& csrRowIndeces_u_w = args.
array<
int>(
"csrRowIndeces_u_w");
3040 xt::pyarray<int>& csrColumnOffsets_u_w = args.
array<
int>(
"csrColumnOffsets_u_w");
3041 xt::pyarray<int>& csrRowIndeces_v_p = args.
array<
int>(
"csrRowIndeces_v_p");
3042 xt::pyarray<int>& csrColumnOffsets_v_p = args.
array<
int>(
"csrColumnOffsets_v_p");
3043 xt::pyarray<int>& csrRowIndeces_v_u = args.
array<
int>(
"csrRowIndeces_v_u");
3044 xt::pyarray<int>& csrColumnOffsets_v_u = args.
array<
int>(
"csrColumnOffsets_v_u");
3045 xt::pyarray<int>& csrRowIndeces_v_v = args.
array<
int>(
"csrRowIndeces_v_v");
3046 xt::pyarray<int>& csrColumnOffsets_v_v = args.
array<
int>(
"csrColumnOffsets_v_v");
3047 xt::pyarray<int>& csrRowIndeces_v_w = args.
array<
int>(
"csrRowIndeces_v_w");
3048 xt::pyarray<int>& csrColumnOffsets_v_w = args.
array<
int>(
"csrColumnOffsets_v_w");
3049 xt::pyarray<int>& csrRowIndeces_w_p = args.
array<
int>(
"csrRowIndeces_w_p");
3050 xt::pyarray<int>& csrColumnOffsets_w_p = args.
array<
int>(
"csrColumnOffsets_w_p");
3051 xt::pyarray<int>& csrRowIndeces_w_u = args.
array<
int>(
"csrRowIndeces_w_u");
3052 xt::pyarray<int>& csrColumnOffsets_w_u = args.
array<
int>(
"csrColumnOffsets_w_u");
3053 xt::pyarray<int>& csrRowIndeces_w_v = args.
array<
int>(
"csrRowIndeces_w_v");
3054 xt::pyarray<int>& csrColumnOffsets_w_v = args.
array<
int>(
"csrColumnOffsets_w_v");
3055 xt::pyarray<int>& csrRowIndeces_w_w = args.
array<
int>(
"csrRowIndeces_w_w");
3056 xt::pyarray<int>& csrColumnOffsets_w_w = args.
array<
int>(
"csrColumnOffsets_w_w");
3057 xt::pyarray<double>& globalJacobian = args.
array<
double>(
"globalJacobian");
3058 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
3059 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
3060 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
3061 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
3062 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
3063 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
3064 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
3065 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
3066 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
3067 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
3068 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
3069 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
3070 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
3071 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
3072 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
3073 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
3074 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
3075 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
3076 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
3077 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
3078 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
3079 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
3080 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
3081 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
3082 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
3083 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
3084 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
3085 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
3086 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
3087 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
3088 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
3089 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
3090 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
3091 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
3092 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
3093 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
3094 xt::pyarray<int>& csrColumnOffsets_eb_p_p = args.
array<
int>(
"csrColumnOffsets_eb_p_p");
3095 xt::pyarray<int>& csrColumnOffsets_eb_p_u = args.
array<
int>(
"csrColumnOffsets_eb_p_u");
3096 xt::pyarray<int>& csrColumnOffsets_eb_p_v = args.
array<
int>(
"csrColumnOffsets_eb_p_v");
3097 xt::pyarray<int>& csrColumnOffsets_eb_p_w = args.
array<
int>(
"csrColumnOffsets_eb_p_w");
3098 xt::pyarray<int>& csrColumnOffsets_eb_u_p = args.
array<
int>(
"csrColumnOffsets_eb_u_p");
3099 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.
array<
int>(
"csrColumnOffsets_eb_u_u");
3100 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.
array<
int>(
"csrColumnOffsets_eb_u_v");
3101 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.
array<
int>(
"csrColumnOffsets_eb_u_w");
3102 xt::pyarray<int>& csrColumnOffsets_eb_v_p = args.
array<
int>(
"csrColumnOffsets_eb_v_p");
3103 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.
array<
int>(
"csrColumnOffsets_eb_v_u");
3104 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.
array<
int>(
"csrColumnOffsets_eb_v_v");
3105 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.
array<
int>(
"csrColumnOffsets_eb_v_w");
3106 xt::pyarray<int>& csrColumnOffsets_eb_w_p = args.
array<
int>(
"csrColumnOffsets_eb_w_p");
3107 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.
array<
int>(
"csrColumnOffsets_eb_w_u");
3108 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.
array<
int>(
"csrColumnOffsets_eb_w_v");
3109 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.
array<
int>(
"csrColumnOffsets_eb_w_w");
3110 xt::pyarray<double>& q_dragBeam1 = args.
array<
double>(
"q_dragBeam1");
3111 xt::pyarray<double>& q_dragBeam2 = args.
array<
double>(
"q_dragBeam2");
3112 xt::pyarray<double>& q_dragBeam3 = args.
array<
double>(
"q_dragBeam3");
3113 xt::pyarray<double>& ebqe_dragBeam1 = args.
array<
double>(
"ebqe_dragBeam1");
3114 xt::pyarray<double>& ebqe_dragBeam2 = args.
array<
double>(
"ebqe_dragBeam2");
3115 xt::pyarray<double>& ebqe_dragBeam3 = args.
array<
double>(
"ebqe_dragBeam3");
3119 for(
int eN=0;eN<nElements_global;eN++)
3121 double eps_rho,eps_mu;
3123 double elementJacobian_p_p[nDOF_test_element][nDOF_trial_element],
3124 elementJacobian_p_u[nDOF_test_element][nDOF_trial_element],
3125 elementJacobian_p_v[nDOF_test_element][nDOF_trial_element],
3126 elementJacobian_p_w[nDOF_test_element][nDOF_trial_element],
3127 elementJacobian_u_p[nDOF_test_element][nDOF_trial_element],
3128 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
3129 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
3130 elementJacobian_u_w[nDOF_test_element][nDOF_trial_element],
3131 elementJacobian_v_p[nDOF_test_element][nDOF_trial_element],
3132 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
3133 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element],
3134 elementJacobian_v_w[nDOF_test_element][nDOF_trial_element],
3135 elementJacobian_w_p[nDOF_test_element][nDOF_trial_element],
3136 elementJacobian_w_u[nDOF_test_element][nDOF_trial_element],
3137 elementJacobian_w_v[nDOF_test_element][nDOF_trial_element],
3138 elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
3139 for (
int i=0;i<nDOF_test_element;i++)
3140 for (
int j=0;j<nDOF_trial_element;j++)
3142 elementJacobian_p_p[i][j]=0.0;
3143 elementJacobian_p_u[i][j]=0.0;
3144 elementJacobian_p_v[i][j]=0.0;
3145 elementJacobian_p_w[i][j]=0.0;
3146 elementJacobian_u_p[i][j]=0.0;
3147 elementJacobian_u_u[i][j]=0.0;
3148 elementJacobian_u_v[i][j]=0.0;
3149 elementJacobian_u_w[i][j]=0.0;
3150 elementJacobian_v_p[i][j]=0.0;
3151 elementJacobian_v_u[i][j]=0.0;
3152 elementJacobian_v_v[i][j]=0.0;
3153 elementJacobian_v_w[i][j]=0.0;
3154 elementJacobian_w_p[i][j]=0.0;
3155 elementJacobian_w_u[i][j]=0.0;
3156 elementJacobian_w_v[i][j]=0.0;
3157 elementJacobian_w_w[i][j]=0.0;
3159 for (
int k=0;k<nQuadraturePoints_element;k++)
3161 int eN_k = eN*nQuadraturePoints_element+k,
3162 eN_k_nSpace = eN_k*nSpace,
3163 eN_nDOF_trial_element = eN*nDOF_trial_element;
3166 double p=0.0,
u=0.0,
v=0.0,
w=0.0,
3167 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
3175 dmass_adv_u[nSpace],
3176 dmass_adv_v[nSpace],
3177 dmass_adv_w[nSpace],
3179 dmom_u_adv_u[nSpace],
3180 dmom_u_adv_v[nSpace],
3181 dmom_u_adv_w[nSpace],
3183 dmom_v_adv_u[nSpace],
3184 dmom_v_adv_v[nSpace],
3185 dmom_v_adv_w[nSpace],
3187 dmom_w_adv_u[nSpace],
3188 dmom_w_adv_v[nSpace],
3189 dmom_w_adv_w[nSpace],
3190 mom_uu_diff_ten[nSpace],
3191 mom_vv_diff_ten[nSpace],
3192 mom_ww_diff_ten[nSpace],
3203 dmom_u_ham_grad_p[nSpace],
3205 dmom_v_ham_grad_p[nSpace],
3207 dmom_w_ham_grad_p[nSpace],
3218 dpdeResidual_p_u[nDOF_trial_element],dpdeResidual_p_v[nDOF_trial_element],dpdeResidual_p_w[nDOF_trial_element],
3219 dpdeResidual_u_p[nDOF_trial_element],dpdeResidual_u_u[nDOF_trial_element],
3220 dpdeResidual_v_p[nDOF_trial_element],dpdeResidual_v_v[nDOF_trial_element],
3221 dpdeResidual_w_p[nDOF_trial_element],dpdeResidual_w_w[nDOF_trial_element],
3222 Lstar_u_p[nDOF_test_element],
3223 Lstar_v_p[nDOF_test_element],
3224 Lstar_w_p[nDOF_test_element],
3225 Lstar_u_u[nDOF_test_element],
3226 Lstar_v_v[nDOF_test_element],
3227 Lstar_w_w[nDOF_test_element],
3228 Lstar_p_u[nDOF_test_element],
3229 Lstar_p_v[nDOF_test_element],
3230 Lstar_p_w[nDOF_test_element],
3235 dsubgridError_p_u[nDOF_trial_element],
3236 dsubgridError_p_v[nDOF_trial_element],
3237 dsubgridError_p_w[nDOF_trial_element],
3238 dsubgridError_u_p[nDOF_trial_element],
3239 dsubgridError_u_u[nDOF_trial_element],
3240 dsubgridError_v_p[nDOF_trial_element],
3241 dsubgridError_v_v[nDOF_trial_element],
3242 dsubgridError_w_p[nDOF_trial_element],
3243 dsubgridError_w_w[nDOF_trial_element],
3244 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
3245 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
3248 jacInv[nSpace*nSpace],
3249 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
3251 p_test_dV[nDOF_test_element],vel_test_dV[nDOF_test_element],
3252 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
3257 dmom_u_source[nSpace],
3258 dmom_v_source[nSpace],
3259 dmom_w_source[nSpace],
3262 G[nSpace*nSpace],G_dd_G,tr_G,h_phi, dmom_adv_star[nSpace], dmom_adv_sge[nSpace];
3264 ck.calculateMapping_element(eN,
3268 mesh_trial_ref.data(),
3269 mesh_grad_trial_ref.data(),
3274 ck.calculateH_element(eN,
3276 nodeDiametersArray.data(),
3278 mesh_trial_ref.data(),
3280 ck.calculateMappingVelocity_element(eN,
3282 mesh_velocity_dof.data(),
3284 mesh_trial_ref.data(),
3289 dV = fabs(jacDet)*dV_ref.data()[k];
3290 ck.calculateG(jacInv,G,G_dd_G,tr_G);
3293 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3294 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3297 ck.gradTrialFromRef(&p_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,p_grad_trial);
3298 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
3300 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_ref.data()[k*nDOF_trial_element],p);
3301 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
u);
3302 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
v);
3303 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
w);
3305 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial,grad_p);
3306 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_u);
3307 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_v);
3308 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial,grad_w);
3310 for (
int j=0;j<nDOF_trial_element;j++)
3312 p_test_dV[j] = p_test_ref.data()[k*nDOF_trial_element+j]*dV;
3313 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
3314 for (
int I=0;I<nSpace;I++)
3316 p_grad_test_dV[j*nSpace+I] = p_grad_trial[j*nSpace+I]*dV;
3317 vel_grad_test_dV[j*nSpace+I] = vel_grad_trial[j*nSpace+I]*dV;
3322 porosity = q_porosity.data()[eN_k];
3334 elementDiameter.data()[eN],
3335 smagorinskyConstant,
3336 turbulenceClosureModel,
3341 &normal_phi.data()[eN_k_nSpace],
3342 kappa_phi.data()[eN_k],
3394 q_dragBeam1.data()[eN_k],
3395 q_dragBeam2.data()[eN_k],
3396 q_dragBeam3.data()[eN_k]);
3398 mass_source = q_mass_source.data()[eN_k];
3404 q_dragAlpha.data()[eN_k],
3405 q_dragBeta.data()[eN_k],
3418 eps_solid.data()[0],
3419 phi_solid.data()[eN_k],
3420 q_velocity_solid.data()[eN_k_nSpace+0],
3421 q_velocity_solid.data()[eN_k_nSpace+1],
3422 q_velocity_solid.data()[eN_k_nSpace+2],
3430 if (turbulenceClosureModel >= 3)
3432 const double c_mu = 0.09;
3445 q_turb_var_0.data()[eN_k],
3446 q_turb_var_1.data()[eN_k],
3447 &q_turb_var_grad_0.data()[eN_k_nSpace],
3466 mass_adv[0] -= MOVING_DOMAIN*
xt;
3467 mass_adv[1] -= MOVING_DOMAIN*yt;
3468 mass_adv[2] -= MOVING_DOMAIN*zt;
3470 mom_u_adv[0] -= MOVING_DOMAIN*mom_u_acc*
xt;
3471 mom_u_adv[1] -= MOVING_DOMAIN*mom_u_acc*yt;
3472 mom_u_adv[2] -= MOVING_DOMAIN*mom_u_acc*zt;
3473 dmom_u_adv_u[0] -= MOVING_DOMAIN*dmom_u_acc_u*
xt;
3474 dmom_u_adv_u[1] -= MOVING_DOMAIN*dmom_u_acc_u*yt;
3475 dmom_u_adv_u[2] -= MOVING_DOMAIN*dmom_u_acc_u*zt;
3477 mom_v_adv[0] -= MOVING_DOMAIN*mom_v_acc*
xt;
3478 mom_v_adv[1] -= MOVING_DOMAIN*mom_v_acc*yt;
3479 mom_v_adv[2] -= MOVING_DOMAIN*mom_v_acc*zt;
3480 dmom_v_adv_v[0] -= MOVING_DOMAIN*dmom_v_acc_v*
xt;
3481 dmom_v_adv_v[1] -= MOVING_DOMAIN*dmom_v_acc_v*yt;
3482 dmom_v_adv_v[2] -= MOVING_DOMAIN*dmom_v_acc_v*zt;
3484 mom_w_adv[0] -= MOVING_DOMAIN*mom_w_acc*
xt;
3485 mom_w_adv[1] -= MOVING_DOMAIN*mom_w_acc*yt;
3486 mom_w_adv[2] -= MOVING_DOMAIN*mom_w_acc*zt;
3487 dmom_w_adv_w[0] -= MOVING_DOMAIN*dmom_w_acc_w*
xt;
3488 dmom_w_adv_w[1] -= MOVING_DOMAIN*dmom_w_acc_w*yt;
3489 dmom_w_adv_w[2] -= MOVING_DOMAIN*dmom_w_acc_w*zt;
3494 q_mom_u_acc_beta_bdf.data()[eN_k],
3500 q_mom_v_acc_beta_bdf.data()[eN_k],
3506 q_mom_w_acc_beta_bdf.data()[eN_k],
3514 dmom_adv_sge[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt);
3515 dmom_adv_sge[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt);
3516 dmom_adv_sge[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt);
3520 pdeResidual_p =
ck.Advection_strong(dmass_adv_u,grad_u) +
3521 ck.Advection_strong(dmass_adv_v,grad_v) +
3522 ck.Advection_strong(dmass_adv_w,grad_w);
3524 pdeResidual_u =
ck.Mass_strong(mom_u_acc_t) +
3525 ck.Advection_strong(dmom_adv_sge,grad_u) +
3526 ck.Hamiltonian_strong(dmom_u_ham_grad_p,grad_p) +
3527 ck.Reaction_strong(mom_u_source);
3529 pdeResidual_v =
ck.Mass_strong(mom_v_acc_t) +
3530 ck.Advection_strong(dmom_adv_sge,grad_v) +
3531 ck.Hamiltonian_strong(dmom_v_ham_grad_p,grad_p) +
3532 ck.Reaction_strong(mom_v_source);
3534 pdeResidual_w =
ck.Mass_strong(mom_w_acc_t) +
3535 ck.Advection_strong(dmom_adv_sge,grad_w) +
3536 ck.Hamiltonian_strong(dmom_w_ham_grad_p,grad_p) +
3537 ck.Reaction_strong(mom_w_source);
3540 for (
int j=0;j<nDOF_trial_element;j++)
3542 int j_nSpace = j*nSpace;
3543 dpdeResidual_p_u[j]=
ck.AdvectionJacobian_strong(dmass_adv_u,&vel_grad_trial[j_nSpace]);
3544 dpdeResidual_p_v[j]=
ck.AdvectionJacobian_strong(dmass_adv_v,&vel_grad_trial[j_nSpace]);
3545 dpdeResidual_p_w[j]=
ck.AdvectionJacobian_strong(dmass_adv_w,&vel_grad_trial[j_nSpace]);
3547 dpdeResidual_u_p[j]=
ck.HamiltonianJacobian_strong(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace]);
3548 dpdeResidual_u_u[j]=
ck.MassJacobian_strong(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3549 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]);
3551 dpdeResidual_v_p[j]=
ck.HamiltonianJacobian_strong(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace]);
3552 dpdeResidual_v_v[j]=
ck.MassJacobian_strong(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3553 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]);
3555 dpdeResidual_w_p[j]=
ck.HamiltonianJacobian_strong(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace]);
3556 dpdeResidual_w_w[j]=
ck.MassJacobian_strong(dmom_w_acc_w_t,vel_trial_ref.data()[k*nDOF_trial_element+j]) +
3557 ck.AdvectionJacobian_strong(dmom_adv_sge,&vel_grad_trial[j_nSpace]);
3560 dpdeResidual_u_u[j]+=
ck.ReactionJacobian_strong(dmom_u_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3561 dpdeResidual_v_v[j]+=
ck.ReactionJacobian_strong(dmom_v_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3562 dpdeResidual_w_w[j]+=
ck.ReactionJacobian_strong(dmom_w_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j]);
3567 double tmpR=dmom_u_acc_u_t + dmom_u_source[0];
3569 elementDiameter.data()[eN],
3574 dmom_u_ham_grad_p[0],
3577 q_cfl.data()[eN_k]);
3584 dmom_u_ham_grad_p[0],
3587 q_cfl.data()[eN_k]);
3590 tau_v = useMetrics*tau_v1+(1.0-useMetrics)*tau_v0;
3591 tau_p = useMetrics*tau_p1+(1.0-useMetrics)*tau_p0;
3625 dmom_adv_star[0] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+0] - MOVING_DOMAIN*
xt + useRBLES*subgridError_u);
3626 dmom_adv_star[1] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+1] - MOVING_DOMAIN*yt + useRBLES*subgridError_v);
3627 dmom_adv_star[2] = dmom_u_acc_u*(q_velocity_sge.data()[eN_k_nSpace+2] - MOVING_DOMAIN*zt + useRBLES*subgridError_w);
3630 for (
int i=0;i<nDOF_test_element;i++)
3632 int i_nSpace = i*nSpace;
3633 Lstar_u_p[i]=
ck.Advection_adjoint(dmass_adv_u,&p_grad_test_dV[i_nSpace]);
3634 Lstar_v_p[i]=
ck.Advection_adjoint(dmass_adv_v,&p_grad_test_dV[i_nSpace]);
3635 Lstar_w_p[i]=
ck.Advection_adjoint(dmass_adv_w,&p_grad_test_dV[i_nSpace]);
3636 Lstar_u_u[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3637 Lstar_v_v[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3638 Lstar_w_w[i]=
ck.Advection_adjoint(dmom_adv_star,&vel_grad_test_dV[i_nSpace]);
3639 Lstar_p_u[i]=
ck.Hamiltonian_adjoint(dmom_u_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
3640 Lstar_p_v[i]=
ck.Hamiltonian_adjoint(dmom_v_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
3641 Lstar_p_w[i]=
ck.Hamiltonian_adjoint(dmom_w_ham_grad_p,&vel_grad_test_dV[i_nSpace]);
3643 Lstar_u_u[i]+=
ck.Reaction_adjoint(dmom_u_source[0],vel_test_dV[i]);
3644 Lstar_v_v[i]+=
ck.Reaction_adjoint(dmom_v_source[1],vel_test_dV[i]);
3645 Lstar_w_w[i]+=
ck.Reaction_adjoint(dmom_w_source[2],vel_test_dV[i]);
3649 dmom_u_adv_u[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3650 dmom_u_adv_u[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3651 dmom_u_adv_u[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
3653 dmom_v_adv_v[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3654 dmom_v_adv_v[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3655 dmom_v_adv_v[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
3657 dmom_w_adv_w[0] += dmom_u_acc_u*(useRBLES*subgridError_u);
3658 dmom_w_adv_w[1] += dmom_u_acc_u*(useRBLES*subgridError_v);
3659 dmom_w_adv_w[2] += dmom_u_acc_u*(useRBLES*subgridError_w);
3663 for(
int i=0;i<nDOF_test_element;i++)
3665 int i_nSpace = i*nSpace;
3666 for(
int j=0;j<nDOF_trial_element;j++)
3668 int j_nSpace = j*nSpace;
3669 elementJacobian_p_p[i][j] +=
ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_p[i]) +
3670 ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_p[i]) +
3671 ck.SubgridErrorJacobian(dsubgridError_w_p[j],Lstar_w_p[i]);
3673 elementJacobian_p_u[i][j] +=
ck.AdvectionJacobian_weak(dmass_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) +
3674 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_p[i]);
3675 elementJacobian_p_v[i][j] +=
ck.AdvectionJacobian_weak(dmass_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) +
3676 ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_p[i]);
3677 elementJacobian_p_w[i][j] +=
ck.AdvectionJacobian_weak(dmass_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&p_grad_test_dV[i_nSpace]) +
3678 ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_p[i]);
3680 elementJacobian_u_p[i][j] +=
ck.HamiltonianJacobian_weak(dmom_u_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) +
3681 ck.SubgridErrorJacobian(dsubgridError_u_p[j],Lstar_u_u[i]);
3682 elementJacobian_u_u[i][j] +=
ck.MassJacobian_weak(dmom_u_acc_u_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3683 ck.AdvectionJacobian_weak(dmom_u_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3684 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]) +
3686 ck.ReactionJacobian_weak(dmom_u_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3688 ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_u[i]) +
3689 ck.SubgridErrorJacobian(dsubgridError_u_u[j],Lstar_u_u[i]) +
3690 ck.NumericalDiffusionJacobian(q_numDiff_u_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3691 elementJacobian_u_v[i][j] +=
ck.AdvectionJacobian_weak(dmom_u_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3692 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]) +
3694 ck.ReactionJacobian_weak(dmom_u_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3696 ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_u[i]);
3697 elementJacobian_u_w[i][j] +=
ck.AdvectionJacobian_weak(dmom_u_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3698 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]) +
3700 ck.ReactionJacobian_weak(dmom_u_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3702 ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_u[i]);
3704 elementJacobian_v_p[i][j] +=
ck.HamiltonianJacobian_weak(dmom_v_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) +
3705 ck.SubgridErrorJacobian(dsubgridError_v_p[j],Lstar_v_v[i]);
3706 elementJacobian_v_u[i][j] +=
ck.AdvectionJacobian_weak(dmom_v_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3707 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]) +
3709 ck.ReactionJacobian_weak(dmom_v_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3711 ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_v[i]);
3712 elementJacobian_v_v[i][j] +=
ck.MassJacobian_weak(dmom_v_acc_v_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3713 ck.AdvectionJacobian_weak(dmom_v_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3714 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]) +
3716 ck.ReactionJacobian_weak(dmom_v_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3718 ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_v[i]) +
3719 ck.SubgridErrorJacobian(dsubgridError_v_v[j],Lstar_v_v[i]) +
3720 ck.NumericalDiffusionJacobian(q_numDiff_v_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3721 elementJacobian_v_w[i][j] +=
ck.AdvectionJacobian_weak(dmom_v_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3722 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]) +
3724 ck.ReactionJacobian_weak(dmom_v_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3726 ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_v[i]);
3728 elementJacobian_w_p[i][j] +=
ck.HamiltonianJacobian_weak(dmom_w_ham_grad_p,&p_grad_trial[j_nSpace],vel_test_dV[i]) +
3729 ck.SubgridErrorJacobian(dsubgridError_w_p[j],Lstar_w_w[i]);
3730 elementJacobian_w_u[i][j] +=
ck.AdvectionJacobian_weak(dmom_w_adv_u,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3731 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]) +
3733 ck.ReactionJacobian_weak(dmom_w_source[0],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3735 ck.SubgridErrorJacobian(dsubgridError_p_u[j],Lstar_p_w[i]);
3736 elementJacobian_w_v[i][j] +=
ck.AdvectionJacobian_weak(dmom_w_adv_v,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3737 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]) +
3739 ck.ReactionJacobian_weak(dmom_w_source[1],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3741 ck.SubgridErrorJacobian(dsubgridError_p_v[j],Lstar_p_w[i]);
3742 elementJacobian_w_w[i][j] +=
ck.MassJacobian_weak(dmom_w_acc_w_t,vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3743 ck.AdvectionJacobian_weak(dmom_w_adv_w,vel_trial_ref.data()[k*nDOF_trial_element+j],&vel_grad_test_dV[i_nSpace]) +
3744 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]) +
3746 ck.ReactionJacobian_weak(dmom_w_source[2],vel_trial_ref.data()[k*nDOF_trial_element+j],vel_test_dV[i]) +
3748 ck.SubgridErrorJacobian(dsubgridError_p_w[j],Lstar_p_w[i]) +
3749 ck.SubgridErrorJacobian(dsubgridError_w_w[j],Lstar_w_w[i]) +
3750 ck.NumericalDiffusionJacobian(q_numDiff_w_last.data()[eN_k],&vel_grad_trial[j_nSpace],&vel_grad_test_dV[i_nSpace]);
3757 for (
int i=0;i<nDOF_test_element;i++)
3759 int eN_i = eN*nDOF_test_element+i;
3760 for (
int j=0;j<nDOF_trial_element;j++)
3762 int eN_i_j = eN_i*nDOF_trial_element+j;
3763 globalJacobian.data()[csrRowIndeces_p_p.data()[eN_i] + csrColumnOffsets_p_p.data()[eN_i_j]] += elementJacobian_p_p[i][j];
3764 globalJacobian.data()[csrRowIndeces_p_u.data()[eN_i] + csrColumnOffsets_p_u.data()[eN_i_j]] += elementJacobian_p_u[i][j];
3765 globalJacobian.data()[csrRowIndeces_p_v.data()[eN_i] + csrColumnOffsets_p_v.data()[eN_i_j]] += elementJacobian_p_v[i][j];
3766 globalJacobian.data()[csrRowIndeces_p_w.data()[eN_i] + csrColumnOffsets_p_w.data()[eN_i_j]] += elementJacobian_p_w[i][j];
3768 globalJacobian.data()[csrRowIndeces_u_p.data()[eN_i] + csrColumnOffsets_u_p.data()[eN_i_j]] += elementJacobian_u_p[i][j];
3769 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
3770 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_u_v.data()[eN_i_j]] += elementJacobian_u_v[i][j];
3771 globalJacobian.data()[csrRowIndeces_u_w.data()[eN_i] + csrColumnOffsets_u_w.data()[eN_i_j]] += elementJacobian_u_w[i][j];
3773 globalJacobian.data()[csrRowIndeces_v_p.data()[eN_i] + csrColumnOffsets_v_p.data()[eN_i_j]] += elementJacobian_v_p[i][j];
3774 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_v_u.data()[eN_i_j]] += elementJacobian_v_u[i][j];
3775 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_v_v.data()[eN_i_j]] += elementJacobian_v_v[i][j];
3776 globalJacobian.data()[csrRowIndeces_v_w.data()[eN_i] + csrColumnOffsets_v_w.data()[eN_i_j]] += elementJacobian_v_w[i][j];
3778 globalJacobian.data()[csrRowIndeces_w_p.data()[eN_i] + csrColumnOffsets_w_p.data()[eN_i_j]] += elementJacobian_w_p[i][j];
3779 globalJacobian.data()[csrRowIndeces_w_u.data()[eN_i] + csrColumnOffsets_w_u.data()[eN_i_j]] += elementJacobian_w_u[i][j];
3780 globalJacobian.data()[csrRowIndeces_w_v.data()[eN_i] + csrColumnOffsets_w_v.data()[eN_i_j]] += elementJacobian_w_v[i][j];
3781 globalJacobian.data()[csrRowIndeces_w_w.data()[eN_i] + csrColumnOffsets_w_w.data()[eN_i_j]] += elementJacobian_w_w[i][j];
3788 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
3790 int ebN = exteriorElementBoundariesArray.data()[ebNE],
3791 eN = elementBoundaryElementsArray.data()[ebN*2+0],
3792 eN_nDOF_trial_element = eN*nDOF_trial_element,
3793 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
3794 double eps_rho,eps_mu;
3795 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
3797 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
3798 ebNE_kb_nSpace = ebNE_kb*nSpace,
3799 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
3800 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
3811 dmom_u_acc_u_ext=0.0,
3813 dmom_v_acc_v_ext=0.0,
3815 dmom_w_acc_w_ext=0.0,
3816 mass_adv_ext[nSpace],
3817 dmass_adv_u_ext[nSpace],
3818 dmass_adv_v_ext[nSpace],
3819 dmass_adv_w_ext[nSpace],
3820 mom_u_adv_ext[nSpace],
3821 dmom_u_adv_u_ext[nSpace],
3822 dmom_u_adv_v_ext[nSpace],
3823 dmom_u_adv_w_ext[nSpace],
3824 mom_v_adv_ext[nSpace],
3825 dmom_v_adv_u_ext[nSpace],
3826 dmom_v_adv_v_ext[nSpace],
3827 dmom_v_adv_w_ext[nSpace],
3828 mom_w_adv_ext[nSpace],
3829 dmom_w_adv_u_ext[nSpace],
3830 dmom_w_adv_v_ext[nSpace],
3831 dmom_w_adv_w_ext[nSpace],
3832 mom_uu_diff_ten_ext[nSpace],
3833 mom_vv_diff_ten_ext[nSpace],
3834 mom_ww_diff_ten_ext[nSpace],
3835 mom_uv_diff_ten_ext[1],
3836 mom_uw_diff_ten_ext[1],
3837 mom_vu_diff_ten_ext[1],
3838 mom_vw_diff_ten_ext[1],
3839 mom_wu_diff_ten_ext[1],
3840 mom_wv_diff_ten_ext[1],
3841 mom_u_source_ext=0.0,
3842 mom_v_source_ext=0.0,
3843 mom_w_source_ext=0.0,
3845 dmom_u_ham_grad_p_ext[nSpace],
3847 dmom_v_ham_grad_p_ext[nSpace],
3849 dmom_w_ham_grad_p_ext[nSpace],
3850 dmom_u_adv_p_ext[nSpace],
3851 dmom_v_adv_p_ext[nSpace],
3852 dmom_w_adv_p_ext[nSpace],
3853 dflux_mass_u_ext=0.0,
3854 dflux_mass_v_ext=0.0,
3855 dflux_mass_w_ext=0.0,
3856 dflux_mom_u_adv_p_ext=0.0,
3857 dflux_mom_u_adv_u_ext=0.0,
3858 dflux_mom_u_adv_v_ext=0.0,
3859 dflux_mom_u_adv_w_ext=0.0,
3860 dflux_mom_v_adv_p_ext=0.0,
3861 dflux_mom_v_adv_u_ext=0.0,
3862 dflux_mom_v_adv_v_ext=0.0,
3863 dflux_mom_v_adv_w_ext=0.0,
3864 dflux_mom_w_adv_p_ext=0.0,
3865 dflux_mom_w_adv_u_ext=0.0,
3866 dflux_mom_w_adv_v_ext=0.0,
3867 dflux_mom_w_adv_w_ext=0.0,
3872 bc_mom_u_acc_ext=0.0,
3873 bc_dmom_u_acc_u_ext=0.0,
3874 bc_mom_v_acc_ext=0.0,
3875 bc_dmom_v_acc_v_ext=0.0,
3876 bc_mom_w_acc_ext=0.0,
3877 bc_dmom_w_acc_w_ext=0.0,
3878 bc_mass_adv_ext[nSpace],
3879 bc_dmass_adv_u_ext[nSpace],
3880 bc_dmass_adv_v_ext[nSpace],
3881 bc_dmass_adv_w_ext[nSpace],
3882 bc_mom_u_adv_ext[nSpace],
3883 bc_dmom_u_adv_u_ext[nSpace],
3884 bc_dmom_u_adv_v_ext[nSpace],
3885 bc_dmom_u_adv_w_ext[nSpace],
3886 bc_mom_v_adv_ext[nSpace],
3887 bc_dmom_v_adv_u_ext[nSpace],
3888 bc_dmom_v_adv_v_ext[nSpace],
3889 bc_dmom_v_adv_w_ext[nSpace],
3890 bc_mom_w_adv_ext[nSpace],
3891 bc_dmom_w_adv_u_ext[nSpace],
3892 bc_dmom_w_adv_v_ext[nSpace],
3893 bc_dmom_w_adv_w_ext[nSpace],
3894 bc_mom_uu_diff_ten_ext[nSpace],
3895 bc_mom_vv_diff_ten_ext[nSpace],
3896 bc_mom_ww_diff_ten_ext[nSpace],
3897 bc_mom_uv_diff_ten_ext[1],
3898 bc_mom_uw_diff_ten_ext[1],
3899 bc_mom_vu_diff_ten_ext[1],
3900 bc_mom_vw_diff_ten_ext[1],
3901 bc_mom_wu_diff_ten_ext[1],
3902 bc_mom_wv_diff_ten_ext[1],
3903 bc_mom_u_source_ext=0.0,
3904 bc_mom_v_source_ext=0.0,
3905 bc_mom_w_source_ext=0.0,
3906 bc_mom_u_ham_ext=0.0,
3907 bc_dmom_u_ham_grad_p_ext[nSpace],
3908 bc_mom_v_ham_ext=0.0,
3909 bc_dmom_v_ham_grad_p_ext[nSpace],
3910 bc_mom_w_ham_ext=0.0,
3911 bc_dmom_w_ham_grad_p_ext[nSpace],
3912 fluxJacobian_p_p[nDOF_trial_element],
3913 fluxJacobian_p_u[nDOF_trial_element],
3914 fluxJacobian_p_v[nDOF_trial_element],
3915 fluxJacobian_p_w[nDOF_trial_element],
3916 fluxJacobian_u_p[nDOF_trial_element],
3917 fluxJacobian_u_u[nDOF_trial_element],
3918 fluxJacobian_u_v[nDOF_trial_element],
3919 fluxJacobian_u_w[nDOF_trial_element],
3920 fluxJacobian_v_p[nDOF_trial_element],
3921 fluxJacobian_v_u[nDOF_trial_element],
3922 fluxJacobian_v_v[nDOF_trial_element],
3923 fluxJacobian_v_w[nDOF_trial_element],
3924 fluxJacobian_w_p[nDOF_trial_element],
3925 fluxJacobian_w_u[nDOF_trial_element],
3926 fluxJacobian_w_v[nDOF_trial_element],
3927 fluxJacobian_w_w[nDOF_trial_element],
3928 jac_ext[nSpace*nSpace],
3930 jacInv_ext[nSpace*nSpace],
3931 boundaryJac[nSpace*(nSpace-1)],
3932 metricTensor[(nSpace-1)*(nSpace-1)],
3933 metricTensorDetSqrt,
3934 p_grad_trial_trace[nDOF_trial_element*nSpace],
3935 vel_grad_trial_trace[nDOF_trial_element*nSpace],
3937 p_test_dS[nDOF_test_element],
3938 vel_test_dS[nDOF_test_element],
3940 x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
3941 vel_grad_test_dS[nDOF_trial_element*nSpace],
3945 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty;
3946 ck.calculateMapping_elementBoundary(eN,
3952 mesh_trial_trace_ref.data(),
3953 mesh_grad_trial_trace_ref.data(),
3954 boundaryJac_ref.data(),
3960 metricTensorDetSqrt,
3964 ck.calculateMappingVelocity_elementBoundary(eN,
3968 mesh_velocity_dof.data(),
3970 mesh_trial_trace_ref.data(),
3971 xt_ext,yt_ext,zt_ext,
3978 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
3979 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
3980 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
3982 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3983 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
3987 ck.gradTrialFromRef(&p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
3988 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
3990 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_ext);
3991 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
3992 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
3993 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext);
3994 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
3995 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
3996 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
3997 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
3999 for (
int j=0;j<nDOF_trial_element;j++)
4001 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4002 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4003 for (
int I=0;I<nSpace;I++)
4004 vel_grad_test_dS[j*nSpace+I] = vel_grad_trial_trace[j*nSpace+I]*dS;
4009 bc_p_ext = isDOFBoundary_p.data()[ebNE_kb]*ebqe_bc_p_ext.data()[ebNE_kb]+(1-isDOFBoundary_p.data()[ebNE_kb])*p_ext;
4010 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
4011 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
4012 bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext;
4014 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
4025 elementDiameter.data()[eN],
4026 smagorinskyConstant,
4027 turbulenceClosureModel,
4030 ebqe_vf_ext.data()[ebNE_kb],
4031 ebqe_phi_ext.data()[ebNE_kb],
4032 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
4033 ebqe_kappa_phi_ext.data()[ebNE_kb],
4067 mom_uu_diff_ten_ext,
4068 mom_vv_diff_ten_ext,
4069 mom_ww_diff_ten_ext,
4070 mom_uv_diff_ten_ext,
4071 mom_uw_diff_ten_ext,
4072 mom_vu_diff_ten_ext,
4073 mom_vw_diff_ten_ext,
4074 mom_wu_diff_ten_ext,
4075 mom_wv_diff_ten_ext,
4080 dmom_u_ham_grad_p_ext,
4082 dmom_v_ham_grad_p_ext,
4084 dmom_w_ham_grad_p_ext,
4085 ebqe_dragBeam1.data()[ebNE_kb],
4086 ebqe_dragBeam2.data()[ebNE_kb],
4087 ebqe_dragBeam3.data()[ebNE_kb]);
4095 elementDiameter.data()[eN],
4096 smagorinskyConstant,
4097 turbulenceClosureModel,
4100 bc_ebqe_vf_ext.data()[ebNE_kb],
4101 bc_ebqe_phi_ext.data()[ebNE_kb],
4102 &ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],
4103 ebqe_kappa_phi_ext.data()[ebNE_kb],
4116 bc_dmom_u_acc_u_ext,
4118 bc_dmom_v_acc_v_ext,
4120 bc_dmom_w_acc_w_ext,
4126 bc_dmom_u_adv_u_ext,
4127 bc_dmom_u_adv_v_ext,
4128 bc_dmom_u_adv_w_ext,
4130 bc_dmom_v_adv_u_ext,
4131 bc_dmom_v_adv_v_ext,
4132 bc_dmom_v_adv_w_ext,
4134 bc_dmom_w_adv_u_ext,
4135 bc_dmom_w_adv_v_ext,
4136 bc_dmom_w_adv_w_ext,
4137 bc_mom_uu_diff_ten_ext,
4138 bc_mom_vv_diff_ten_ext,
4139 bc_mom_ww_diff_ten_ext,
4140 bc_mom_uv_diff_ten_ext,
4141 bc_mom_uw_diff_ten_ext,
4142 bc_mom_vu_diff_ten_ext,
4143 bc_mom_vw_diff_ten_ext,
4144 bc_mom_wu_diff_ten_ext,
4145 bc_mom_wv_diff_ten_ext,
4146 bc_mom_u_source_ext,
4147 bc_mom_v_source_ext,
4148 bc_mom_w_source_ext,
4150 bc_dmom_u_ham_grad_p_ext,
4152 bc_dmom_v_ham_grad_p_ext,
4154 bc_dmom_w_ham_grad_p_ext,
4155 ebqe_dragBeam1.data()[ebNE_kb],
4156 ebqe_dragBeam2.data()[ebNE_kb],
4157 ebqe_dragBeam3.data()[ebNE_kb]);
4159 if (turbulenceClosureModel >= 3)
4161 const double turb_var_grad_0_dummy[3] = {0.,0.,0.};
4162 const double c_mu = 0.09;
4171 ebqe_vf_ext.data()[ebNE_kb],
4172 ebqe_phi_ext.data()[ebNE_kb],
4175 ebqe_turb_var_0.data()[ebNE_kb],
4176 ebqe_turb_var_1.data()[ebNE_kb],
4177 turb_var_grad_0_dummy,
4178 mom_uu_diff_ten_ext,
4179 mom_vv_diff_ten_ext,
4180 mom_ww_diff_ten_ext,
4181 mom_uv_diff_ten_ext,
4182 mom_uw_diff_ten_ext,
4183 mom_vu_diff_ten_ext,
4184 mom_vw_diff_ten_ext,
4185 mom_wu_diff_ten_ext,
4186 mom_wv_diff_ten_ext,
4199 ebqe_vf_ext.data()[ebNE_kb],
4200 ebqe_phi_ext.data()[ebNE_kb],
4203 ebqe_turb_var_0.data()[ebNE_kb],
4204 ebqe_turb_var_1.data()[ebNE_kb],
4205 turb_var_grad_0_dummy,
4206 bc_mom_uu_diff_ten_ext,
4207 bc_mom_vv_diff_ten_ext,
4208 bc_mom_ww_diff_ten_ext,
4209 bc_mom_uv_diff_ten_ext,
4210 bc_mom_uw_diff_ten_ext,
4211 bc_mom_vu_diff_ten_ext,
4212 bc_mom_vw_diff_ten_ext,
4213 bc_mom_wu_diff_ten_ext,
4214 bc_mom_wv_diff_ten_ext,
4215 bc_mom_u_source_ext,
4216 bc_mom_v_source_ext,
4217 bc_mom_w_source_ext);
4222 mass_adv_ext[0] -= MOVING_DOMAIN*xt_ext;
4223 mass_adv_ext[1] -= MOVING_DOMAIN*yt_ext;
4224 mass_adv_ext[2] -= MOVING_DOMAIN*zt_ext;
4226 mom_u_adv_ext[0] -= MOVING_DOMAIN*mom_u_acc_ext*xt_ext;
4227 mom_u_adv_ext[1] -= MOVING_DOMAIN*mom_u_acc_ext*yt_ext;
4228 mom_u_adv_ext[2] -= MOVING_DOMAIN*mom_u_acc_ext*zt_ext;
4229 dmom_u_adv_u_ext[0] -= MOVING_DOMAIN*dmom_u_acc_u_ext*xt_ext;
4230 dmom_u_adv_u_ext[1] -= MOVING_DOMAIN*dmom_u_acc_u_ext*yt_ext;
4231 dmom_u_adv_u_ext[2] -= MOVING_DOMAIN*dmom_u_acc_u_ext*zt_ext;
4233 mom_v_adv_ext[0] -= MOVING_DOMAIN*mom_v_acc_ext*xt_ext;
4234 mom_v_adv_ext[1] -= MOVING_DOMAIN*mom_v_acc_ext*yt_ext;
4235 mom_v_adv_ext[2] -= MOVING_DOMAIN*mom_v_acc_ext*zt_ext;
4236 dmom_v_adv_v_ext[0] -= MOVING_DOMAIN*dmom_v_acc_v_ext*xt_ext;
4237 dmom_v_adv_v_ext[1] -= MOVING_DOMAIN*dmom_v_acc_v_ext*yt_ext;
4238 dmom_v_adv_v_ext[2] -= MOVING_DOMAIN*dmom_v_acc_v_ext*zt_ext;
4240 mom_w_adv_ext[0] -= MOVING_DOMAIN*mom_w_acc_ext*xt_ext;
4241 mom_w_adv_ext[1] -= MOVING_DOMAIN*mom_w_acc_ext*yt_ext;
4242 mom_w_adv_ext[2] -= MOVING_DOMAIN*mom_w_acc_ext*zt_ext;
4243 dmom_w_adv_w_ext[0] -= MOVING_DOMAIN*dmom_w_acc_w_ext*xt_ext;
4244 dmom_w_adv_w_ext[1] -= MOVING_DOMAIN*dmom_w_acc_w_ext*yt_ext;
4245 dmom_w_adv_w_ext[2] -= MOVING_DOMAIN*dmom_w_acc_w_ext*zt_ext;
4248 bc_mom_u_adv_ext[0] -= MOVING_DOMAIN*bc_mom_u_acc_ext*xt_ext;
4249 bc_mom_u_adv_ext[1] -= MOVING_DOMAIN*bc_mom_u_acc_ext*yt_ext;
4250 bc_mom_u_adv_ext[2] -= MOVING_DOMAIN*bc_mom_u_acc_ext*zt_ext;
4252 bc_mom_v_adv_ext[0] -= MOVING_DOMAIN*bc_mom_v_acc_ext*xt_ext;
4253 bc_mom_v_adv_ext[1] -= MOVING_DOMAIN*bc_mom_v_acc_ext*yt_ext;
4254 bc_mom_v_adv_ext[2] -= MOVING_DOMAIN*bc_mom_v_acc_ext*zt_ext;
4256 bc_mom_w_adv_ext[0] -= MOVING_DOMAIN*bc_mom_w_acc_ext*xt_ext;
4257 bc_mom_w_adv_ext[1] -= MOVING_DOMAIN*bc_mom_w_acc_ext*yt_ext;
4258 bc_mom_w_adv_ext[2] -= MOVING_DOMAIN*bc_mom_w_acc_ext*zt_ext;
4263 isDOFBoundary_u.data()[ebNE_kb],
4264 isDOFBoundary_v.data()[ebNE_kb],
4265 isDOFBoundary_w.data()[ebNE_kb],
4266 isAdvectiveFluxBoundary_p.data()[ebNE_kb],
4267 isAdvectiveFluxBoundary_u.data()[ebNE_kb],
4268 isAdvectiveFluxBoundary_v.data()[ebNE_kb],
4269 isAdvectiveFluxBoundary_w.data()[ebNE_kb],
4270 dmom_u_ham_grad_p_ext[0],
4277 ebqe_bc_flux_mass_ext.data()[ebNE_kb],
4278 ebqe_bc_flux_mom_u_adv_ext.data()[ebNE_kb],
4279 ebqe_bc_flux_mom_v_adv_ext.data()[ebNE_kb],
4280 ebqe_bc_flux_mom_w_adv_ext.data()[ebNE_kb],
4304 dflux_mom_u_adv_p_ext,
4305 dflux_mom_u_adv_u_ext,
4306 dflux_mom_u_adv_v_ext,
4307 dflux_mom_u_adv_w_ext,
4308 dflux_mom_v_adv_p_ext,
4309 dflux_mom_v_adv_u_ext,
4310 dflux_mom_v_adv_v_ext,
4311 dflux_mom_v_adv_w_ext,
4312 dflux_mom_w_adv_p_ext,
4313 dflux_mom_w_adv_u_ext,
4314 dflux_mom_w_adv_v_ext,
4315 dflux_mom_w_adv_w_ext);
4319 ck.calculateGScale(G,normal,h_penalty);
4320 penalty = useMetrics*C_b*h_penalty + (1.0-useMetrics)*ebqe_penalty_ext.data()[ebNE_kb];
4321 for (
int j=0;j<nDOF_trial_element;j++)
4323 int j_nSpace = j*nSpace,ebN_local_kb_j=ebN_local_kb*nDOF_trial_element+j;
4324 fluxJacobian_p_p[j]=0.0;
4325 fluxJacobian_p_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
4326 fluxJacobian_p_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
4327 fluxJacobian_p_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mass_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]);
4329 fluxJacobian_u_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
4330 fluxJacobian_u_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4332 ebqe_phi_ext.data()[ebNE_kb],
4333 sdInfo_u_u_rowptr.data(),
4334 sdInfo_u_u_colind.data(),
4335 isDOFBoundary_u.data()[ebNE_kb],
4336 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4338 mom_uu_diff_ten_ext,
4339 vel_trial_trace_ref.data()[ebN_local_kb_j],
4340 &vel_grad_trial_trace[j_nSpace],
4342 fluxJacobian_u_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4344 ebqe_phi_ext.data()[ebNE_kb],
4345 sdInfo_u_v_rowptr.data(),
4346 sdInfo_u_v_colind.data(),
4347 isDOFBoundary_v.data()[ebNE_kb],
4348 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4350 mom_uv_diff_ten_ext,
4351 vel_trial_trace_ref.data()[ebN_local_kb_j],
4352 &vel_grad_trial_trace[j_nSpace],
4354 fluxJacobian_u_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_u_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4356 ebqe_phi_ext.data()[ebNE_kb],
4357 sdInfo_u_w_rowptr.data(),
4358 sdInfo_u_w_colind.data(),
4359 isDOFBoundary_w.data()[ebNE_kb],
4360 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4362 mom_uw_diff_ten_ext,
4363 vel_trial_trace_ref.data()[ebN_local_kb_j],
4364 &vel_grad_trial_trace[j_nSpace],
4367 fluxJacobian_v_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
4368 fluxJacobian_v_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4370 ebqe_phi_ext.data()[ebNE_kb],
4371 sdInfo_v_u_rowptr.data(),
4372 sdInfo_v_u_colind.data(),
4373 isDOFBoundary_u.data()[ebNE_kb],
4374 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4376 mom_vu_diff_ten_ext,
4377 vel_trial_trace_ref.data()[ebN_local_kb_j],
4378 &vel_grad_trial_trace[j_nSpace],
4380 fluxJacobian_v_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4382 ebqe_phi_ext.data()[ebNE_kb],
4383 sdInfo_v_v_rowptr.data(),
4384 sdInfo_v_v_colind.data(),
4385 isDOFBoundary_v.data()[ebNE_kb],
4386 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4388 mom_vv_diff_ten_ext,
4389 vel_trial_trace_ref.data()[ebN_local_kb_j],
4390 &vel_grad_trial_trace[j_nSpace],
4392 fluxJacobian_v_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_v_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4394 ebqe_phi_ext.data()[ebNE_kb],
4395 sdInfo_v_w_rowptr.data(),
4396 sdInfo_v_w_colind.data(),
4397 isDOFBoundary_w.data()[ebNE_kb],
4398 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4400 mom_vw_diff_ten_ext,
4401 vel_trial_trace_ref.data()[ebN_local_kb_j],
4402 &vel_grad_trial_trace[j_nSpace],
4405 fluxJacobian_w_p[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_p_ext,p_trial_trace_ref.data()[ebN_local_kb_j]);
4406 fluxJacobian_w_u[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_u_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4408 ebqe_phi_ext.data()[ebNE_kb],
4409 sdInfo_w_u_rowptr.data(),
4410 sdInfo_w_u_colind.data(),
4411 isDOFBoundary_u.data()[ebNE_kb],
4412 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4414 mom_wu_diff_ten_ext,
4415 vel_trial_trace_ref.data()[ebN_local_kb_j],
4416 &vel_grad_trial_trace[j_nSpace],
4418 fluxJacobian_w_v[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_v_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4420 ebqe_phi_ext.data()[ebNE_kb],
4421 sdInfo_w_v_rowptr.data(),
4422 sdInfo_w_v_colind.data(),
4423 isDOFBoundary_v.data()[ebNE_kb],
4424 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4426 mom_wv_diff_ten_ext,
4427 vel_trial_trace_ref.data()[ebN_local_kb_j],
4428 &vel_grad_trial_trace[j_nSpace],
4430 fluxJacobian_w_w[j]=
ck.ExteriorNumericalAdvectiveFluxJacobian(dflux_mom_w_adv_w_ext,vel_trial_trace_ref.data()[ebN_local_kb_j]) +
4432 ebqe_phi_ext.data()[ebNE_kb],
4433 sdInfo_w_w_rowptr.data(),
4434 sdInfo_w_w_colind.data(),
4435 isDOFBoundary_w.data()[ebNE_kb],
4436 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4438 mom_ww_diff_ten_ext,
4439 vel_trial_trace_ref.data()[ebN_local_kb_j],
4440 &vel_grad_trial_trace[j_nSpace],
4446 for (
int i=0;i<nDOF_test_element;i++)
4448 int eN_i = eN*nDOF_test_element+i;
4449 for (
int j=0;j<nDOF_trial_element;j++)
4451 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;
4453 globalJacobian.data()[csrRowIndeces_p_p.data()[eN_i] + csrColumnOffsets_eb_p_p.data()[ebN_i_j]] += fluxJacobian_p_p[j]*p_test_dS[i];
4454 globalJacobian.data()[csrRowIndeces_p_u.data()[eN_i] + csrColumnOffsets_eb_p_u.data()[ebN_i_j]] += fluxJacobian_p_u[j]*p_test_dS[i];
4455 globalJacobian.data()[csrRowIndeces_p_v.data()[eN_i] + csrColumnOffsets_eb_p_v.data()[ebN_i_j]] += fluxJacobian_p_v[j]*p_test_dS[i];
4456 globalJacobian.data()[csrRowIndeces_p_w.data()[eN_i] + csrColumnOffsets_eb_p_w.data()[ebN_i_j]] += fluxJacobian_p_w[j]*p_test_dS[i];
4458 globalJacobian.data()[csrRowIndeces_u_p.data()[eN_i] + csrColumnOffsets_eb_u_p.data()[ebN_i_j]] += fluxJacobian_u_p[j]*vel_test_dS[i];
4459 globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*vel_test_dS[i]+
4460 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4461 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4463 vel_trial_trace_ref.data()[ebN_local_kb_j],
4465 sdInfo_u_u_rowptr.data(),
4466 sdInfo_u_u_colind.data(),
4467 mom_uu_diff_ten_ext,
4468 &vel_grad_test_dS[i*nSpace]);
4469 globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_eb_u_v.data()[ebN_i_j]] += fluxJacobian_u_v[j]*vel_test_dS[i]+
4470 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4471 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4473 vel_trial_trace_ref.data()[ebN_local_kb_j],
4475 sdInfo_u_v_rowptr.data(),
4476 sdInfo_u_v_colind.data(),
4477 mom_uv_diff_ten_ext,
4478 &vel_grad_test_dS[i*nSpace]);
4479 globalJacobian.data()[csrRowIndeces_u_w.data()[eN_i] + csrColumnOffsets_eb_u_w.data()[ebN_i_j]] += fluxJacobian_u_w[j]*vel_test_dS[i]+
4480 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb],
4481 isDiffusiveFluxBoundary_u.data()[ebNE_kb],
4483 vel_trial_trace_ref.data()[ebN_local_kb_j],
4485 sdInfo_u_w_rowptr.data(),
4486 sdInfo_u_w_colind.data(),
4487 mom_uw_diff_ten_ext,
4488 &vel_grad_test_dS[i*nSpace]);
4490 globalJacobian.data()[csrRowIndeces_v_p.data()[eN_i] + csrColumnOffsets_eb_v_p.data()[ebN_i_j]] += fluxJacobian_v_p[j]*vel_test_dS[i];
4491 globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_eb_v_u.data()[ebN_i_j]] += fluxJacobian_v_u[j]*vel_test_dS[i]+
4492 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4493 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4495 vel_trial_trace_ref.data()[ebN_local_kb_j],
4497 sdInfo_v_u_rowptr.data(),
4498 sdInfo_v_u_colind.data(),
4499 mom_vu_diff_ten_ext,
4500 &vel_grad_test_dS[i*nSpace]);
4501 globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_eb_v_v.data()[ebN_i_j]] += fluxJacobian_v_v[j]*vel_test_dS[i]+
4502 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4503 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4505 vel_trial_trace_ref.data()[ebN_local_kb_j],
4507 sdInfo_v_v_rowptr.data(),
4508 sdInfo_v_v_colind.data(),
4509 mom_vv_diff_ten_ext,
4510 &vel_grad_test_dS[i*nSpace]);
4511 globalJacobian.data()[csrRowIndeces_v_w.data()[eN_i] + csrColumnOffsets_eb_v_w.data()[ebN_i_j]] += fluxJacobian_v_w[j]*vel_test_dS[i]+
4512 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb],
4513 isDiffusiveFluxBoundary_v.data()[ebNE_kb],
4515 vel_trial_trace_ref.data()[ebN_local_kb_j],
4517 sdInfo_v_w_rowptr.data(),
4518 sdInfo_v_w_colind.data(),
4519 mom_vw_diff_ten_ext,
4520 &vel_grad_test_dS[i*nSpace]);
4522 globalJacobian.data()[csrRowIndeces_w_p.data()[eN_i] + csrColumnOffsets_eb_w_p.data()[ebN_i_j]] += fluxJacobian_w_p[j]*vel_test_dS[i];
4523 globalJacobian.data()[csrRowIndeces_w_u.data()[eN_i] + csrColumnOffsets_eb_w_u.data()[ebN_i_j]] += fluxJacobian_w_u[j]*vel_test_dS[i]+
4524 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_u.data()[ebNE_kb],
4525 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4527 vel_trial_trace_ref.data()[ebN_local_kb_j],
4529 sdInfo_w_u_rowptr.data(),
4530 sdInfo_w_u_colind.data(),
4531 mom_wu_diff_ten_ext,
4532 &vel_grad_test_dS[i*nSpace]);
4533 globalJacobian.data()[csrRowIndeces_w_v.data()[eN_i] + csrColumnOffsets_eb_w_v.data()[ebN_i_j]] += fluxJacobian_w_v[j]*vel_test_dS[i]+
4534 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_v.data()[ebNE_kb],
4535 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4537 vel_trial_trace_ref.data()[ebN_local_kb_j],
4539 sdInfo_w_v_rowptr.data(),
4540 sdInfo_w_v_colind.data(),
4541 mom_wv_diff_ten_ext,
4542 &vel_grad_test_dS[i*nSpace]);
4543 globalJacobian.data()[csrRowIndeces_w_w.data()[eN_i] + csrColumnOffsets_eb_w_w.data()[ebN_i_j]] += fluxJacobian_w_w[j]*vel_test_dS[i]+
4544 ck.ExteriorElementBoundaryDiffusionAdjointJacobian(isDOFBoundary_w.data()[ebNE_kb],
4545 isDiffusiveFluxBoundary_w.data()[ebNE_kb],
4547 vel_trial_trace_ref.data()[ebN_local_kb_j],
4549 sdInfo_w_w_rowptr.data(),
4550 sdInfo_w_w_colind.data(),
4551 mom_ww_diff_ten_ext,
4552 &vel_grad_test_dS[i*nSpace]);
4712 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
4713 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
4714 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
4715 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
4716 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
4717 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
4718 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
4719 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
4720 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
4721 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
4722 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
4723 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
4724 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
4725 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
4726 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
4727 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
4728 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
4729 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
4730 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
4731 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
4732 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
4733 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
4734 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
4735 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
4736 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
4737 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
4738 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
4739 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
4740 double hFactor = args.
scalar<
double>(
"hFactor");
4741 int nElements_global = args.
scalar<
int>(
"nElements_global");
4742 double useRBLES = args.
scalar<
double>(
"useRBLES");
4743 double useMetrics = args.
scalar<
double>(
"useMetrics");
4744 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
4745 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
4746 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
4747 double sigma = args.
scalar<
double>(
"sigma");
4752 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
4753 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
4754 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
4755 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
4756 double C_dc = args.
scalar<
double>(
"C_dc");
4757 double C_b = args.
scalar<
double>(
"C_b");
4758 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
4759 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
4760 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
4761 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
4762 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
4763 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
4764 xt::pyarray<double>& g = args.
array<
double>(
"g");
4765 xt::pyarray<double>& rho_init = args.
array<
double>(
"rho_init");
4766 const double useVF = args.
scalar<
double>(
"useVF");
4767 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
4768 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
4769 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
4770 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
4771 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
4772 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
4773 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
4774 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
4775 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
4776 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
4777 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
4778 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
4779 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
4780 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
4781 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
4782 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
4783 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
4784 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
4785 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
4786 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
4787 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
4788 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
4789 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
4790 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
4791 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
4792 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
4793 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
4794 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
4795 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
4796 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
4797 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
4798 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
4799 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
4800 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
4801 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
4802 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
4803 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
4804 int offset_p = args.
scalar<
int>(
"offset_p");
4805 int offset_u = args.
scalar<
int>(
"offset_u");
4806 int offset_v = args.
scalar<
int>(
"offset_v");
4807 int offset_w = args.
scalar<
int>(
"offset_w");
4808 int stride_p = args.
scalar<
int>(
"stride_p");
4809 int stride_u = args.
scalar<
int>(
"stride_u");
4810 int stride_v = args.
scalar<
int>(
"stride_v");
4811 int stride_w = args.
scalar<
int>(
"stride_w");
4812 xt::pyarray<double>& cg = args.
array<
double>(
"cg");
4813 xt::pyarray<double>& force = args.
array<
double>(
"force");
4814 xt::pyarray<double>& moment = args.
array<
double>(
"moment");
4815 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
4816 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
4817 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
4818 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
4819 xt::pyarray<int>& forceExtractionFaces = args.
array<
int>(
"forceExtractionFaces");
4820 int nForceExtractionFaces = args.
scalar<
int>(
"nForceExtractionFaces");
4821 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
4822 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
4823 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
4824 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
4825 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
4826 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
4827 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
4828 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
4829 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
4830 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
4831 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
4832 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
4833 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
4834 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
4835 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
4836 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
4837 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
4838 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
4839 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
4840 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
4841 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
4842 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
4843 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
4844 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
4845 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
4846 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
4847 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
4848 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
4849 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
4850 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
4851 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
4859 double tmp1[nSpace], tmp2[nSpace*nSpace];
4861 for (
int fef = 0; fef < nForceExtractionFaces; fef++)
4863 int ebNE = forceExtractionFaces.data()[fef];
4866 int ebN = exteriorElementBoundariesArray.data()[ebNE],
4867 eN = elementBoundaryElementsArray.data()[ebN*2+0],
4868 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
4869 eN_nDOF_trial_element = eN*nDOF_trial_element;
4870 double eps_rho,eps_mu;
4874 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
4877 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
4878 ebNE_kb_nSpace = ebNE_kb*nSpace,
4879 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
4880 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
4889 jac_ext[nSpace*nSpace],
4891 jacInv_ext[nSpace*nSpace],
4892 boundaryJac[nSpace*(nSpace-1)],
4893 metricTensor[(nSpace-1)*(nSpace-1)],
4894 metricTensorDetSqrt,
4895 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],vel_grad_test_dS[nDOF_test_element*nSpace],
4896 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
4897 normal[3],x_ext,y_ext,z_ext,
4898 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,gi[nSpace], unormal,gnormal,uneg,
4899 H_rho,d_rho, H_mu,d_mu, rho, mu,nu;
4901 ck.calculateMapping_elementBoundary(eN,
4907 mesh_trial_trace_ref.data(),
4908 mesh_grad_trial_trace_ref.data(),
4909 boundaryJac_ref.data(),
4915 metricTensorDetSqrt,
4919 dS = metricTensorDetSqrt*dS_ref.data()[kb];
4922 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
4923 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
4925 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
4926 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
4928 H_rho = (1.0-useVF)*
RANS2P_IB::smoothedHeaviside(eps_rho,ebqe_phi_ext.data()[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,ebqe_vf_ext.data()[ebNE_kb]));
4930 H_mu = (1.0-useVF)*
RANS2P_IB::smoothedHeaviside(eps_mu,ebqe_phi_ext.data()[ebNE_kb]) + useVF*fmin(1.0,fmax(0.0,ebqe_vf_ext.data()[ebNE_kb]));
4936 ck.gradTrialFromRef(& p_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,p_grad_trial_trace);
4937 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
4939 ck.valFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],&p_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],p_ext);
4940 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
4941 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
4942 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext);
4943 ck.gradFromDOF(p_dof.data(),&p_l2g.data()[eN_nDOF_trial_element],p_grad_trial_trace,grad_p_ext);
4944 ck.gradFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_u_ext);
4945 ck.gradFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_v_ext);
4946 ck.gradFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],vel_grad_trial_trace,grad_w_ext);
4948 for (
int j=0;j<nDOF_trial_element;j++)
4950 int j_nSpace = j*nSpace;
4951 p_test_dS[j] = p_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4952 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
4954 vel_grad_test_dS[j_nSpace+0] = vel_grad_trial_trace[j_nSpace+0]*dS;
4955 vel_grad_test_dS[j_nSpace+1] = vel_grad_trial_trace[j_nSpace+1]*dS;
4956 vel_grad_test_dS[j_nSpace+2] = vel_grad_trial_trace[j_nSpace+2]*dS;
4959 ebqe_velocity.data()[ebNE_kb_nSpace + 0 ] = u_ext;
4960 ebqe_velocity.data()[ebNE_kb_nSpace + 1 ] = v_ext;
4961 ebqe_velocity.data()[ebNE_kb_nSpace + 2 ] = w_ext;
4964 double norm_S,h_e=elementDiameter.data()[eN],t_nu_0,t_nu_1;
4967 switch (turbulenceClosureModel)
4971 norm_S = sqrt(2.0*(grad_u_ext[0]*grad_u_ext[0] + grad_v_ext[1]*grad_v_ext[1] + grad_w_ext[2]*grad_w_ext[2] +
4972 0.5*(grad_u_ext[1]+grad_v_ext[0])*(grad_u_ext[1]+grad_v_ext[0]) +
4973 0.5*(grad_u_ext[2]+grad_w_ext[0])*(grad_u_ext[2]+grad_w_ext[0]) +
4974 0.5*(grad_v_ext[2]+grad_w_ext[1])*(grad_v_ext[2]+grad_w_ext[1])));
4975 t_nu_0 += smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
4976 t_nu_1 += smagorinskyConstant*smagorinskyConstant*h_e*h_e*norm_S;
4980 double re_0,cs_0,re_1,cs_1;
4981 norm_S = sqrt(2.0*(grad_u_ext[0]*grad_u_ext[0] + grad_v_ext[1]*grad_v_ext[1] + grad_w_ext[2]*grad_w_ext[2] +
4982 0.5*(grad_u_ext[1]+grad_v_ext[0])*(grad_u_ext[1]+grad_v_ext[0]) +
4983 0.5*(grad_u_ext[2]+grad_w_ext[0])*(grad_u_ext[2]+grad_w_ext[0]) +
4984 0.5*(grad_v_ext[2]+grad_w_ext[1])*(grad_v_ext[2]+grad_w_ext[1])));
4985 re_0 = h_e*h_e*norm_S/
nu_0;
4986 cs_0=0.027*pow(10.0,-3.23*pow(re_0,-0.92));
4987 t_nu_0 += cs_0*h_e*h_e*norm_S;
4988 re_1 = h_e*h_e*norm_S/
nu_1;
4989 cs_1=0.027*pow(10.0,-3.23*pow(re_1,-0.92));
4990 t_nu_1 += cs_1*h_e*h_e*norm_S;
4995 nu = t_nu_0*(1.0-H_mu)+t_nu_1*H_mu;
4996 mu =
rho_0*t_nu_0*(1.0-H_mu)+
rho_1*t_nu_1*H_mu;
5006 gi[0] = ebqe_bc_u_ext.data()[ebNE_kb];
5007 gi[1] = ebqe_bc_v_ext.data()[ebNE_kb];
5008 gi[2] = ebqe_bc_w_ext.data()[ebNE_kb];
5010 unormal = normal[0]*u_ext + normal[1]*v_ext + normal[2]*w_ext;
5011 gnormal = normal[0]*gi[0] + normal[1]*gi[1] + normal[2]*gi[2];
5012 uneg = 0.5*(unormal - fabs(unormal));
5017 ck.calculateGScale(G,normal,h_penalty);
5018 penalty = C_b*mu/h_penalty;
5023 tmp1[0] = -mu*(grad_u_ext[0]*normal[0] + grad_u_ext[0]*normal[0]+
5024 grad_u_ext[1]*normal[1] + grad_v_ext[0]*normal[1]+
5025 grad_u_ext[2]*normal[2] + grad_w_ext[0]*normal[2])
5026 + penalty*(u_ext-gi[0]) + p_ext*normal[0] - uneg*rho*(u_ext-gi[0]);
5028 tmp1[1] = -mu*(grad_v_ext[0]*normal[0] + grad_u_ext[1]*normal[0]+
5029 grad_v_ext[1]*normal[1] + grad_v_ext[1]*normal[1]+
5030 grad_v_ext[2]*normal[2] + grad_w_ext[1]*normal[2])
5031 + penalty*(v_ext-gi[1]) + p_ext*normal[1] - uneg*rho*(v_ext-gi[1]);
5033 tmp1[2] = -mu*(grad_w_ext[0]*normal[0] + grad_u_ext[2]*normal[0]+
5034 grad_w_ext[1]*normal[1] + grad_v_ext[2]*normal[1]+
5035 grad_w_ext[2]*normal[2] + grad_w_ext[2]*normal[2])
5036 + penalty*(w_ext-gi[2]) + p_ext*normal[2] - uneg*rho*(w_ext-gi[2]);
5038 tmp2[0*nSpace + 0] = 0.0;
5039 tmp2[0*nSpace + 1] = mu*( (v_ext - gi[1])*normal[0] + (u_ext - gi[0])*normal[1] );
5040 tmp2[0*nSpace + 2] = mu*( (w_ext - gi[2])*normal[0] + (u_ext - gi[0])*normal[2] );
5042 tmp2[1*nSpace + 0] = mu*( (u_ext - gi[0])*normal[1] + (v_ext - gi[1])*normal[0] );
5043 tmp2[1*nSpace + 1] = 0.0;
5044 tmp2[1*nSpace + 2] = mu*( (w_ext - gi[2])*normal[1] + (v_ext - gi[1])*normal[2] );
5046 tmp2[2*nSpace + 0] = mu*( (u_ext - gi[0])*normal[2] + (w_ext - gi[2])*normal[0] );
5047 tmp2[2*nSpace + 1] = mu*( (v_ext - gi[1])*normal[2] + (w_ext - gi[2])*normal[1] );
5048 tmp2[2*nSpace + 2] = 0.0;
5053 force.data()[0] += tmp1[0]*dS;
5054 force.data()[1] += tmp1[1]*dS;
5055 force.data()[2] += tmp1[2]*dS;
5060 moment.data()[0] += ( (y_ext-cg.data()[1])*tmp1[2] - (z_ext-cg.data()[2])*tmp1[1] + tmp2[2*nSpace + 1] - tmp2[1*nSpace + 2] )*dS;
5061 moment.data()[1] += ( (z_ext-cg.data()[2])*tmp1[0] - (x_ext-cg.data()[0])*tmp1[2] + tmp2[0*nSpace + 2] - tmp2[2*nSpace + 0] )*dS;
5062 moment.data()[2] += ( (x_ext-cg.data()[0])*tmp1[1] - (y_ext-cg.data()[1])*tmp1[0] + tmp2[1*nSpace + 0] - tmp2[0*nSpace + 1] )*dS;
5072 xt::pyarray<double>& mesh_trial_ref = args.
array<
double>(
"mesh_trial_ref");
5073 xt::pyarray<double>& mesh_grad_trial_ref = args.
array<
double>(
"mesh_grad_trial_ref");
5074 xt::pyarray<double>& mesh_dof = args.
array<
double>(
"mesh_dof");
5075 xt::pyarray<double>& mesh_velocity_dof = args.
array<
double>(
"mesh_velocity_dof");
5076 double MOVING_DOMAIN = args.
scalar<
double>(
"MOVING_DOMAIN");
5077 xt::pyarray<int>& mesh_l2g = args.
array<
int>(
"mesh_l2g");
5078 xt::pyarray<double>& dV_ref = args.
array<
double>(
"dV_ref");
5079 xt::pyarray<double>& p_trial_ref = args.
array<
double>(
"p_trial_ref");
5080 xt::pyarray<double>& p_grad_trial_ref = args.
array<
double>(
"p_grad_trial_ref");
5081 xt::pyarray<double>& p_test_ref = args.
array<
double>(
"p_test_ref");
5082 xt::pyarray<double>& p_grad_test_ref = args.
array<
double>(
"p_grad_test_ref");
5083 xt::pyarray<double>& vel_trial_ref = args.
array<
double>(
"vel_trial_ref");
5084 xt::pyarray<double>& vel_grad_trial_ref = args.
array<
double>(
"vel_grad_trial_ref");
5085 xt::pyarray<double>& vel_test_ref = args.
array<
double>(
"vel_test_ref");
5086 xt::pyarray<double>& vel_grad_test_ref = args.
array<
double>(
"vel_grad_test_ref");
5087 xt::pyarray<double>& mesh_trial_trace_ref = args.
array<
double>(
"mesh_trial_trace_ref");
5088 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.
array<
double>(
"mesh_grad_trial_trace_ref");
5089 xt::pyarray<double>& dS_ref = args.
array<
double>(
"dS_ref");
5090 xt::pyarray<double>& p_trial_trace_ref = args.
array<
double>(
"p_trial_trace_ref");
5091 xt::pyarray<double>& p_grad_trial_trace_ref = args.
array<
double>(
"p_grad_trial_trace_ref");
5092 xt::pyarray<double>& p_test_trace_ref = args.
array<
double>(
"p_test_trace_ref");
5093 xt::pyarray<double>& p_grad_test_trace_ref = args.
array<
double>(
"p_grad_test_trace_ref");
5094 xt::pyarray<double>& vel_trial_trace_ref = args.
array<
double>(
"vel_trial_trace_ref");
5095 xt::pyarray<double>& vel_grad_trial_trace_ref = args.
array<
double>(
"vel_grad_trial_trace_ref");
5096 xt::pyarray<double>& vel_test_trace_ref = args.
array<
double>(
"vel_test_trace_ref");
5097 xt::pyarray<double>& vel_grad_test_trace_ref = args.
array<
double>(
"vel_grad_test_trace_ref");
5098 xt::pyarray<double>& normal_ref = args.
array<
double>(
"normal_ref");
5099 xt::pyarray<double>& boundaryJac_ref = args.
array<
double>(
"boundaryJac_ref");
5100 double eb_adjoint_sigma = args.
scalar<
double>(
"eb_adjoint_sigma");
5101 xt::pyarray<double>& elementDiameter = args.
array<
double>(
"elementDiameter");
5102 xt::pyarray<double>& nodeDiametersArray = args.
array<
double>(
"nodeDiametersArray");
5103 double hFactor = args.
scalar<
double>(
"hFactor");
5104 int nElements_global = args.
scalar<
int>(
"nElements_global");
5105 int nElementBoundaries_owned = args.
scalar<
int>(
"nElementBoundaries_owned");
5106 double useRBLES = args.
scalar<
double>(
"useRBLES");
5107 double useMetrics = args.
scalar<
double>(
"useMetrics");
5108 double alphaBDF = args.
scalar<
double>(
"alphaBDF");
5109 double epsFact_rho = args.
scalar<
double>(
"epsFact_rho");
5110 double epsFact_mu = args.
scalar<
double>(
"epsFact_mu");
5111 double sigma = args.
scalar<
double>(
"sigma");
5116 double smagorinskyConstant = args.
scalar<
double>(
"smagorinskyConstant");
5117 int turbulenceClosureModel = args.
scalar<
int>(
"turbulenceClosureModel");
5118 double Ct_sge = args.
scalar<
double>(
"Ct_sge");
5119 double Cd_sge = args.
scalar<
double>(
"Cd_sge");
5120 double C_dc = args.
scalar<
double>(
"C_dc");
5121 double C_b = args.
scalar<
double>(
"C_b");
5122 const xt::pyarray<double>& eps_solid = args.
array<
double>(
"eps_solid");
5123 const xt::pyarray<double>& phi_solid = args.
array<
double>(
"phi_solid");
5124 const xt::pyarray<double>& q_velocity_solid = args.
array<
double>(
"q_velocity_solid");
5125 const xt::pyarray<double>& q_porosity = args.
array<
double>(
"q_porosity");
5126 const xt::pyarray<double>& q_dragAlpha = args.
array<
double>(
"q_dragAlpha");
5127 const xt::pyarray<double>& q_dragBeta = args.
array<
double>(
"q_dragBeta");
5128 const xt::pyarray<double>& q_mass_source = args.
array<
double>(
"q_mass_source");
5129 const xt::pyarray<double>& q_turb_var_0 = args.
array<
double>(
"q_turb_var_0");
5130 const xt::pyarray<double>& q_turb_var_1 = args.
array<
double>(
"q_turb_var_1");
5131 const xt::pyarray<double>& q_turb_var_grad_0 = args.
array<
double>(
"q_turb_var_grad_0");
5132 xt::pyarray<int>& p_l2g = args.
array<
int>(
"p_l2g");
5133 xt::pyarray<int>& vel_l2g = args.
array<
int>(
"vel_l2g");
5134 xt::pyarray<double>& p_dof = args.
array<
double>(
"p_dof");
5135 xt::pyarray<double>& u_dof = args.
array<
double>(
"u_dof");
5136 xt::pyarray<double>& v_dof = args.
array<
double>(
"v_dof");
5137 xt::pyarray<double>& w_dof = args.
array<
double>(
"w_dof");
5138 xt::pyarray<double>& g = args.
array<
double>(
"g");
5139 const double useVF = args.
scalar<
double>(
"useVF");
5140 xt::pyarray<double>& vf = args.
array<
double>(
"vf");
5141 xt::pyarray<double>&
phi = args.
array<
double>(
"phi");
5142 xt::pyarray<double>& normal_phi = args.
array<
double>(
"normal_phi");
5143 xt::pyarray<double>& kappa_phi = args.
array<
double>(
"kappa_phi");
5144 xt::pyarray<double>& q_mom_u_acc = args.
array<
double>(
"q_mom_u_acc");
5145 xt::pyarray<double>& q_mom_v_acc = args.
array<
double>(
"q_mom_v_acc");
5146 xt::pyarray<double>& q_mom_w_acc = args.
array<
double>(
"q_mom_w_acc");
5147 xt::pyarray<double>& q_mass_adv = args.
array<
double>(
"q_mass_adv");
5148 xt::pyarray<double>& q_mom_u_acc_beta_bdf = args.
array<
double>(
"q_mom_u_acc_beta_bdf");
5149 xt::pyarray<double>& q_mom_v_acc_beta_bdf = args.
array<
double>(
"q_mom_v_acc_beta_bdf");
5150 xt::pyarray<double>& q_mom_w_acc_beta_bdf = args.
array<
double>(
"q_mom_w_acc_beta_bdf");
5151 xt::pyarray<double>& q_velocity_sge = args.
array<
double>(
"q_velocity_sge");
5152 xt::pyarray<double>& q_cfl = args.
array<
double>(
"q_cfl");
5153 xt::pyarray<double>& q_numDiff_u = args.
array<
double>(
"q_numDiff_u");
5154 xt::pyarray<double>& q_numDiff_v = args.
array<
double>(
"q_numDiff_v");
5155 xt::pyarray<double>& q_numDiff_w = args.
array<
double>(
"q_numDiff_w");
5156 xt::pyarray<double>& q_numDiff_u_last = args.
array<
double>(
"q_numDiff_u_last");
5157 xt::pyarray<double>& q_numDiff_v_last = args.
array<
double>(
"q_numDiff_v_last");
5158 xt::pyarray<double>& q_numDiff_w_last = args.
array<
double>(
"q_numDiff_w_last");
5159 xt::pyarray<int>& sdInfo_u_u_rowptr = args.
array<
int>(
"sdInfo_u_u_rowptr");
5160 xt::pyarray<int>& sdInfo_u_u_colind = args.
array<
int>(
"sdInfo_u_u_colind");
5161 xt::pyarray<int>& sdInfo_u_v_rowptr = args.
array<
int>(
"sdInfo_u_v_rowptr");
5162 xt::pyarray<int>& sdInfo_u_v_colind = args.
array<
int>(
"sdInfo_u_v_colind");
5163 xt::pyarray<int>& sdInfo_u_w_rowptr = args.
array<
int>(
"sdInfo_u_w_rowptr");
5164 xt::pyarray<int>& sdInfo_u_w_colind = args.
array<
int>(
"sdInfo_u_w_colind");
5165 xt::pyarray<int>& sdInfo_v_v_rowptr = args.
array<
int>(
"sdInfo_v_v_rowptr");
5166 xt::pyarray<int>& sdInfo_v_v_colind = args.
array<
int>(
"sdInfo_v_v_colind");
5167 xt::pyarray<int>& sdInfo_v_u_rowptr = args.
array<
int>(
"sdInfo_v_u_rowptr");
5168 xt::pyarray<int>& sdInfo_v_u_colind = args.
array<
int>(
"sdInfo_v_u_colind");
5169 xt::pyarray<int>& sdInfo_v_w_rowptr = args.
array<
int>(
"sdInfo_v_w_rowptr");
5170 xt::pyarray<int>& sdInfo_v_w_colind = args.
array<
int>(
"sdInfo_v_w_colind");
5171 xt::pyarray<int>& sdInfo_w_w_rowptr = args.
array<
int>(
"sdInfo_w_w_rowptr");
5172 xt::pyarray<int>& sdInfo_w_w_colind = args.
array<
int>(
"sdInfo_w_w_colind");
5173 xt::pyarray<int>& sdInfo_w_u_rowptr = args.
array<
int>(
"sdInfo_w_u_rowptr");
5174 xt::pyarray<int>& sdInfo_w_u_colind = args.
array<
int>(
"sdInfo_w_u_colind");
5175 xt::pyarray<int>& sdInfo_w_v_rowptr = args.
array<
int>(
"sdInfo_w_v_rowptr");
5176 xt::pyarray<int>& sdInfo_w_v_colind = args.
array<
int>(
"sdInfo_w_v_colind");
5177 int offset_p = args.
scalar<
int>(
"offset_p");
5178 int offset_u = args.
scalar<
int>(
"offset_u");
5179 int offset_v = args.
scalar<
int>(
"offset_v");
5180 int offset_w = args.
scalar<
int>(
"offset_w");
5181 int stride_p = args.
scalar<
int>(
"stride_p");
5182 int stride_u = args.
scalar<
int>(
"stride_u");
5183 int stride_v = args.
scalar<
int>(
"stride_v");
5184 int stride_w = args.
scalar<
int>(
"stride_w");
5185 xt::pyarray<double>& globalResidual = args.
array<
double>(
"globalResidual");
5186 int nExteriorElementBoundaries_global = args.
scalar<
int>(
"nExteriorElementBoundaries_global");
5187 xt::pyarray<int>& exteriorElementBoundariesArray = args.
array<
int>(
"exteriorElementBoundariesArray");
5188 xt::pyarray<int>& elementBoundaryElementsArray = args.
array<
int>(
"elementBoundaryElementsArray");
5189 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.
array<
int>(
"elementBoundaryLocalElementBoundariesArray");
5190 xt::pyarray<double>& ebqe_vf_ext = args.
array<
double>(
"ebqe_vf_ext");
5191 xt::pyarray<double>& bc_ebqe_vf_ext = args.
array<
double>(
"bc_ebqe_vf_ext");
5192 xt::pyarray<double>& ebqe_phi_ext = args.
array<
double>(
"ebqe_phi_ext");
5193 xt::pyarray<double>& bc_ebqe_phi_ext = args.
array<
double>(
"bc_ebqe_phi_ext");
5194 xt::pyarray<double>& ebqe_normal_phi_ext = args.
array<
double>(
"ebqe_normal_phi_ext");
5195 xt::pyarray<double>& ebqe_kappa_phi_ext = args.
array<
double>(
"ebqe_kappa_phi_ext");
5196 const xt::pyarray<double>& ebqe_porosity_ext = args.
array<
double>(
"ebqe_porosity_ext");
5197 const xt::pyarray<double>& ebqe_turb_var_0 = args.
array<
double>(
"ebqe_turb_var_0");
5198 const xt::pyarray<double>& ebqe_turb_var_1 = args.
array<
double>(
"ebqe_turb_var_1");
5199 xt::pyarray<int>& isDOFBoundary_p = args.
array<
int>(
"isDOFBoundary_p");
5200 xt::pyarray<int>& isDOFBoundary_u = args.
array<
int>(
"isDOFBoundary_u");
5201 xt::pyarray<int>& isDOFBoundary_v = args.
array<
int>(
"isDOFBoundary_v");
5202 xt::pyarray<int>& isDOFBoundary_w = args.
array<
int>(
"isDOFBoundary_w");
5203 xt::pyarray<int>& isAdvectiveFluxBoundary_p = args.
array<
int>(
"isAdvectiveFluxBoundary_p");
5204 xt::pyarray<int>& isAdvectiveFluxBoundary_u = args.
array<
int>(
"isAdvectiveFluxBoundary_u");
5205 xt::pyarray<int>& isAdvectiveFluxBoundary_v = args.
array<
int>(
"isAdvectiveFluxBoundary_v");
5206 xt::pyarray<int>& isAdvectiveFluxBoundary_w = args.
array<
int>(
"isAdvectiveFluxBoundary_w");
5207 xt::pyarray<int>& isDiffusiveFluxBoundary_u = args.
array<
int>(
"isDiffusiveFluxBoundary_u");
5208 xt::pyarray<int>& isDiffusiveFluxBoundary_v = args.
array<
int>(
"isDiffusiveFluxBoundary_v");
5209 xt::pyarray<int>& isDiffusiveFluxBoundary_w = args.
array<
int>(
"isDiffusiveFluxBoundary_w");
5210 xt::pyarray<double>& ebqe_bc_p_ext = args.
array<
double>(
"ebqe_bc_p_ext");
5211 xt::pyarray<double>& ebqe_bc_flux_mass_ext = args.
array<
double>(
"ebqe_bc_flux_mass_ext");
5212 xt::pyarray<double>& ebqe_bc_flux_mom_u_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_u_adv_ext");
5213 xt::pyarray<double>& ebqe_bc_flux_mom_v_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_v_adv_ext");
5214 xt::pyarray<double>& ebqe_bc_flux_mom_w_adv_ext = args.
array<
double>(
"ebqe_bc_flux_mom_w_adv_ext");
5215 xt::pyarray<double>& ebqe_bc_u_ext = args.
array<
double>(
"ebqe_bc_u_ext");
5216 xt::pyarray<double>& ebqe_bc_flux_u_diff_ext = args.
array<
double>(
"ebqe_bc_flux_u_diff_ext");
5217 xt::pyarray<double>& ebqe_penalty_ext = args.
array<
double>(
"ebqe_penalty_ext");
5218 xt::pyarray<double>& ebqe_bc_v_ext = args.
array<
double>(
"ebqe_bc_v_ext");
5219 xt::pyarray<double>& ebqe_bc_flux_v_diff_ext = args.
array<
double>(
"ebqe_bc_flux_v_diff_ext");
5220 xt::pyarray<double>& ebqe_bc_w_ext = args.
array<
double>(
"ebqe_bc_w_ext");
5221 xt::pyarray<double>& ebqe_bc_flux_w_diff_ext = args.
array<
double>(
"ebqe_bc_flux_w_diff_ext");
5222 xt::pyarray<double>& q_x = args.
array<
double>(
"q_x");
5223 xt::pyarray<double>& q_velocity = args.
array<
double>(
"q_velocity");
5224 xt::pyarray<double>& ebqe_velocity = args.
array<
double>(
"ebqe_velocity");
5225 xt::pyarray<double>& flux = args.
array<
double>(
"flux");
5226 xt::pyarray<double>& elementResidual_p_save = args.
array<
double>(
"elementResidual_p_save");
5227 xt::pyarray<int>& boundaryFlags = args.
array<
int>(
"boundaryFlags");
5228 xt::pyarray<double>& barycenters = args.
array<
double>(
"barycenters");
5229 xt::pyarray<double>& wettedAreas = args.
array<
double>(
"wettedAreas");
5230 xt::pyarray<double>& netForces_p = args.
array<
double>(
"netForces_p");
5231 xt::pyarray<double>& netForces_v = args.
array<
double>(
"netForces_v");
5232 xt::pyarray<double>& netMoments = args.
array<
double>(
"netMoments");
5233 xt::pyarray<double>& q_dragBeam1 = args.
array<
double>(
"q_dragBeam1");
5234 xt::pyarray<double>& q_dragBeam2 = args.
array<
double>(
"q_dragBeam2");
5235 xt::pyarray<double>& q_dragBeam3 = args.
array<
double>(
"q_dragBeam3");
5236 xt::pyarray<double>& ebqe_dragBeam1 = args.
array<
double>(
"ebqe_dragBeam1");
5237 xt::pyarray<double>& ebqe_dragBeam2 = args.
array<
double>(
"ebqe_dragBeam2");
5238 xt::pyarray<double>& ebqe_dragBeam3 = args.
array<
double>(
"ebqe_dragBeam3");
5239 int nBeams = args.
scalar<
int>(
"nBeams");
5240 int nBeamElements = args.
scalar<
int>(
"nBeamElements");
5241 int beam_quadOrder = args.
scalar<
int>(
"beam_quadOrder");
5242 double beam_Cd = args.
scalar<
double>(
"beam_Cd");
5243 xt::pyarray<double>& beamRadius = args.
array<
double>(
"beamRadius");
5244 xt::pyarray<double>& xq = args.
array<
double>(
"xq");
5245 xt::pyarray<double>& yq = args.
array<
double>(
"yq");
5246 xt::pyarray<double>& zq = args.
array<
double>(
"zq");
5247 xt::pyarray<double>& Beam_h = args.
array<
double>(
"Beam_h");
5248 xt::pyarray<double>& dV_beam = args.
array<
double>(
"dV_beam");
5249 xt::pyarray<double>& q1 = args.
array<
double>(
"q1");
5250 xt::pyarray<double>& q2 = args.
array<
double>(
"q2");
5251 xt::pyarray<double>& q3 = args.
array<
double>(
"q3");
5252 xt::pyarray<double>& vel_avg = args.
array<
double>(
"vel_avg");
5253 xt::pyarray<double>& netBeamDrag = args.
array<
double>(
"netBeamDrag");
5257 double globalConservationError=0.0;
5258 for(
int eN=0;eN<nElements_global;eN++)
5260 double eps_rho,eps_mu;
5264 for(
int k=0;k<nQuadraturePoints_element;k++)
5267 int eN_k = eN*nQuadraturePoints_element+k,
5268 eN_k_nSpace = eN_k*nSpace,
5269 eN_nDOF_trial_element = eN*nDOF_trial_element;
5270 double p=0.0,
u=0.0,
v=0.0,
w=0.0,
5271 grad_p[nSpace],grad_u[nSpace],grad_v[nSpace],grad_w[nSpace],
5279 dmass_adv_u[nSpace],
5280 dmass_adv_v[nSpace],
5281 dmass_adv_w[nSpace],
5283 dmom_u_adv_u[nSpace],
5284 dmom_u_adv_v[nSpace],
5285 dmom_u_adv_w[nSpace],
5287 dmom_v_adv_u[nSpace],
5288 dmom_v_adv_v[nSpace],
5289 dmom_v_adv_w[nSpace],
5291 dmom_w_adv_u[nSpace],
5292 dmom_w_adv_v[nSpace],
5293 dmom_w_adv_w[nSpace],
5294 mom_uu_diff_ten[nSpace],
5295 mom_vv_diff_ten[nSpace],
5296 mom_ww_diff_ten[nSpace],
5307 dmom_u_ham_grad_p[nSpace],
5309 dmom_v_ham_grad_p[nSpace],
5311 dmom_w_ham_grad_p[nSpace],
5322 Lstar_u_p[nDOF_test_element],
5323 Lstar_v_p[nDOF_test_element],
5324 Lstar_w_p[nDOF_test_element],
5325 Lstar_u_u[nDOF_test_element],
5326 Lstar_v_v[nDOF_test_element],
5327 Lstar_w_w[nDOF_test_element],
5328 Lstar_p_u[nDOF_test_element],
5329 Lstar_p_v[nDOF_test_element],
5330 Lstar_p_w[nDOF_test_element],
5335 tau_p=0.0,tau_p0=0.0,tau_p1=0.0,
5336 tau_v=0.0,tau_v0=0.0,tau_v1=0.0,
5339 jacInv[nSpace*nSpace],
5340 p_grad_trial[nDOF_trial_element*nSpace],vel_grad_trial[nDOF_trial_element*nSpace],
5341 p_test_dV[nDOF_trial_element],vel_test_dV[nDOF_trial_element],
5342 p_grad_test_dV[nDOF_test_element*nSpace],vel_grad_test_dV[nDOF_test_element*nSpace],
5348 dmom_u_source[nSpace],
5349 dmom_v_source[nSpace],
5350 dmom_w_source[nSpace],
5352 G[nSpace*nSpace],G_dd_G,tr_G,norm_Rv,h_phi, dmom_adv_star[nSpace],dmom_adv_sge[nSpace];
5354 ck.calculateMapping_element(eN,
5358 mesh_trial_ref.data(),
5359 mesh_grad_trial_ref.data(),
5364 ck.calculateH_element(eN,
5366 nodeDiametersArray.data(),
5368 mesh_trial_ref.data(),
5371 ck.calculateMappingVelocity_element(eN,
5373 mesh_velocity_dof.data(),
5375 mesh_trial_ref.data(),
5380 dV = fabs(jacDet)*dV_ref.data()[k];
5381 ck.calculateG(jacInv,G,G_dd_G,tr_G);
5384 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5385 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5388 ck.gradTrialFromRef(&vel_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,vel_grad_trial);
5390 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
u);
5391 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
v);
5392 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_ref.data()[k*nDOF_trial_element],
w);
5395 for (
int j=0;j<nDOF_trial_element;j++)
5397 vel_test_dV[j] = vel_test_ref.data()[k*nDOF_trial_element+j]*dV;
5400 porosity = q_porosity.data()[eN_k];
5404 q_velocity.data()[eN_k_nSpace+0]=
u;
5405 q_velocity.data()[eN_k_nSpace+1]=
v;
5406 q_velocity.data()[eN_k_nSpace+2]=
w;
5407 q_x.data()[eN_k_nSpace+0]=x;
5408 q_x.data()[eN_k_nSpace+1]=y;
5409 q_x.data()[eN_k_nSpace+2]=
z;
5421 q_mom_u_acc_beta_bdf.data()[eN_k],
5427 q_mom_v_acc_beta_bdf.data()[eN_k],
5433 q_mom_w_acc_beta_bdf.data()[eN_k],
5459 q_dragBeam1.data()[eN_k],
5460 q_dragBeam2.data()[eN_k],
5461 q_dragBeam3.data()[eN_k],
5463 netBeamDrag.data());
5500 for (
int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
5502 int ebN = exteriorElementBoundariesArray.data()[ebNE],
5503 eN = elementBoundaryElementsArray.data()[ebN*2+0],
5504 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
5505 eN_nDOF_trial_element = eN*nDOF_trial_element;
5506 double eps_rho, eps_mu;
5507 for (
int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
5509 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
5510 ebNE_kb_nSpace = ebNE_kb*nSpace,
5511 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
5512 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
5522 dmom_u_acc_u_ext=0.0,
5524 dmom_v_acc_v_ext=0.0,
5526 dmom_w_acc_w_ext=0.0,
5527 mass_adv_ext[nSpace],
5528 dmass_adv_u_ext[nSpace],
5529 dmass_adv_v_ext[nSpace],
5530 dmass_adv_w_ext[nSpace],
5531 mom_u_adv_ext[nSpace],
5532 dmom_u_adv_u_ext[nSpace],
5533 dmom_u_adv_v_ext[nSpace],
5534 dmom_u_adv_w_ext[nSpace],
5535 mom_v_adv_ext[nSpace],
5536 dmom_v_adv_u_ext[nSpace],
5537 dmom_v_adv_v_ext[nSpace],
5538 dmom_v_adv_w_ext[nSpace],
5539 mom_w_adv_ext[nSpace],
5540 dmom_w_adv_u_ext[nSpace],
5541 dmom_w_adv_v_ext[nSpace],
5542 dmom_w_adv_w_ext[nSpace],
5543 mom_uu_diff_ten_ext[nSpace],
5544 mom_vv_diff_ten_ext[nSpace],
5545 mom_ww_diff_ten_ext[nSpace],
5546 mom_uv_diff_ten_ext[1],
5547 mom_uw_diff_ten_ext[1],
5548 mom_vu_diff_ten_ext[1],
5549 mom_vw_diff_ten_ext[1],
5550 mom_wu_diff_ten_ext[1],
5551 mom_wv_diff_ten_ext[1],
5552 mom_u_source_ext=0.0,
5553 mom_v_source_ext=0.0,
5554 mom_w_source_ext=0.0,
5556 dmom_u_ham_grad_p_ext[nSpace],
5558 dmom_v_ham_grad_p_ext[nSpace],
5560 dmom_w_ham_grad_p_ext[nSpace],
5561 dmom_u_adv_p_ext[nSpace],
5562 dmom_v_adv_p_ext[nSpace],
5563 dmom_w_adv_p_ext[nSpace],
5565 flux_mom_u_adv_ext=0.0,
5566 flux_mom_v_adv_ext=0.0,
5567 flux_mom_w_adv_ext=0.0,
5568 flux_mom_uu_diff_ext=0.0,
5569 flux_mom_uv_diff_ext=0.0,
5570 flux_mom_uw_diff_ext=0.0,
5571 flux_mom_vu_diff_ext=0.0,
5572 flux_mom_vv_diff_ext=0.0,
5573 flux_mom_vw_diff_ext=0.0,
5574 flux_mom_wu_diff_ext=0.0,
5575 flux_mom_wv_diff_ext=0.0,
5576 flux_mom_ww_diff_ext=0.0,
5581 bc_mom_u_acc_ext=0.0,
5582 bc_dmom_u_acc_u_ext=0.0,
5583 bc_mom_v_acc_ext=0.0,
5584 bc_dmom_v_acc_v_ext=0.0,
5585 bc_mom_w_acc_ext=0.0,
5586 bc_dmom_w_acc_w_ext=0.0,
5587 bc_mass_adv_ext[nSpace],
5588 bc_dmass_adv_u_ext[nSpace],
5589 bc_dmass_adv_v_ext[nSpace],
5590 bc_dmass_adv_w_ext[nSpace],
5591 bc_mom_u_adv_ext[nSpace],
5592 bc_dmom_u_adv_u_ext[nSpace],
5593 bc_dmom_u_adv_v_ext[nSpace],
5594 bc_dmom_u_adv_w_ext[nSpace],
5595 bc_mom_v_adv_ext[nSpace],
5596 bc_dmom_v_adv_u_ext[nSpace],
5597 bc_dmom_v_adv_v_ext[nSpace],
5598 bc_dmom_v_adv_w_ext[nSpace],
5599 bc_mom_w_adv_ext[nSpace],
5600 bc_dmom_w_adv_u_ext[nSpace],
5601 bc_dmom_w_adv_v_ext[nSpace],
5602 bc_dmom_w_adv_w_ext[nSpace],
5603 bc_mom_uu_diff_ten_ext[nSpace],
5604 bc_mom_vv_diff_ten_ext[nSpace],
5605 bc_mom_ww_diff_ten_ext[nSpace],
5606 bc_mom_uv_diff_ten_ext[1],
5607 bc_mom_uw_diff_ten_ext[1],
5608 bc_mom_vu_diff_ten_ext[1],
5609 bc_mom_vw_diff_ten_ext[1],
5610 bc_mom_wu_diff_ten_ext[1],
5611 bc_mom_wv_diff_ten_ext[1],
5612 bc_mom_u_source_ext=0.0,
5613 bc_mom_v_source_ext=0.0,
5614 bc_mom_w_source_ext=0.0,
5615 bc_mom_u_ham_ext=0.0,
5616 bc_dmom_u_ham_grad_p_ext[nSpace],
5617 bc_mom_v_ham_ext=0.0,
5618 bc_dmom_v_ham_grad_p_ext[nSpace],
5619 bc_mom_w_ham_ext=0.0,
5620 bc_dmom_w_ham_grad_p_ext[nSpace],
5621 jac_ext[nSpace*nSpace],
5623 jacInv_ext[nSpace*nSpace],
5624 boundaryJac[nSpace*(nSpace-1)],
5625 metricTensor[(nSpace-1)*(nSpace-1)],
5626 metricTensorDetSqrt,
5627 dS,p_test_dS[nDOF_test_element],vel_test_dS[nDOF_test_element],
5628 p_grad_trial_trace[nDOF_trial_element*nSpace],vel_grad_trial_trace[nDOF_trial_element*nSpace],
5629 vel_grad_test_dS[nDOF_trial_element*nSpace],
5630 normal[3],x_ext,y_ext,z_ext,xt_ext,yt_ext,zt_ext,integralScaling,
5634 G[nSpace*nSpace],G_dd_G,tr_G,h_phi,h_penalty,penalty,
5635 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, junk[nSpace];
5637 ck.calculateMapping_elementBoundary(eN,
5643 mesh_trial_trace_ref.data(),
5644 mesh_grad_trial_trace_ref.data(),
5645 boundaryJac_ref.data(),
5651 metricTensorDetSqrt,
5655 ck.calculateMappingVelocity_elementBoundary(eN,
5659 mesh_velocity_dof.data(),
5661 mesh_trial_trace_ref.data(),
5662 xt_ext,yt_ext,zt_ext,
5672 dS = ((1.0-MOVING_DOMAIN)*metricTensorDetSqrt + MOVING_DOMAIN*integralScaling)*dS_ref.data()[kb];
5675 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
5676 ck.calculateGScale(G,&ebqe_normal_phi_ext.data()[ebNE_kb_nSpace],h_phi);
5678 eps_rho = epsFact_rho*(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5679 eps_mu = epsFact_mu *(useMetrics*h_phi+(1.0-useMetrics)*elementDiameter.data()[eN]);
5683 ck.gradTrialFromRef(&vel_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,vel_grad_trial_trace);
5686 ck.valFromDOF(u_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
5687 ck.valFromDOF(v_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
5688 ck.valFromDOF(w_dof.data(),&vel_l2g.data()[eN_nDOF_trial_element],&vel_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext);
5691 for (
int j=0;j<nDOF_trial_element;j++)
5693 vel_test_dS[j] = vel_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
5696 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
5697 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
5698 bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext;
5700 porosity_ext = ebqe_porosity_ext.data()[ebNE_kb];
5718 phi.data()[ebNE_kb],
5721 ebqe_dragBeam1.data()[ebNE_kb],
5722 ebqe_dragBeam2.data()[ebNE_kb],
5723 ebqe_dragBeam3.data()[ebNE_kb],