Simulation with TGS can blow up and create NaNs

[ikkah]

Library and Version

v5.4.2-106.1 (latest available)

Operating System

Windows 10

Steps to Trigger Behavior

Run the attached test code in place of SnippetContactReport.cpp in the snippets framework.

Code Snippet to Reproduce Behavior

// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions
// are met:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of NVIDIA CORPORATION nor the names of its
//    contributors may be used to endorse or promote products derived
//    from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
// OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Copyright (c) 2008-2024 NVIDIA Corporation. All rights reserved.
// Copyright (c) 2004-2008 AGEIA Technologies, Inc. All rights reserved.
// Copyright (c) 2001-2004 NovodeX AG. All rights reserved.  

// ****************************************************************************
// This snippet illustrates the use of simple contact reports.
//
// It defines a filter shader function that requests touch reports for 
// all pairs, and a contact callback function that saves the contact points.  
// It configures the scene to use this filter and callback, and prints the 
// number of contact reports each frame. If rendering, it renders each 
// contact as a line whose length and direction are defined by the contact 
// impulse.
// 
// ****************************************************************************

#include <vector>
#include "PxPhysicsAPI.h"
#include "../snippetutils/SnippetUtils.h"
#include "../snippetcommon/SnippetPrint.h"
#include "../snippetcommon/SnippetPVD.h"

using namespace physx;

static PxDefaultAllocator		gAllocator;
static PxDefaultErrorCallback	gErrorCallback;
static PxFoundation* gFoundation = NULL;
static PxPhysics* gPhysics = NULL;
static PxDefaultCpuDispatcher* gDispatcher = NULL;
static PxScene* gScene = NULL;
static PxPvd* gPvd = NULL;

std::vector<PxVec3> gContactPositions;
std::vector<PxVec3> gContactImpulses;

static PxFilterFlags contactReportFilterShader( PxFilterObjectAttributes attributes0, PxFilterData filterData0,
	PxFilterObjectAttributes attributes1, PxFilterData filterData1,
	PxPairFlags& pairFlags, const void* constantBlock, PxU32 constantBlockSize )
{
	PX_UNUSED( attributes0 );
	PX_UNUSED( attributes1 );
	PX_UNUSED( filterData0 );
	PX_UNUSED( filterData1 );
	PX_UNUSED( constantBlockSize );
	PX_UNUSED( constantBlock );

	// all initial and persisting reports for everything, with per-point data
	pairFlags = PxPairFlag::eSOLVE_CONTACT | PxPairFlag::eDETECT_DISCRETE_CONTACT
		| PxPairFlag::eNOTIFY_TOUCH_FOUND
		| PxPairFlag::eNOTIFY_TOUCH_PERSISTS
		| PxPairFlag::eNOTIFY_CONTACT_POINTS;
	return PxFilterFlag::eDEFAULT;
}

class ContactReportCallback : public PxSimulationEventCallback
{
	void onConstraintBreak( PxConstraintInfo* constraints, PxU32 count ) { PX_UNUSED( constraints ); PX_UNUSED( count ); }
	void onWake( PxActor** actors, PxU32 count ) { PX_UNUSED( actors ); PX_UNUSED( count ); }
	void onSleep( PxActor** actors, PxU32 count ) { PX_UNUSED( actors ); PX_UNUSED( count ); }
	void onTrigger( PxTriggerPair* pairs, PxU32 count ) { PX_UNUSED( pairs ); PX_UNUSED( count ); }
	void onAdvance( const PxRigidBody* const*, const PxTransform*, const PxU32 ) {}
	void onContact( const PxContactPairHeader& pairHeader, const PxContactPair* pairs, PxU32 nbPairs )
	{
		PX_UNUSED( ( pairHeader ) );
		std::vector<PxContactPairPoint> contactPoints;

		for( PxU32 i = 0; i < nbPairs; i++ )
		{
			PxU32 contactCount = pairs[ i ].contactCount;
			if( contactCount )
			{
				contactPoints.resize( contactCount );
				pairs[ i ].extractContacts( &contactPoints[ 0 ], contactCount );

				for( PxU32 j = 0; j < contactCount; j++ )
				{
					gContactPositions.push_back( contactPoints[ j ].position );
					gContactImpulses.push_back( contactPoints[ j ].impulse );
				}
			}
		}
	}
};

ContactReportCallback gContactReportCallback;

PxRigidDynamic* rootBody;

void initPhysics( bool /*interactive*/ )
{
	gFoundation = PxCreateFoundation( PX_PHYSICS_VERSION, gAllocator, gErrorCallback );
	gPvd = PxCreatePvd( *gFoundation );
	PxPvdTransport* transport = PxDefaultPvdSocketTransportCreate( PVD_HOST, 5425, 10 );
	gPvd->connect( *transport, PxPvdInstrumentationFlag::eALL );
	gPhysics = PxCreatePhysics( PX_PHYSICS_VERSION, *gFoundation, PxTolerancesScale(), true, gPvd );
	PxInitExtensions( *gPhysics, gPvd );
	PxU32 numCores = SnippetUtils::getNbPhysicalCores();
	gDispatcher = PxDefaultCpuDispatcherCreate( numCores == 0 ? 0 : numCores - 1 );
	PxSceneDesc sceneDesc( gPhysics->getTolerancesScale() );
	sceneDesc.cpuDispatcher = gDispatcher;
	sceneDesc.solverType = physx::PxSolverType::eTGS;
	sceneDesc.filterShader = contactReportFilterShader;
	sceneDesc.simulationEventCallback = &gContactReportCallback;
	gScene = gPhysics->createScene( sceneDesc );

	PxPvdSceneClient* pvdClient = gScene->getScenePvdClient();
	if( pvdClient )
	{
		pvdClient->setScenePvdFlag( PxPvdSceneFlag::eTRANSMIT_CONTACTS, true );
	}

	PxTransform localTm( PxVec3( 0, 0, 0 ) );

	rootBody = gPhysics->createRigidDynamic( localTm );
	rootBody->setMassSpaceInertiaTensor( physx::PxVec3( 1.0f ) );
	rootBody->setLinearVelocity( physx::PxVec3( 100.0f ) );
	rootBody->setAngularVelocity( physx::PxVec3( 100.0f ) );
	gScene->addActor( *rootBody );

	auto cogBody = gPhysics->createRigidDynamic( localTm );
	cogBody->setMassSpaceInertiaTensor( physx::PxVec3( 0.01f ) );
	cogBody->setLinearVelocity( physx::PxVec3( 100.0f ) );
	cogBody->setAngularVelocity( physx::PxVec3( 100.0f ) );
	gScene->addActor( *cogBody );

	auto seg0Body = gPhysics->createRigidDynamic( localTm );
	seg0Body->setMassSpaceInertiaTensor( physx::PxVec3( 0.1f ) );
	seg0Body->setLinearVelocity( physx::PxVec3( 100.0f ) );
	seg0Body->setAngularVelocity( physx::PxVec3( 100.0f ) );
	gScene->addActor( *seg0Body );

	{
		auto joint = PxD6JointCreate( *gPhysics, rootBody, localTm, cogBody, localTm );
		joint->setMotion( physx::PxD6Axis::eX, physx::PxD6Motion::eLOCKED );
		joint->setMotion( physx::PxD6Axis::eY, physx::PxD6Motion::eLOCKED );
		joint->setMotion( physx::PxD6Axis::eZ, physx::PxD6Motion::eLOCKED );
		joint->setMotion( physx::PxD6Axis::eTWIST, physx::PxD6Motion::eLOCKED );
		joint->setMotion( physx::PxD6Axis::eSWING1, physx::PxD6Motion::eLIMITED );
		joint->setMotion( physx::PxD6Axis::eSWING2, physx::PxD6Motion::eLIMITED );
		joint->setSwingLimit( physx::PxJointLimitCone( PxDegToRad( 5.0f ), PxDegToRad( 5.0f ), PxSpring( 30.0f, 10.0f ) ) );
	}

	{
		auto joint = PxD6JointCreate( *gPhysics, cogBody, localTm, seg0Body, localTm );
		joint->setMotion( physx::PxD6Axis::eX, physx::PxD6Motion::eLIMITED );
		joint->setMotion( physx::PxD6Axis::eY, physx::PxD6Motion::eLIMITED );
		joint->setMotion( physx::PxD6Axis::eZ, physx::PxD6Motion::eLIMITED );
		joint->setMotion( physx::PxD6Axis::eTWIST, physx::PxD6Motion::eLIMITED );
		joint->setMotion( physx::PxD6Axis::eSWING1, physx::PxD6Motion::eLIMITED );
		joint->setMotion( physx::PxD6Axis::eSWING2, physx::PxD6Motion::eLIMITED );
		joint->setLinearLimit( physx::PxJointLinearLimit( 0.05f, PxSpring( 20.0f, 5.0f ) ) );
		joint->setTwistLimit( physx::PxJointAngularLimitPair( PxDegToRad( -2.0f ), PxDegToRad( 2.0f ), PxSpring( 30.0f, 10.0f ) ) );
		joint->setSwingLimit( physx::PxJointLimitCone( PxDegToRad( 5.0f ), PxDegToRad( 5.0f ), PxSpring( 30.0f, 10.0f ) ) );
	}
}

void stepPhysics( bool /*interactive*/ )
{
	gContactPositions.clear();
	gContactImpulses.clear();

	static int step = 0;

	printf( "vel: %g - pos %g\n", rootBody->getAngularVelocity().magnitude(), rootBody->getGlobalPose().p.y );

	if( step == 0 )
	{
		rootBody->setRigidBodyFlag( physx::PxRigidBodyFlag::eKINEMATIC, true );
		step++;
	}
	else if( step == 1 )
	{
		rootBody->setGlobalPose( physx::PxTransform( physx::PxIDENTITY() ) );
		step++;
	}
	else if( step == 2 )
	{
		rootBody->setGlobalPose( physx::PxTransform( physx::PxIDENTITY() ) );
		rootBody->setRigidBodyFlag( physx::PxRigidBodyFlag::eKINEMATIC, false );
		step++;
	}

	gScene->simulate( 1.0f / 60.0f );
	gScene->fetchResults( true );
}

void cleanupPhysics( bool /*interactive*/ )
{
	PX_RELEASE( gScene );
	PX_RELEASE( gDispatcher );
	PxCloseExtensions();
	PX_RELEASE( gPhysics );
	if( gPvd )
	{
		PxPvdTransport* transport = gPvd->getTransport();
		PX_RELEASE( gPvd );
		PX_RELEASE( transport );
	}
	PX_RELEASE( gFoundation );

	printf( "SnippetContactReport done.\n" );
}

int snippetMain( int, const char* const* )
{
#ifdef RENDER_SNIPPET
	extern void renderLoop();
	renderLoop();
#else
	initPhysics( false );
	for( PxU32 i = 0; i < 250; i++ )
		stepPhysics( false );
	cleanupPhysics( false );
#endif

	return 0;
}

Expected Behavior

I'd expect the simulation to end up in a stable state without blowing up to inf/NaN.

Actual Behavior

Printing the velocities and positions with the printf in the code gives me this:

vel: 173.205 - pos 0
vel: 0 - pos 0
vel: 0 - pos 0
vel: 220.588 - pos 0.783922
vel: 109.336 - pos 3.7014
vel: 89.0139 - pos 6.59901
vel: 29506.7 - pos 5.62387
vel: 413746 - pos -43.4753
vel: 191324 - pos -240.349
vel: 2.49619e+06 - pos -435.468
vel: 5.64549e+12 - pos 4.68051e+06
vel: -nan(ind) - pos -nan(ind)
vel: -nan(ind) - pos -nan(ind)
vel: -nan(ind) - pos -nan(ind)
vel: -nan(ind) - pos -nan(ind)

Running the code in the debug config also triggers some internal asserts in DyTGSDynamics.cpp.

Commenting out the line that enables physx::PxSolverType::eTGS fixes the issue.

[vreutskyy]

Hi @ikkah! Thanks for reporting.